% LMcs140301_4th.m
% Case Study 14.3.1 from 
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% page 815. A case study using Wilcoxon signed rank test
% Written by Eugene.Gallagher@umb.edu 11/16/10 Revised 11/16/10, 3/1/11
% 
D=[13.32 13.06 14.02 11.86 13.58 13.77 13.51 14.42 14.44 15.43];
expected=14.6;
[P,H,STATS] = signrank(D,expected,'alpha',0.05,'method','exact');
fprintf(...
'\nThe signed rank test statistic is %4.1f with exact 2-tailed p=%6.4g\n',...
 STATS.signedrank,P);
[P,H,STATS] = signrank(D,expected,'alpha',0.05,'method','approximate');
fprintf('The signed rank test approximate 2-tailed p=%6.4g for z=%6.4f\n',P,STATS.zval);
[H,P,CI,STATS] = ttest(D,expected,0.05,'both');
fprintf('The one-sample t test 2-tailed p=%6.4g\n',P);
fprintf(...
'The mean TL/HDl = %4.2f with 95%% CI: [%4.2f %4.2f]\n',...
    mean(D), CI(1),CI(2));
[p,h,stats] = signtest(D,expected,0.05,'method','exact');
fprintf('The sign test statistic is %4.1f with exact 2-tailed p=%6.4g\n',stats.sign,p);
[p,h,stats] = signtest(D,expected,'method','approximate');
fprintf('The sign test z=%5.3f with approximate 2-tailed p=%6.4g\n',stats.zval,p);

% Plot histogram and check for symmetry
binsize=0.25;  % Needed in order to properly scale the normal pdf
edges=11:binsize:16;
% hist(D);
[N,BIN] = histc(D,edges);
bar(edges,N,'histc')
axis([10.9 16.1 0 max(N)+0.2])
set(get(gca,'Children'),'FaceColor',[.8 .8 1]);
xlabel('TL/HDl','FontSize',20);
ylabel('Number of Cases','FontSize',20);
ax1=gca;
set(ax1,'Ytick',[0:max(N)],'Xtick',edges)
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=10.9:0.01:16.1;
fy_j=1/(sqrt(2*pi)*sigma_j)*exp(-1/2*((y_j-mu_j)./sigma_j).^2);
fy_j=n*binsize*fy_j;
cutoff=expected;
fy_cutoff=n*binsize/(sqrt(2*pi)*sigma_j)*exp(-1/2*((cutoff-mu_j)./sigma_j).^2);
% 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.
v=axis;
h1=line(y_j,fy_j ,'Color','r','Parent',ax1);
set(h1,'linestyle','--','color','r','linewidth',2)
h1=line([cutoff cutoff]',[0 fy_cutoff]','Color','b','Parent',ax1);
set(h1,'linestyle','-.','color','b','linewidth',3)
h2=line([cutoff cutoff]',[fy_cutoff v(4)]','Color','b','Parent',ax1);
set(h2,'linestyle','-.','color','b','linewidth',3)
s=sprintf('Case Study 14.3.1, %2.0f samples untransformed',n);
title(s,'FontSize',22)
figure(gcf);pause
