Maximum Likelihood Estimation

UPDATED VERSION for MMA 6.0 Copyright 2008, Crooke, Froeb & Tschantz.  This may be used for educational purposes free of charge.

ORIGINAL VERSION at http://library.wolfram.com/infocenter/Articles/1689/
Phillip S. Crooke, Luke Froeb, Steven T. Tschantz 
Organization: Vanderbilt University
Department: Owen Graduate School of Management

Journal: Mathematica in Education and Research
Year: 1999; Volume: 8; Issue: 1; Page range: 17-23

Philip (SuperDagger[Crooke]), Luke (Froeb^‡), and Steven (SuperDagger[Tschantz])
Department of (SuperDagger[Mathematics])
Owen Graduate School of (Management^‡)
Vanderbilt University
Nashville, Tennessee 37240

Introduction

Parameter estimation from data is a problem that is common to many disciplines.  The objective is to infer the value of a parameter or a set of parameters from observed data.  In this paper we implement, using Mathematica, an efficient method for parameter estimation called maximum likelihood.

Consider a random sample (Subscript[X, 1]),(Subscript[X, 2]),…,(Subscript[X, n]) from a probability distribution having the probability density function f(x;θ) where θ is a parameter that lies in some set Ω⊂ℝ (called the parameter space).  The joint probability density function of (Subscript[X, 1]),(Subscript[X, 2]),…,(Subscript[X, n]) is f((Subscript[x, 1]);θ)f((Subscript[x, 2]);θ)⋯ f((Subscript[x, n]);θ).  Regarding the joint probability density function as a function of θ, we define the likelihood function L of the random sample as

	L(θ;(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n]))=f((Subscript[x, 1]);θ)f((Subscript[x, 2]);θ)⋯f((Subscript[x, n]);θ)=(∏_(i=1)^n)f((Subscript[x, i]),θ).

Suppose that there exists a function g((Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n])) such that if (OverHat[θ])=g((Subscript[X, 1]),(Subscript[X, 2]),…,(Subscript[X, n])), then 

	 L((OverHat[θ]);(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n]))=L(g((Subscript[X, 1]),(Subscript[X, 2]),…,(Subscript[X, n]));(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n]))=(max_(θ∈Ω))L(θ;(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n])).
	 
The function (statistic) g((Subscript[X, 1]),(Subscript[X, 2]),…,(Subscript[X, n])) is called the maximum likelihood estimator of θ.  This function may or may not be unique.  If the likelihood function is a differentiable function of θ, then the maximum likelihood estimator can often be obtained by examining the critical points of  L(θ;(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n]))  i.e., ((d L)/(d θ))=0.  If the probability density function depends on k parameters i.e. f(x;(Subscript[θ, 1]),(Subscript[θ, 2]),…,(Subscript[θ, k])), then the likelihood function is defined in an analogous way as

	L((Subscript[θ, 1]),(Subscript[θ, 2]),…,(Subscript[θ, k]);(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n]))=(∏_(i=1)^n)f((Subscript[x, i]);(Subscript[θ, 1]),(Subscript[θ, 2]),…,(Subscript[θ, k]))
	
and maximum likelihood estimators of the parameters	 (Subscript[θ, 1]),(Subscript[θ, 2]),…,(Subscript[θ, k]) are functions (OverHat[(Subscript[θ, i])])=(Subscript[g, i])((Subscript[X, 1]),(Subscript[X, 2]),…,(Subscript[X, n])),i=1,2,…,k such that 

	L((Subscript[(OverHat[θ]), 1]),(Subscript[(OverHat[θ]), 2]),…(Subscript[(OverHat[θ]), k]);(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n]))=(max_(((Subscript[θ, 1]),(Subscript[θ, 2]),…,(Subscript[θ, k]))∈Ω⊂(ℝ^k)))L((Subscript[θ, 1]),(Subscript[θ, 2]),…,(Subscript[θ, k]);(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n])).

Again if L is sufficiently smooth, then estimators are obtained by solving the system of k equations:

		((∂ L((Subscript[θ, 1]),(Subscript[θ, 2]),…,(Subscript[θ, k]);(Subscript[x, 1]),(Subscript[x, 2]),…,(Subscript[x, n])))/(∂ (Subscript[θ, i])))=0, i=1,2,…,k.
		
