% LMEx050202
% An example of customized pdf fitting.
% From Larsen & Marx (2006). Introduction to Mathematical Statistics,
% Fourth Edition. page 350
% Dept. Environmental, Earth & Ocean Sciences
% Written by Eugene.Gallagher@umb.edu 10/29/10
% 
% Fit the custom pdf using a Matlab anonymous function @
%
X=[9.2 5.6 18.4 12.1 10.7];
fyytheta=@(y,theta)(1./theta.^2.*y.*exp(-y./theta));
theta=10; % This is an initial guess; 
[theta, thetaCI]= mle(X,'pdf',fyytheta,'start',theta,'lowerbound',0)
% Note that must specify the lower bound for theta is 0 as specified
% in the problem or negative pdf's are produced
% Here is a plot of the pdf for the maximum likelihood estimator.
ezplot('1./5.6.^2.*y.*exp(-y./5.6)',[0.001 100])
figure(gcf);pause

y=0.001:.01:101;
fyy=1./theta.^2.*y.*exp(-y./theta);
fX=1./theta.^2.*X.*exp(-X./theta);
plot(y,fyy,'LineWidth',2);
axis([0 51 -0.002 0.07])
ax1=gca;
hold on
h1=stem(X,fX,':k','filled');
set(h1,'LineWidth',2);
s=sprintf('Example 5.2.2, theta=%3.1f, Likelihood=%7.3g',theta,prod(fX));
title(s,'FontSize',22)
figure(gcf);pause
hold off

% Plot a different value of theta
thetan=theta*2;
y=0.001:.01:101;
fyy=1./thetan.^2.*y.*exp(-y./thetan);
X=[9.2 5.6 18.4 12.1 10.7];
fX=1./thetan.^2.*X.*exp(-X./thetan);
plot(y,fyy,'LineWidth',2);
axis([0 101 -0.002 0.04])
ax1=gca;
hold on
h1=stem(X,fX,':k','filled');
set(h1,'LineWidth',2);
s=sprintf('Example 5.2.2, theta=%3.1f, Likelihood=%7.3g',thetan,prod(fX));
title(s,'FontSize',22)
figure(gcf);pause
hold off
% Plot a different value of theta
thetan=theta/2;
y=0.001:.01:26;
fyy=1./thetan.^2.*y.*exp(-y./thetan);
X=[9.2 5.6 18.4 12.1 10.7];
fX=1./thetan.^2.*X.*exp(-X./thetan);
plot(y,fyy,'LineWidth',2);
axis([0 26 -0.002 0.14])
ax1=gca;
hold on
h1=stem(X,fX,':k','filled');
set(h1,'LineWidth',2);
s=sprintf('Example 5.2.2, theta=%3.1f, Likelihood=%7.3g',thetan,prod(fX));
title(s,'FontSize',22)
figure(gcf);pause
hold off
%---------------------------
% What would happen if an additional observation of 30 were observed?
X=[X 30]; 
[theta, thetaCI]= mle(X,'pdf',fyytheta,'start',theta,'lowerbound',0)
% Note that must specify the lower bound for theta is 0 as specified
% in the problem or negative pdf's are produced
y=0.001:.01:101;
fyy=1./theta.^2.*y.*exp(-y./theta);
fX=1./theta.^2.*X.*exp(-X./theta);
plot(y,fyy,'LineWidth',2);
axis([0 51 -0.002 0.07])
ax1=gca;
hold on
h1=stem(X,fX,':k','filled');
set(h1,'LineWidth',2);
s=sprintf('Example 5.2.2, theta=%3.1f, Likelihood=%7.3g',theta,prod(fX));
title(s,'FontSize',22)
figure(gcf);pause
hold off
% Plot old theta of 5.6
thetan=5.6;
fyy=1./thetan.^2.*y.*exp(-y./thetan);
fX=1./thetan.^2.*X.*exp(-X./thetan);
plot(y,fyy,'LineWidth',2);
axis([0 51 -0.002 0.07])
ax1=gca;
hold on
h1=stem(X,fX,':k','filled');
set(h1,'LineWidth',2);
s=sprintf('Example 5.2.2, theta=%3.1f, Likelihood=%7.3g',thetan,prod(fX));
title(s,'FontSize',22)
figure(gcf);pause
hold off
% Plot a different value of theta
thetan=theta*2;
y=0.001:.01:101;
fyy=1./thetan.^2.*y.*exp(-y./thetan);
fX=1./thetan.^2.*X.*exp(-X./thetan);
plot(y,fyy,'LineWidth',2);
axis([0 101 -0.002 0.04])
ax1=gca;
hold on
h1=stem(X,fX,':k','filled');
set(h1,'LineWidth',2);
s=sprintf('Example 5.2.2, theta=%3.1f, Likelihood=%7.3g',thetan,prod(fX));
title(s,'FontSize',22)
figure(gcf);pause
hold off
% Plot a different value of theta
thetan=theta/2;
y=0.001:.01:51;
fyy=1./thetan.^2.*y.*exp(-y./thetan);
fX=1./thetan.^2.*X.*exp(-X./thetan);
plot(y,fyy,'LineWidth',2);
axis([0 51 -0.002 0.12])
ax1=gca;
hold on
h1=stem(X,fX,':k','filled');
set(h1,'LineWidth',2);
s=sprintf('Example 5.2.2, theta=%3.1f, Likelihood=%7.3g',thetan,prod(fX));
title(s,'FontSize',22)
figure(gcf);pause
hold off
