|
Maximum Likelihood Estimation Examples
This work gives MAPLE replicates of ML-estimation examples from Charles H. Franklin lecture notes . Two examples, for Gaussian and Poisson distributions, are included. Gaussian model has two parameters  and Poisson model has one parameter  . Also included the symbolic example for binomial disribution, Gaussian distribution. Finally use of MAPLE's Statistics package for numerical experimentation is illustrated.
|
| |
|
|
Define Data and Models
|
| |
Define two sets sample points to experiment with
| > |
S1 := [1.364, 0.235,-0.846,-0.285,-1.646];
S2 := [5, 0, 1, 1, 0, 3, 2, 3, 4, 1];
|
|
|
(2.1) |
| > |
plots[listplot](S1,axes=boxed,symbol=solidcircle,title="sample vector #1",gridlines=true);
plots[listplot](S2,axes=boxed,symbol=solidcircle,title="sample vector #2",gridlines=true);
|
Define Poisson and Gaussian distributions (PDF)
| > |
PoissonDist := (x,l) -> exp(-l)*l^x/GAMMA(x+1);
|
|
|
(2.2) |
| > |
GaussianDist := (x,sigma,mu) -> (1/sqrt(2*Pi*sigma^2))*exp(-((x-mu)/sigma)^2/2);
|
|
|
(2.3) |
|
|
Compute Likelihood and Solve for Most Likely Parameters
|
| |
|
Gaussian likelihood for both data sets
|
| |
| > |
LG1 := product( GaussianDist(S1[i],s,m),i=1..nops(S1) );
LG2 := product( GaussianDist(S2[i],s,m),i=1..nops(S2) );
|
|
|
(3.1.1) |
For standart deviation  the likelihood is plotted by ...
| > |
plot( subs(s=1,LG1), m=-1..3,axes=boxed, title="likelihood for Gaussian PDF sigma=1");
|
mean that maximizes the likelihood is found as
| > |
extrema( subs(s=1,LG1), {}, {m}, 'muML'); muML;
|
|
|
(3.1.2) |
Likelihood variation for for both parameters is depicted below, peak of this graph gives the most likely parameters
| > |
plot3d( LG1, m=-4..4, s=0..3,axes=boxed,style=patchcontour,labels=[mu,sigma,likelihood],title="Gaussian PDF likelihood for 1st data set");
plot3d( LG2, m=-4..4, s=0..3,axes=boxed,style=patchcontour,labels=[mu,sigma,likelihood],title="Gaussian PDF likelihood for 2nd data set");
|
Most likely parameters are computed using command "extrema"
| > |
extrema(LG1, {}, {s,m}, 'GaussianParameters1ML'); GaussianParameters1ML;
extrema(LG2, {}, {s,m}, 'GaussianParameters2ML'); extrema(LG2, {}, {s,m}, 'GaussianParameters2ML'); GaussianParameters2ML;;
|
|
|
(3.1.3) |
Select the solution with positive sigma
| > |
MLE1 := select( x -> verify(eval(s,x),0,'greater_than'), GaussianParameters1ML)[1];
MLE2 := select( x -> verify(eval(s,x),0,'greater_than'), GaussianParameters2ML)[1];
|
|
|
(3.1.4) |
plot the likelihood functions and ML estimates of parameters
| > |
infolevel[densityplot]:=1;
|
|
|
(3.1.5) |
| > |
Likelihood1DensityPlot := plots[densityplot]( (LG1), m=-4..4, s=0..3,axes=boxed,labels=[mu,sigma],color=yellow,style=patchnogrid,title="Gaussian PDF likelihood for 1st data set",restricttoranges = true, transparency=0.37):
Likelihood2DensityPlot := plots[densityplot]( (LG2), m=-4..4, s=0..3,axes=boxed,labels=[mu,sigma],color=yellow,style=patchnogrid,title="Gaussian PDF likelihood for 2nd data set",restricttoranges = true,scaletorange=0..5e-9, transparency=0.37):
|
|
plots/densityplot: Adjust to range: 0. .. .310390800181162448e-3
plots/densityplot: Adjust to range: 0. .. .500000000000000010e-8
|
|
| > |
Likelihood1ContourPlot := plots[contourplot]( (LG1), m=-4..4, s=0..3,axes=boxed,labels=[mu,sigma],grid=[56,56]):
Likelihood2ContourPlot := plots[contourplot]( (LG2), m=-4..4, s=0..3,axes=boxed,labels=[mu,sigma],grid=[56,56]):
|
| > |
plots[display]({Likelihood1DensityPlot , Likelihood1ContourPlot , plots[pointplot]( eval([m,s], MLE1), color=red, symbol=solidcircle,symbolsize=30, color=blue,transparency=0. ) }, title="Gaussian PDF likelihood for 1st data set and ML estimate (as big dot)");
plots[display]({Likelihood2DensityPlot , Likelihood2ContourPlot , plots[pointplot]( eval([m,s], MLE2), color=red, symbol=solidcircle,symbolsize=30, color=blue,transparency=0. ) }, title="Gaussian PDF likelihood for 2nd data set and ML estimate (as big dot)" );
|
l
|
|
Poisson likelihood for 2nd data set
|
| |
| > |
#LP1 := product( PoissonDist(S1[i],l),i=1..nops(S1));
LP2 := product( PoissonDist(S2[i],l),i=1..nops(S2));
|
|
|
(3.2.1) |
| > |
extrema(LP2, {}, 'l', 'lML'); evalf(lML);
lML := select( x -> verify( eval(l,x), 0, 'greater_than'), lML)[1];
lML := eval(l,lML);
|
|
|
(3.2.2) |
| > |
plots[display]({
plot( LP2, l=0..7, axes=boxed, labels=[lambda, likelihood],gridlines=true),
plots[pointplot]( [[lML, subs( l=lML,LP2)]], symbol=solidcircle, symbolsize=20,color=blue )
});
plots[display]({
plot( log(LP2), l=0..7, axes=boxed, labels=[lambda, log-likelihood],gridlines=true),
plots[pointplot]( [lML, subs(l=lML, log(LP2))], symbol=solidcircle, symbolsize=20,color=blue )
});
|
|
|
|
Symbolic Results
|
| |
This section gives the symbolic results for Gaussian and Binomial distributions
|
ML for Gaussian Distribution
|
| |
PDF as model and log of model
| > |
model := GaussianDist;
logmodel := (x,sigma,mu) -> log(1/sqrt(2*Pi*sigma^2))-((x-mu)/sigma)^2/2;
|
|
|
(4.1.1) |
Likelihood and log-likelihood
| > |
L := product(model(x[i],sigma,mu),i=1..N);
LL := sum( expand(logmodel(x[i],sigma,mu)),i=1..N);
|
|
|
(4.1.2) |
obtain extremum parameter equations
| > |
equ1 := diff(LL,mu);
equ2 := diff(LL,sigma);
|
|
|
(4.1.3) |
solve for parameters that maximizes likelihood
| > |
pML := solve({equ1,equ2},{sigma,mu});
|
|
|
(4.1.4) |
This result corresponds to common definiton of mean and standard deviation for Normal distribution
 