Without formally defining the concepts, maximum likelihood estimators are consistent (approach the true parameters as the sample size increases), asymptotically efficient (there is no other consistent estimator with a smaller variance), and have a normal limiting distribution, making it possible to construct confidence intervals for the parameters.

Mathematica's analytic and numeric capabilities make it an attractive platform for teaching and doing maximum likelihood estimation.  The programs presented in this paper can handle both discrete and continuous random variables.

The MaximumLikelihood Package (execute these cells only once)

<<Histograms`

<<MultivariateStatistics`

The Mathematica program to calculate the maximum likelihood estimator uses some of the standard add-on packages.

We have combined several functions that compute and illustrate the maximum likelihood algorithm into a package that is called MaximumLikelihood.  The primary function in the package is called ML.  Other functions of the package  analyze and graph the estimates.  The ML function computes parameter estimates that maximize the logarithim of the likelihood function.  This is done using the FindMinimum  function of Mathematica. The program also computes the variance of the estimators and presents the results of the estimation in a standard format.

© 1998, Crooke, Froeb & Tschantz.  This program may be freely distributed, as long as the copyright notice is not changed.  It may not be sold.  If you modify or extend the package with new examples, please send a copy to luke.m.froeb@vanderbilt.edu.

BeginPackage["MaximumLikelihood`"];

ML::usage=
"ML[data,f,startingvalues] returns the maximum likelihood estimates and their estimated variance for the probability density f and corresponding data.  The probability density is a pure function that is mapped over the data.";

LogLikelihood::usage=
"LogLikelihood is the negative of the log likelihood at the values of the estimated parameters, returned by the ML function.";

vcov::usage=
"vcov is the variance-covariance matrix of the parameter estimator returned by the ML function.";

PlotDensities::usage=
"PlotDensities[data,f,g] plots the theoretical (f) and estimated densities (g) for a continuous distribution.";

ListPlotDensities::usage=
"ListPlotDensities[data,f,g] plots the theoretical (f) and estimated densities (g) for a discrete distribution.";

Plot2DConfidenceRegion::usage=
"Plot2DConfidenceRegion[estimate] gives a two-dimensional contour plot of the 99%, 95%, and 90% confidence regions centered on the maximum likelihood estimates.";

ML::nodata="Empty data list in ML routine.";
ML::failure="Failure to find maximum of likelihood in ML routine.";

