% LMcs140203_4th.m
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% page 807. 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=[15 13;12 8;12 12.5;14 12;13 12;13 12.5;13 12.5;12 14;12.5 12;12 11;
    12.5 10];
[p,h,stats] = signtest(D(:,2),D(:,1),0.05,'method','exact');
fprintf('\nThe sign test statistic is %4.1f with exact 1-tailed p=%6.4g\n',stats.sign,p/2);
[p,h,stats] = signtest(D(:,2),D(:,1),'method','approximate');
fprintf('The sign test z=%5.3f with approximate 1-tailed p=%6.4g\n',stats.zval,p/2);
% It is a 1-tailed p test to the left since the expectation is that mean
% circulation time is reduced by 4 months of cyclandelate.
[H,P,CI,STATS] = ttest(D(:,2),D(:,1),0.05,'left');
fprintf('The t statistic was %5.3f with %2.0f df.\n',STATS.tstat,STATS.df)
fprintf('The paired t test 1-tailed p=%6.4g\n',P);
fprintf(...
'The mean circulation time = %4.2f secs with 95%% CI: [%4.2f %4.2f]\n',...
    mean(D(:,2)-D(:,1)), CI(1),CI(2));
[P,H,STATS] = signrank(D(:,2),D(:,1),'alpha',0.05,'method','exact');
fprintf(...
'The sign rank test statistic is %4.1f with exact 1-tailed p=%6.4g\n',...
 STATS.signedrank,P/2);
[P,H,STATS] = signrank(D(:,2),D(:,1),'alpha',0.05,'method','approximate');
fprintf('The sign rank test approximate 1-tailed p=%6.4g for z=%6.4f\n',P/2,STATS.zval);

% Plot histogram and check for symmetry
binsize=.5;  % Needed in order to properly scale the normal pdf
edges=-4:binsize:2.5;
% hist(D);
[N,BIN] = histc(D(:,2)-D(:,1),edges);
bar(edges,N,'histc')
axis([-4.1 2.6 0 3.3])
set(get(gca,'Children'),'FaceColor',[.8 .8 1]);
xlabel('''After-Before'' Mean Circulation Time (secs)','FontSize',20);
ylabel('Number of Cases','FontSize',20);
ax1=gca;
set(ax1,'Ytick',[0:3],'Xtick',[-4:0.5:2.5])
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,c]=size(D);
mu_j=mean(D(:,2)-D(:,1));
sigma_j=std(D(:,2)-D(:,1)); % sigmaj is the standard deviation; = 1 after Z transform
y_j=-4.1:0.01:2.6;
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=0;
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.2.2, %2.0f samples untransformed',n);
title(s,'FontSize',22)
figure(gcf);pause