% LMcs100402_4th.m
% Pp. 618-621 in 
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% Written by Eugene.Gallagher@umb.edu 11/21/10, revised 3/15/11
DATA=[zeros(162,1);ones(267,1);repmat(2,271,1);repmat(3,185,1)
    repmat(4,111,1);repmat(5,61,1);repmat(6,27,1);repmat(7,8,1)
    repmat(8,3,1);repmat(9,1,1)];
% These expanded data needed only for poissfit.m
Deaths=[0:10]';
obs=[162 267 271 185 111 61 27 8 3 1 0]'; n=sum(obs);
% this will fit the Poisson parameter:
lambda=poissfit(DATA)
% Here's a simpler way to calculate the Poisson parameter
lambda=sum(Deaths.*obs)/n
expc=poisspdf(Deaths,lambda)'*n;
% Call Matlab's goodness-of-fit test (see help file for Poisson fit)
[h,p,st] = chi2gof(Deaths,'ctrs',Deaths,'frequency',obs, ...
           'expected',expc,'nparams',1)
% or solve it explicitly
EXP=expc;
EXP(8)=n-sum(expc(1:7));
EXP(9:end)=[];
OBS=obs;
OBS(8)=n-sum(obs(1:7));
OBS(9:end)=[];
d1=sum((OBS'-EXP).^2./EXP);
df=length(OBS)-1-1;
P = 1-chi2cdf(d1,df);
fprintf(...
   'The chi-square statistic,%6.1f, with %2.0f d.f. has pvalue=%6.4f\n',...
   d1,df,P)
alpha=0.05;
criticalval=chi2inv(1-alpha,df);
fprintf('The critical value for alpha=%4.2f is %5.3f\n',alpha,criticalval)
if d1 >=       criticalval
    disp('Reject null')
else
    disp('Fail to reject null')
end
% Plot the results:
bar(0:7,[st.O;st.E]',1,'grouped')
axis([-.75 7.75 0 310]); 
set(gca,'ytick',0:100:300,...
    'xticklabel',{'0 ';'1 ';'2 ';'3 ';'4 ';'5 ';'6';'7+'},'FontSize',20);
legend('Observed','Poisson Expectation','NorthEast','FontSize',20)
xlabel('Number of Deaths','FontSize',20)
ylabel('Frequency','FontSize',20)
title('Case Study 10.4.2','FontSize',22);figure(gcf);