% LMEx050201_4th.m
% An example of binomial MLE fitting & a difference in the way
% Mathworks & Larsen & Marx define the geometric distribution
% From Larsen & Marx (2006). Introduction to Mathematical Statistics,
% Fourth Edition. page 348-349
% Written by Eugene.Gallagher@umb.edu 10/28/10; revised 3/3/11
% Tom Lane on 10/29/10: "unfortunately it looks like your text and MATLAB
% use different definitions for the [geometric] distribution. Our version
% has positive probability on 0,1,2,.... Yours starts at 1. The version we
% use is the one described in "Univariate Discrete Distributions" by 
% Johnson, Kotz, and Kemp. Wikipedia shows both versions.

X=[3 2 1 3]; [Phat,PCI] =mle(X,'distribution','geometric')
% Matlab gives the wrong answer, because it uses a different definition 
% of the geometric distribution, defined for k=0, 1, 2, 3 ... inf
% Larsen & Marx use a form defined for positive k: k=1, 2, 3, ... inf
% Larsen and Marx define the geometric distribution as the number of
% trials before a success, so k=0 is not part of the domain, but Mathworks
% defines the geometric distribution as the number of failures
% before a success, allowing k=0 to be defined.
% The correct MLE for L & M is 4/9
% 10/29 email from Tom Lane, Mathworks, says that I should call
% mle using

[Phat,PCI] =mle(X-1,'distribution','geometric');
format rat;
fprintf('The maximum likelihood estimate for p is\n')
disp(Phat)
format
fprintf('with 95%% CIs: [%5.3f %5.3f]\n',PCI);
% Just following along the text's derivation (p. 348) of the MLE:
syms p; s=diff(5*log(1-p)+4*log(p),p)
 solve(s,p)
 % find the second derivative and plot;
 s2=diff(diff(5*log(1-p)+4*log(p),p))
 % This plot shows that the second derivative is negative for all 0<p<1
 ezplot(s2,[0 1]);figure(gcf);pause
% From Larsen & Marx (2006) p. 349 provides the MLE formula
n=length(X);Phat=n/sum(X)
