Maximum Likelihood Estimation

The Statistics Toolbox™ function mle is a convenient front end to the individual distribution fitting functions, and more. The function computes maximum likelihood estimates (MLEs) for distributions beyond those for which Statistics Toolbox software provides specific pdf functions.

For some pdfs, MLEs can be given in closed form and computed directly. For other pdfs, a search for the maximum likelihood must be employed. The search can be controlled with an options input argument, created using the statset function. For efficient searches, it is important to choose a reasonable distribution model and set appropriate convergence tolerances.

MLEs can be heavily biased, especially for small samples. As sample size increases, however, MLEs become unbiased minimum variance estimators with approximate normal distributions. This is used to compute confidence bounds for the estimates.

For example, consider the following distribution of means from repeated random samples of an exponential distribution:

mu = 1; % Population parameter
n = 1e3; % Sample size
ns = 1e4; % Number of samples

rng default  % For reproducibility
samples = exprnd(mu,n,ns); % Population samples
means = mean(samples); % Sample means

The Central Limit Theorem says that the means will be approximately normally distributed, regardless of the distribution of the data in the samples. The normfit function can be used to find the normal distribution that best fits the means:

[muhat,sigmahat,muci,sigmaci] = normfit(means)
muhat =

    1.0000


sigmahat =

    0.0315


muci =

    0.9994
    1.0006


sigmaci =

    0.0311
    0.0319

The function returns MLEs for the mean and standard deviation and their 95% confidence intervals.

To visualize the distribution of sample means together with the fitted normal distribution, you must scale the fitted pdf, with area = 1, to the area of the histogram being used to display the means:

numbins = 50;
hist(means,numbins)
h = findobj(gca,'Type','patch');
h.FaceColor = [.8 .8 1];
hold on
[bincounts,binpositions] = hist(means,numbins);
binwidth = binpositions(2) - binpositions(1);
histarea = binwidth*sum(bincounts);
x = binpositions(1):0.001:binpositions(end);
y = normpdf(x,muhat,sigmahat);
plot(x,histarea*y,'r','LineWidth',2)

Was this topic helpful?