% LMcs140201_4th.m
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% page 804. A case study solved by the sign test
% Written by Eugene.Gallagher@umb.edu 11/16/10 Revised 11/16/10, 3/1/11
% 
D=[7.02 7.35 7.32 7.33 7.15 7.26 7.25 7.35 7.38 7.20 7.31 7.24 7.34 ...
    7.32 7.34 7.14 7.20 7.41 7.77 7.12 7.45 7.28 7.34 7.22 7.32 7.4 ...
    6.99 7.1 7.3 7.21 7.33 7.28 7.35 7.24 7.36 7.09 7.32 6.95 7.35 ...
    7.36 6.6 7.29 7.31];
[p,h,stats] = signtest(D,7.39,0.05,'method','exact')
fprintf('\nThe sign test exact p=%6.4g\n',p);
[p,h,stats] = signtest(D,7.39,'method','approximate');
fprintf('The sign test approximate p=%6.4g;z=%6.4f\n',p,stats.zval);
[H,P,CI,STATS] = ttest(D,7.39);
fprintf('The one-sample t test 2-tailed p=%6.4g\n',P);
fprintf('The mean pH = %4.2f with 95%% CI: [%4.2f %4.2f]\n',mean(D),...
    CI(1),CI(2));
[P,H,STATS] = signrank(D,7.39,'alpha',0.05,'method','exact')
fprintf('The sign rank test exact p=%6.4g\n',P);
[P,H,STATS] = signrank(D,7.39,'alpha',0.05,'method','approximate');
fprintf('The sign rank test approximate p=%6.4g\n',P);
% Plot histogram and check for symmetry
binsize=.1;  % Needed in order to properly scale the normal pdf
edges=6.4:binsize:8;
% hist(D);
[N,BIN] = histc(D,edges);
bar(edges,N,'histc')
axis([6.35 8.05 0 21])
set(get(gca,'Children'),'FaceColor',[.8 .8 1]);
xlabel('pH','FontSize',20);
ylabel('Number of Cases','FontSize',20);
ax1=gca;
set(ax1,'Ytick',[0:5:25],'Xtick',[6.4:0.2:8])
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=6.35:0.01:8.05;
fy_j=normpdf(y_j,mu_j,sigma_j)*n*binsize;
cutoff=7.39;
fy_cutoff=n*binsize*normpdf(cutoff,mu_j,sigma_j);
% will properly scale the height of the pdf 
% Plot using ax1 handle, saved above,to save this graph
% on top of the previous graph.
v=axis;
h1=line(y_j,fy_j ,'Color','r','Parent',ax1);
set(h1,'linestyle','--','color','r','linewidth',3)
h2=line([cutoff cutoff]',[0 fy_cutoff]','Color','b','Parent',ax1);
set(h2,'linestyle','-.','color','r','linewidth',3)
h3=line([cutoff cutoff]',[fy_cutoff v(4)]','Color','b','Parent',ax1);
set(h3,'linestyle','-.','color','b','linewidth',3)
s=sprintf('Case Study 14.2.1, %2.0f samples untransformed',n);
title(s,'FontSize',22)
figure(gcf);pause
