5 % LMcs050202 
% An example of pdf fitting.
% Based on Gallagher's LMEx050202.m
% From Larsen & Marx (2006). Introduction to Mathematical Statistics,
% Fourth Edition. page 359-362
% Note that LM use this example to demonstrate method of moments
% Dept. Environmental, Earth & Ocean Sciences
% Written by Eugene.Gallagher@umb.edu 1/1/2010
% 
% Fit the Gamma pdf to maximum 24-h precipitation levels for 36 storms from
% 1900-1969.
%
hold off; clf
RAIN=[31 2.82 3.98 4.02 9.5 4.5 11.4 10.71 6.31 4.95 5.64 5.51 13.4 9.72 ...
    6.47 10.16 4.21 11.6 4.75 6.85 6.25 3.42 11.8 0.8 3.69 3.1 22.22 ...
    7.43 5.0 4.58 4.46 8 3.73 3.5 6.2 0.67];
[PARMHAT,PARMCI] = gamfit(RAIN)
fprintf(...
'The MLE estimate of r is %5.3f with 95%% CI: [%5.3f %5.3f]\n',...
PARMHAT(1),PARMCI(1,1),PARMCI(2,1))
fprintf('Larsen & Marx find r=1.60, so barely within the 95%% CI\n')
fprintf('Matlab''s gamma b is Larsen & Marx''s 1/lambda\n')
fprintf(...
'The MLE estimate of lambda is %5.3f with 95%% CI: [%5.3f %5.3f]\n',...
1/PARMHAT(2),1/PARMCI(2,2),1/PARMCI(1,2))
fprintf('Larsen & Marx find lambda=0.22, so within the 95%% CI\n')

Binsize=4; % maximum of 16;
X=0:32;y=gampdf(X,PARMHAT(1),PARMHAT(2));plot(X,y,'--b');figure(gcf)
% Histogram
EDGES=0:Binsize:32;N = histc(RAIN,EDGES);
bar(EDGES,N./(sum(N)*Binsize),'histc')
set(get(gca,'Children'),'FaceColor',[.8 .8 1])
;title('Figure 5.2.3','FontSize',22)
xlabel('Maximum 24-hour rainfall (in.)','FontSize',20); 
ylabel('Density','FontSize',20);
figure(gcf);
hold on
x=0:0.1:32;
y=gampdf(x,PARMHAT(1),PARMHAT(2));
% Need to multiply the gamma pdf by binsize to get frequencies, but not
% here
plot(x,y,'--r','Linewidth',3)
figure(gcf);pause
hold off
% generate a random rainfall sequence
RandRain=gamrnd(PARMHAT(1),PARMHAT(2),36,1);
Binsize=2;
EDGES=0:Binsize:ceil(max(RandRain));N = histc(RandRain,EDGES);
bar(EDGES,N./(sum(N)*Binsize),'histc');
set(get(gca,'Children'),'FaceColor',[.8 .8 1])
xlabel('Maximum 24-hour rainfall (in.)','FontSize',20);
ylabel('Density','FontSize',20);
title('Randomly generated rainfall','FontSize',22); figure(gcf);
hold on
x=0:0.1:ceil(max(RandRain));
y=gampdf(x,PARMHAT(1),PARMHAT(2));
% Need to multiply the gamma pdf by binsize
plot(x,y,'--r','Linewidth',3)
figure(gcf);pause
hold off

% Why are Matlab's fits so much different from Larsen & Marx?
y=0:0.1:32;
fyy = gampdf(y,PARMHAT(1),PARMHAT(2));
LMY=gampdf(y,1.60,1/0.22);
fX=gampdf(RAIN,PARMHAT(1),PARMHAT(2));
plot(y,fyy,'LineWidth',2);
xlabel('Maximum 24-hr rainfall (in.)','FontSize',20)
ylabel('Density','FontSize',20)
axis([0 32 0 0.12])
ax1=gca;
hold on
h1=stem(RAIN,fX,':k','filled');
set(h1,'LineWidth',2);
s=sprintf('Case Study 5.2.2, r=%6.4f, lambda=%6.4f, Likelihood=%7.3g',...
    PARMHAT(1),1/PARMHAT(2),prod(fX));
title(s,'FontSize',18)
figure(gcf);pause
hold off
% Why are Matlab's fits so much different from Larsen & Marx?
plot(y,LMY,'LineWidth',2);
xlabel('Maximum 24-hr rainfall (in.)','FontSize',20)
ylabel('Density','FontSize',20)
fX=gampdf(RAIN,1.6,1/0.22);
axis([0 32 0 0.12])
ax1=gca;
hold on
h1=stem(RAIN,fX,':k','filled');
set(h1,'LineWidth',2);
s=sprintf('Case Study 5.2.2, r=%6.4f, lambda=%6.4f, Likelihood=%7.3g',...
    1.6,.22,prod(fX));
title(s,'FontSize',18)
figure(gcf);pause
plot(y,fyy,'--r','LineWidth',3);
hold off
