% LMEx040302_4th.m
% Larsen & Marx Example 4.3.2; Example of the central limit theorem
% Page 302-303
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% Integration to find the variance of a uniform pdf using the symbolic math
% toolbox
% Written by Eugene.Gallagher@umb.edu as LMex040302_3rd.m in 2001
% revised 1/8/11
% See also LMEx040304_4th.m  on rounding for variance of uniform pdf
% See also LMEx040304_3rd.m  on rounding for variance of uniform pdf
% Problems in program. The histogram generates a line extending to the
% right of the axis boundary, despite the axis command. This can be edited
% out in a graphics editor.

% Generate 40 sets of 5 uniform random numbers as in Table 4.3.2, p 302
% Figure in Table 4.3.2 uses 7 bins and 40 samples;
% Try 1000 samples and numbins=20 for a prettier plot
 n=5;samples=40; numbins=7; % Numbins for the histogram, multiple of 7
%  n=30;samples=40; numbins=7; % Numbins for the histogram, multiple of 7
% n=5;samples=1000; numbins=14; % Numbins for the histogram, multiple of 7
% n=100;samples=1000; numbins=14; % Numbins for the histogram, multiple of 7
% n=1000;samples=500; numbins=14; % For the histogram
% n=1000;samples=40; numbins=7; % For the histogram
C=rand(samples,n);
C6=sum(C')';  % sums of n replicates
% Use symbolic math toolbox to calculate the variance of a uniform
% random distribution
mu=0.5;
syms f z y
disp('Definite integral (0,1) of y^2')
int(y^2,0,1)
EYsq=int(sym(y^2),0,1);
VARY=EYsq-mu^2;
disp('Var of a uniformly distributed random number on the interval (0,1):')
disp(VARY)
vary=eval(VARY); % convert from symbolic to numeric
% The variance is 1/12
% The Z ratios of the sums are in column C7 of Table 4.3.2
C7=(C6-n*mu)./(sqrt(n)*sqrt(vary));
Table040302=[[1:samples]' C C6 C7];
if n==5   % Reproduce Table 4.3.2
 fprintf('                     Table 4.3.2\n')
 fprintf('      y1      y2      y3      y4      y5      y    Z-ratio\n');
 fprintf('%2.0f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f %7.4f\n',Table040302');
end

% Use the code on page 5-115 of Statistics toolbox to draw histogram:
binwidthscaling=7/numbins; % This will produce 7 bins (3.5-(-3.5))
edges=-3.5:binwidthscaling:3.5;N=histc(C7,edges);
axis([-3.6 3.6 -0.005 0.53])
bar(edges,N/sum(N)/binwidthscaling,1,'histc');
ax1=gca;
set(get(gca,'Children'),'FaceColor',[.8 .8 1])
axis([-3.6 3.6 -0.005 0.53])
ylabel('Density','FontSize',16);xlabel('Z-ratio','FontSize',16)
set(ax1,'xtick',edges,'FontSize',14);
set(ax1,'ytick',0:0.1:0.5,'FontSize',14)
figure(gcf);pause
% Superimpose the normal probability pdf:
% The normal probability equation is provided on p. 293
% This is for mean 0, and unit standard
% deviation. The more general equation (Legendre & Legendre, 1998 p. 147) is:
% f(y_j)=1/(sqrt(2*pi)*sigma_j)*exp(-1/2*((y_j-mu_j)/sigma_j)^2)
mu_j=0;
sigma_j=1; % sigmaj is the standard deviation; = 1 after Z transform
y_j=-3.5:0.1:3.5;
fy_j=1/(sqrt(2*pi)*sigma_j)*exp(-1/2*((y_j-mu_j)./sigma_j).^2);
% fyj=1/(sqrt(2*pi)*sigmaj)*exp(-1/2*((y-muj)/sigmaj).^2); 
% Plot using ax1 handle, saved above,to save this graph
% on top of the previous graph.
h1=line(y_j,fy_j ,'Color','r','Parent',ax1,'Linewidth',2);
set(h1,'linestyle','--','color','r','linewidth',2)
s=sprintf('Example 4.3.2: n=%4.0f and samples=%3.0f.',n,samples);
title(s,'FontSize',20)
figure(gcf);pause