this derivation is also described in wikipedia under the heading "Estimation of parameters" in relation to ML method.
|
|
ML for Binomial Distribution
|
| |
Likelihood for binomial distribution for N success and M fail experiment, with probability of success being p, is given by
|
|
(4.2.1) |
| > |
MLequation := diff(LB,p);
|
|
|
(4.2.2) |
| > |
pML := solve( MLequation, p);
|
|
|
(4.2.3) |
which is obvious result for frequentiest approach
check that obtained solution is maximizes
| > |
ddLL := diff(log(LB), p$2);
|
|
|
(4.2.4) |
| > |
simplify(subs(p=pML,ddLL));
factor(%);
verify(%,0,'less_than') assuming N>0,M>0;
|
|
|
(4.2.5) |
which indicates that found is a maximizing point.
|
|
|
Using MAPLE's Statistics package for Numerical Experiments
|
| |
|
Gaussian case
|
| |
create Gaussian distribution with mean 3 and sigma 1
| > |
GD := Distribution( Normal(3,1) );
|
|
|
(5.1.1) |
| > |
XG := RandomVariable( GD );
|
|
|
(5.1.2) |
| > |
XGV := Sample( XG, 2250);
|
|
|
(5.1.3) |
| > |
PointPlot(XGV,axes=boxed, title="Gaussian distributed samples");
Statistics[Histogram](XGV,axes=boxed, title="Histogram of Gaussian distributed samples");
|
| > |
LL := expand(log(Likelihood( Normal(mu, sigma), XGV))) assuming sigma>0, mu::real;
|
|
|
(5.1.4) |
| > |
#extrema(LL, {}, {mu, sigma}, 'GML');
|
| > |
equs := convert(VectorCalculus[Gradient](LL, [mu, sigma]), 'list');
|
|
|
(5.1.5) |
| > |
solve(equs, {mu, sigma});
|
|
|
(5.1.6) |
| > |
Statistics[Mean]( XGV );
|
|
|
(5.1.7) |
Result indicates that ML estimates of distribution parameters (2.99, 1.001) are close to actual (3,1) and as theoretical results indicates  is equal to mean of sample data
|
|
Poisson case
|
| |
create Gaussian distribution with mean 3 and sigma 1
| > |
PD := Distribution( Poisson(3) );
|
|
|
(5.2.1) |
| > |
XP := RandomVariable( PD );
|
|
|
(5.2.2) |
| > |
XPV := Sample( XP, 2250);
|
|
|
(5.2.3) |
| > |
PointPlot(XPV,axes=boxed, title="Poisson distributed samples");
Statistics[Histogram](XPV,axes=boxed, title="Histogram of Gaussian distributed samples");
|
| > |
LL := expand(log(Likelihood( Poisson(lambda), XPV))) assuming sigma>0, mu::real;
|
|
|
(5.2.4) |
| > |
equs := convert(VectorCalculus[Gradient](LL, [lambda]), 'list');
|
|
|
(5.2.5) |
|
|
(5.2.6) |
|
|
(5.2.7) |
Result indicates that ML estimates of distribution parameters (3.031) are close to actual (3) and as theoretical result indicates that  is equal to mean of sample data
|
|
|
Notes
|
| |
I encountered following problems in this worksheet
|
about simplification
|
| |
failed to simplify log of likelihood function in symbolic form . A number of assumptions are used
| > |
#LL := log(L) assuming N::integer,sigma::real,mu::real;
#LL := sum( expand(log(model(x[i],sigma,mu))),i=1..N);
#LL := expand(LL) assuming N::integer,sigma::real,mu::real;
|
|
|
|
