% LMcs130301_4th.m
% Case Study 13.3.1 p 790-791 in 
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% Written by Eugene.Gallagher@umb.edu 11/14/10, revised 1/21/11, 2/28/11
% Revised: 1/21/11
% X(:,1) Hemoglobin before 60 km walk and X(:,2) after 60-km walk
X=[14.6 13.8
    17.3 15.4
    10.9 11.3
    12.8 11.6
    16.6 16.4
    12.2 12.6
    11.2 11.8
    15.4 15.0       
    14.8 14.4
    16.2 15.0];
D=(X(:,1)-X(:,2));    
[H,P,CI,STATS] = TTEST(X(:,1),X(:,2),0.05,'both');
fprintf(...
 'The mean ''After-Before'' difference in hemoglobin was %5.3f g/dl \n',...
 mean(D));
fprintf(...
 'The sd for the difference in hemoglobin was %5.3f g/dl \n',...
 STATS.sd);
fprintf(...
   'The CI for the difference in hemoglobin (g/dl) was [%5.3f %5.3f]\n',...
 CI);
fprintf('The paired t statistic was %6.4f with %2.0f df\n',...
    STATS.tstat,STATS.df)
fprintf('The paired t test 2-tailed p=%6.4f\n',P);

% Plot histogram and check for symmetry
binsize=0.25;  % Needed in order to properly scale the normal pdf
edges=-3:binsize:3;

% hist(D);
[N,BIN] = histc(D,edges);
bar(edges,N,'histc')
axis([-3.25 3.25 0 2.1])
set(get(gca,'Children'),'FaceColor',[.8 .8 1]);
xlabel('''After-Before'' Difference in Hemoglobin (g/dl)','FontSize',20);
ylabel('Number of Cases','FontSize',20);
ax1=gca;
set(ax1,'Ytick',[0:3])
figure(gcf);pause

% Superimpose the normal probability pdf on a histogram of differences.
% 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);
n=length(D);
mu_j=mean(D);
sigma_j=std(D); % sigmaj is the standard deviation; = 1 after Z transform
y_j=-3.25:0.1:3.25;
fy_j=1/(sqrt(2*pi)*sigma_j)*exp(-1/2*((y_j-mu_j)./sigma_j).^2);
fy_j=n*binsize*fy_j;  % will properly scale the height of the pdf
% 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('Case Study 13.3.1, %1.0f samples untransformed',n);
title(s,'FontSize',22)
figure(gcf);pause