% LMex040305_4th.m
% Example 4.3.5 Memphis earthquakes and the Poisson distribution
% Pp 305-306 in
% Larsen & Marx (2006) Introduction to Mathematical Statistics, 4th edition
% Written Fall 2010 by Eugene Gallagher for EEOS601 Revised: 1/8/11
% Eugene.Gallagher@umb.edu 
% http://alpha.es.umb.edu/faculty/edg/files/edgwebp.html
% Annual earthquakes registering 2.5 on the Richter scale and within 40
% miles of downtown Memphis have a Poisson distribution with lambda=6.5
% Calculate the exact probability that nine or more such earthquakes will
% strike next year using the Poisson distribution and compare that with the
% normal approximation based on the central limit theorem.
lambda=6.5;
N=9;
P=1-poisscdf(N-1,lambda);
fprintf('\nThe exact p value using the Poisson cumulative distribution ')
fprintf('function with lambda=%5.3f is %7.4f.\n',lambda,P)
P2=1-normcdf(((N-1)-lambda)/sqrt(lambda));
fprintf('\nThe approximate p value using the central limit theorem is ')
fprintf('%7.4f.\n',P2)
continuity=0.5;
% "The continuity correction is appropriate whenever a discrete probability
% model is being approximated by the area under a curve.)
P3=1-normcdf(((N+continuity-1)-lambda)/sqrt(lambda));
fprintf('\nThe approximate p value using the central limit theorem with ');
fprintf('the continuity correction is ')
fprintf('%7.4f.\n',P3)

% Optional: plot the Poisson pdf and the normal approximation
% Use LMFig040301_4th.m as a model
X=0:18;
Y = poisspdf(X,lambda);
bar(X,Y,1);axis([-0.6 19 0 0.2]);title('Example 4.3.5','FontSize',14);
set(get(gca,'Children'),'FaceColor',[.8 .8 1])
ylabel('Probability','FontSize',14);
ax1=gca;set(ax1,'xtick',0:2:18,'FontSize',12);
set(ax1,'ytick',0:0.05:0.2,'FontSize',12)
figure(gcf);pause
% Superimpose the normal probability pdf:
% 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):
% f(y_j)=1/(sqrt(2*pi)*sigma_j)*exp(-1/2*((y_j-mu_j)/sigma_j)^2)
mu_j=lambda;
sigma_j=sqrt(lambda); % sigmaj is the standard deviation
y_j=0:0.2:18;
fy_j=1/(sqrt(2*pi)*sigma_j)*exp(-1/2*((y_j-mu_j)./sigma_j).^2);
% 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)
xlabel('Number of Memphis Earthquakes > 2.5'),
figure(gcf);pause
