% LMEx050301_4th.m
% Simulation of Confidence intervals; 50 CIs each of size n = 4 were drawn
% from the normal pdf
% From Larsen & Marx  (2006) Introduction to Mathematical Statistics,
% Fourth Edition. page 364
% Dept. Environmental, Earth & Ocean Sciences
% Written by Eugene.Gallagher@umb.edu 10/29/10, revised 1/24/11, 3/4/11
% http://alpha.es.umb.edu/faculty/edg/files/edgwebp.html
% Call the inverse of the cumulative normal distribution
% for a 95% CI, 0.025 of the tail should be to the right, so find
% the Z value for 97.5% of the cumulative normal distribution, or 1.96
Z=norminv(0.975);
X=[6.5 9.2 9.9 12.4];meanX=mean(X);sigma=0.8;n=length(X);
CI=[meanX-Z*sigma/sqrt(n);meanX+Z*sigma/sqrt(n)]
fprintf('The mean of X is %4.2f and the 95%% CI is [%4.2f %4.2f]\n',...
    meanX,CI(1),CI(2));
 % Plot Figure 5.3.1
mu=0;
sigma=1;
X=-3.5:0.1:3.5;
Y = normpdf(X,mu,sigma);
plot(X,Y,'--r','LineWidth',2);
axis([-3.55 3.55 0 0.41]);
title('Figure 5.3.1','FontSize',20);
ax1=gca;
xlabel('Z','FontSize',16),
ylabel('f_z(z)','FontSize',16);
ax1=gca;
set(ax1,'XTick',[-1.96 0 1.96] ,'FontSize',16,...
    'XTickLabel',{'-1.96','0','1.96'},'FontSize',18)
set(ax1,'ytick',0:.1:0.4,'FontSize',16)
hold on;
xf=-1.96:.01:1.96;yf=normpdf(xf,mu,sigma);
fill([-1.96 xf 1.96]',[0 yf 0]',[.8 .8 1])
text(1.6,0.25,'Area = 0.95','FontSize',20)
title('Figure 5.3.1','FontSize',22)
figure(gcf);pause
hold off;
% Figure 5.3.2
% Generate n random samples of size 4 with mean 10 and sigma=8
% Larsen & Marx use n=50;
MU=10;SIGMA=0.8;         % Don't change unless plot axes are changed too.
fprintf('mu =%5.3f and sigma = %5.3f\n',MU,SIGMA);
n=100;                    % Larsen & Marx used 50
samsize=4;               % Larsen & Marx used 4
R = normrnd(MU,SIGMA,samsize,n);
[r,c]=size(R);
mR=mean(R)';
CIs=[mR+Z*SIGMA/sqrt(r) mR-Z*SIGMA/sqrt(r)];
% first row of CIs is the upper and 2nd row is the lower CI for 50
% random samples of size 4;
i=find(~(CIs(:,1)>=10 & CIs(:,2)<=10));
j=find(CIs(:,1)>=10 & CIs(:,2)<=10);
fprintf('With sample size %d, of %2.0f CI''s, \n',samsize,n)
fprintf('%2.0f (%5.3f%%) didn''t contain mu=10.\n',length(i),...
    100*length(i)/n)
if ~isempty(i) & n<1000
    disp(CIs(i,:))
end
% display 5.3.2
X=[1:n]';Y=repmat(MU,n,1);
plot(X,Y,'--b');hold on
errorbar(X(i),mR(i),CIs(i,1)-mR(i),CIs(i,2)-mR(i),'-r');
errorbar(X(j),mR(j),CIs(j,1)-mR(j),CIs(j,2)-mR(j),'-g');
s=sprintf('Possible 95%% confidence intervals, sample size=%d',samsize);
xlabel(s,'FontSize',18),
ylabel('True mean','FontSize',18);title('Figure 5.3.2','FontSize',22)
if samsize < 4
    axis([0-0.02*n 1.02*n 6.9 13.1])
    set(gca,'Ytick',7:13,'FontSize',18)
elseif samsize<10
    axis([0-0.02*n 1.02*n 7.9 12.1])
    set(gca,'Ytick',8:12,'FontSize',18)
else
    axis([0-0.02*n 1.02*n 8.4 11.6])
    set(gca,'Ytick',8.5:.5:11.5,'FontSize',18)
end
figure(gcf);pause;hold off