Begin["`Private`"];
ML[data_List,f_,startingvalues_List]:=
	Module[{L,parameters,h,SE},
	If[Length[data]==0,Message[ML::nodata];Return[];];
	L=Apply[Plus,Log[Map[f,data]]];
	results=Apply[FindMinimum,FlattenAt[{-L,startingvalues,MaxIterations->500},{2}]];		
	parameters=Transpose[startingvalues][[1]];
	If[Not[Apply[And,Map[NumberQ,parameters/.results[[2]]]]],Message[ML::failure];Return[]];
	h=Transpose[(Map[D[Log[Map[f,data]],#]&,parameters])/.results[[2]]];
	vcovresults=Inverse[Transpose[h].h];
	SE=Table[Sqrt[vcovresults[[i,i]] ],{i,Length[parameters]}];		
	Print["-----Maximum Likelihood with Gradient Variance Computation-----"];	
	Print["Number of Observations: "<>ToString[Length[data]] ];	
	Print["Log(L): "<>ToString[-results[[1]] ] ];	
	Print[NumberForm[TableForm[Transpose[{parameters, SE, parameters/SE}]/.results[[2]],TableHeadings->{ToString[#]& /@ parameters,{"Parameter","Standard Error","z-statistic"}}
				],4]];
	Flatten[{LogLikelihood->-results[[1]],results[[2]],vcov->vcovresults}]
	];

PlotDensities[data_List,f_,g_]:=
Show[{
Histograms`Histogram[data, Histograms`HistogramScale -> 1],
Plot[{f[t],g[t]},{t,Min[data],Max[data]},
Axes->{True,False},PlotStyle->{{Thickness[.005]},			
						{Thickness[.005],Dashing[{.02}]}}]
}];

Plot2DConfidenceRegion[estimate_]:=
	Module[{parameters},
	parameters=Drop[Drop[Map[#[[1]]&,estimate],1],-1];
Graphics[{
{PointSize[0.04],RGBColor[1,0,0],Point[parameters]},
MultivariateStatistics`EllipsoidQuantile[MultivariateStatistics`MultinormalDistribution[parameters,vcov],.99],MultivariateStatistics`EllipsoidQuantile[MultivariateStatistics`MultinormalDistribution[parameters,vcov],.95],MultivariateStatistics`EllipsoidQuantile[MultivariateStatistics`MultinormalDistribution[parameters,vcov],.9]}/.estimate,Axes->True]];

ListPlotDensities[data_List,f_,g_]:=
	Module[{},
Show[
Histograms`Histogram[data, 
Histograms`HistogramCategories->Range[Min[data],Max[data]+1],Histograms`HistogramScale -> 1],
ListPlot[
Map[{#,f[#]}&,Range[Min[data],Max[data]]],
Joined->True],
ListPlot[
Map[{#,g[#]}&,Range[Min[data],Max[data]]],
Joined->True,
PlotStyle->Dashing[{.03}]]
]];

End[];
EndPackage[];

In the following sections we demonstrate the use of this package on data generated by some well-known probability distributions.

Example: Normal Distribution

In this section, we illustrate use of the MaximumLikelihood package on determining the mean μ  and standard deviation σ  from a random sample that has been generated by the Normal distribution:

			f(x;μ,σ)=((e^(-((((x+μ))^2)/(2(σ^2)))))/(σ(Sqrt[2π]))), x∈ℝ.
			
In the ContinuousDistribution package,  this probability distribution density is defined as NormalDistribution.

Our first approach to estimating μ and σ is to construct the likelihood function L, and then use Mathematica's built-in function FindMinimum to find a maximum of the likelihood function on the given data.  We will then perform the calculation using ML.  For demonstration purposes, we draw a set of pseudo random numbers from the Normal distribution with μ=10 and σ=1.

T=30;
exact={μ->10,σ->1};
data=Table[Random[NormalDistribution[μ,σ]/.exact],{T}];

Next, we define the probability density that we would like to fit as a pure function.

Clear[f];
f=PDF[NormalDistribution[μ,σ],#]&;
f[x]

The maximum likelihood estimates are are choosen to maximize the logarithm of likelihood function, whose contour plot is presented below.  The estimated parameter values should be near the center of the contour plot.

likelihood[x_,y_]=Apply[Plus,Log[f /@ data]]/.{μ->x,σ->y};
ContourPlot[likelihood[x,y],{x,9,11},{y,.5,2},
	Axes->True,
	AxesLabel->{"μ","σ"}]

To find the maximum likelihood estimates, we first specify starting values that are fed to an equation solver which in this case is the built-in Mathematica function FindMinimum.  If one starting value for each parameter is used, then the FindMinimum function symbolically differentiates and then uses a steepest descent algorithm to find a minimum.  If two starting values are given, then it uses a different method, avoiding symbolic differentiation.  The ML function performs the estimation and then prints the output in a standard format.  First we recall the syntax for this function.

?ML

To use the function, it is important that the density match the data.  To check this, evaluate the density at the first observation in the data.

f[data[[1]] ]

This is the correct form and we now apply ML to our sample data.  If we take a range as our starting values, the FindMinimum command uses a derivative free algorithm to find the maximum.

startingvalues={{μ,8,12},{σ,.5,2}};
estimate=ML[data,f,startingvalues]

If the ML function fails to find an extremum of the likelihood function or if the parameters are not real-valued (assuming that the parameters are real), then choose new starting estimates for the parameter and recalculate.  Below we plot the original and estimated probability densities.  PlotDensities is a function that is defined in the MaximumLikelihood package.

f/.estimate

PlotDensities[data,f/.exact,f/.estimate]

To draw an inference about the true parameter values from our maximum likelihood estimates, we need an estimate of the variance of the estimators, i.e. some measure of how far they are away from the true parameters.  The variance matrix of the estimator is computed in the ML function by computing the gradient of the log density with respect to the unknown parameters evaluated at the maximum likelihood estimates, and summed over the data.  We can use the estimated variance (along with our knowledge that the limiting distribution of the estimator is a normal distribution) to compute confidence intervals centered on the estimates.

Plot2DConfidenceRegion[estimate]

The contours in the picture above represent three different confidence regions.  Since the confidence regions are themselves random variables, 99% of the samples will yield regions that contain the parameters,  μ and σ, and 1% of the samples will not.  The outer or largest ellipse is a 99% confidence region, the middle ellipse is a 95% confidence region, and the inner or smallest ellipse is a 90% confidence region.

To illustrate the large sample properties of the maximum likelihood estimator, we re-run this example with a sample size of 1000.

T=1000;
data=Table[Random[NormalDistribution[μ,σ]/.exact],{T}];
estimate=ML[data,f,startingvalues];
PlotDensities[data,f/.exact,f/.estimate ];

PlotDensities[data,f/.exact,f/.estimate]

We see that the estimated (dotted) density is closer to the true (solid) density and that the histogram of the data is also much closer to the true density.  The histogram is known as a nonparametric estimator of the density because it does not depend on specified parameters.

Example: Poisson Distribution

Next we consider the estimation of the mean λ for a sample data drawn from the Poisson distribution:

		f(x;λ)=(((e^(-λ))(λ^x))/(x!)),x=0,1,2,… .
		
We use the same approach to illustrate the method, generating a sample with a known mean and then estimating the mean from the sample, using the MaximumLikelihood package.  We draw sample data from a distribution with λ=10.

Clear[f];
T=1000;
exact=λ->10;
data=Table[Random[PoissonDistribution[λ]/.exact],{T}];
f=PDF[PoissonDistribution[λ],#]&;
f[x]

Our starting guess for λ is between 8 and 12.

startingvalues={{λ,8,12}};
estimate=ML[data,f,startingvalues]

To see how close we came to the true probability density function, we plot both the exact density, from which the data was drawn, and the estimated density.  We also plot a histogram of the data to illustrate how data affects the maximum likelihood estimator.

ListPlotDensities[data,f/.exact,f/.estimate]

The solid line is the theoretical distribution and the dashed line is the probabiltiy distribution with the estimated mean.

Example: Beta Distribution

Our next example is the Beta Distribution

		f(x;α,β)={({{((Γ(α+β)(((1-x))^(β-1))(x^(α-1)))/(Γ(α)Γ(β))),0<x<1}, {0, x∈ℝ-(0,1)}})
		
where α,β>0 are parameters and Γ denotes the gamma function.  As above, we draw sample data from the specified distribution.

Clear[f];
T=1000;
exact={α->2,β->5};
data=Table[Random[BetaDistribution[α,β]/.exact],{T}];
f=PDF[BetaDistribution[α,β],#]&;
f[x]

Next we apply the maximum likelihood algorithm with α=1 and β=4 as initial estimates.

startingvalues={{α,.5,2},{β,3,5}};
estimate=ML[data,f,startingvalues]

PlotDensities[data,f/.exact,f/.estimate]

Here is a picture of the confidence intervals for the parameters.

Plot2DConfidenceRegion[estimate]

The positive slope of the contours for this particular example indicates that the estimators of α and β are positively correlated.

Example: Normal Linear Regression Model

Let (Subscript[Y, i]) be a normally-distributed random variable with mean β (Subscript[x, 1i])+δ (Subscript[x, 2i]) and variance (σ^2).  The values (Subscript[x, 1i]) and (Subscript[x, 2i]) are nonstochastic.  Suppose we observe n independent draws from this distribution.  The estimation problem is to recover the unknown parameters, (β,δ,σ), from observed data {((Subscript[y, i]),(Subscript[x, 1i]),(Subscript[x, 2i])):i=1,2,…,n}.  The maximum likelihood estimator turns out to be equivalent to fitting a plane in 3-space through the data using a least squares criterion.  As we have done with the above examples, we generate data from known values of these parameters and then use the MaximumLikelihood package to estimate them.

Clear[f];
T=30;
X=Table[Random[MultinormalDistribution[{0,0},{{1,.8},{.8,1}}]],{T}];
exact={β->2,δ->2,σ->1};ε=Table[Random[NormalDistribution[0,σ]/.exact],{T}];
y=(X.{β,δ}  + ε) /.exact;
data=Transpose[{X,y}];
f=PDF[NormalDistribution[{β,δ}.#[[1]],σ],#[[2]] ]&;
f[data[[1]]]

We choose as starting estimates:  β=1, δ=3, and σ=0.5.

startingvalues={{β,.5,2.},{δ,2.,4.},{σ,.2,.9}};
estimate=ML[data,f,startingvalues];

Note that neither β nor δ is estimated very precisely and their standard errors are very large.  This is due to the high correlation between (Subscript[x, 1]) and (Subscript[x, 2]), and is called collinearity.  Because each variable is large when the other variable is large, it is very difficult to estimate their separate individual effects.

To compute a confidence region centered on the estimated slope coefficients, we note that the marginal distribution of β and δ is normal, and can be computed by dropping the σ terms from the estimated joint distribution. To accomplish this, we define estimateT by dropping all references to σ from the original estimate.

estimateT=Flatten[{Drop[estimate,-2],vcov->(Drop[#,-1]& /@Take[#,2])&[vcov/.estimate]}]

Plot2DConfidenceRegion[estimateT]

Note that the positive correlation between (Subscript[x, 1]) and (Subscript[x, 2]), leads to a negative correlation between their slope estimators. This implies that the sum of their coefficients is estimated very precisely, as determined by the the width of the ellipse, but their difference is not, as determined by the length of the ellipse.

Example: Multinomial Logit with Changing Choice Set

In this example, we estimate parameters in a Logit Discrete Choice Model with a changing choice set.  Ordinarily, setting up this problem is difficult because the data structure is not fixed, namely, the format of the sample data changes from data point to data point.  This does not pose a problem for Mathematica as long as the density is defined flexibly enough to map over whatever data structure you select.

The logit model is used to estimate descrete choices among a set of alternatives.  In this example, suppose that each data point represents an individual choosing among a set of alternatives, where the choice set changes from individual to individual.  Let X={(Subscript[x, i j])} be a set of observable choice attributes, like price or quality and assume that  y={(Subscript[y, i])} is vector of zeros with a single one, indicating which alternative is chosen.  We generate a distribution f(X,y;α,β) that depends on two parameters, α and β.

Clear[f,maxpos];
maxpos[x_]:=Map[If[#==Position[x,Max[x]][[1,1]],1,0]&,Range[Length[x]]];
T=100;
exact={α->2,β->2};
n:=Random[UniformDistribution[{2,8}]];X=Table[Table[Random[MultinormalDistribution[{0,0},{{1,.8},{.8,1}}]],{n}],{T}];
ε=Map[Table[Random[ExtremeValueDistribution[0,1]],{#}]&,Map[Length[#]&,X]];
y=Map[maxpos,Map[Map[({α,β}/. exact).#&,#]&,X]+ε] ;
data=Transpose[{X,y}];
f:=	Apply[Times,(Exp[#]/Apply[Plus,Exp[#]]& @Map[{α,β}.#&,#[[1]] ])^#[[2]]]&;

The probability density function for a three-choice, two-attribute decision is given by the following.

Clear[x,y];
f[{Array[x,{3,2}],Array[y,3]}]

The three expressions are the probabilities of choosing the three alternatives.  Since only one of the alternatives is chosen, two of the y[i] values are zero, and thus drop out of the likelihood.  As can be seen from the liklihood, the probabilities depend on the choice attributes of all the choices--even the alternatives not chosen.  Here are the first three observations of our data set.  Note that the dimensions of the choice set change from observation to observation.

Take[data,3] //TraditionalForm

The pure density function, which was defined above, can be mapped over the data set because it is held in a pure state and is evaluated only when it receives an argument. Next we estimate α and β using two starting values for each parameter.

startingvalues={{α,1,3},{β,1,3}};
estimate=ML[data,f,startingvalues]

Again we depict the the confidence regions as a graphic.

Plot2DConfidenceRegion[estimate]

Conclusion

The previous shortcomings of the MaximumLikelihood package (relatively slow speed and large memory requirements) have largely been overcome with the new functions built into Mathematica version 6.0.  Mathematica has a number of advantages for estimation.  Its fexibility allows it to address a wide range of models and to design new models with its object-oriented language.  This makes the package attractive for simulating and estimating non standard models.  Its ability to graph, simulate, and analyze the estimated distributions makes the package an attractive tool for teaching maximum likelihood estimation.

References

1.  Arthur Goldberger, A Course in Econometrics, Cambridge, Massachusetts:  Harvard University Press, 1991.

2.  William Greene, Econometric Analysis, 3rd ed., New Jersey:  Prentice Hall, 1997.

3.  Harold Larson, Introduction to the Theory of Statistics, New York:  John Wiley & Sons, 1973.

4.  Kenneth Train, Qualitative Choice Analysis, Cambridge, Massachusetts:  The MIT Press, 1986.