:

## Simple examples of Maximum Likelihood Estimation

Maple

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.

 > restart;

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

 > LB := p^N*(1-p)^M;
 (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

 > with(Statistics):

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)
 > solve(equs, {lambda});
 (5.2.6)
 > Statistics[Mean](XPV);
 (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

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;

 about graphics failed to overlay a solid graphics on top of density plot, so i use tranparency for density plot but solid graphics seems to be effected.