Maximum Likelihood Estimation

Get access to whole semester resourcesSemester Pass

Chapter Sixty-Eight

Syllabus topic Module 2, "Maximum Likelihood Estimation"

Pages 395 to 403 of 591

In one line

Choose the parameter value that makes the data you actually observed as probable as possible.

In the wording a student can write in an examination: given a model with unknown parameters and some observed data, the likelihood is the probability of that data as a function of the parameters. The maximum likelihood estimate, or MLE, is the parameter value that maximises it.

Likelihood is not probability

The distinction papers ask for, and the one students state loosely.

	Probability	Likelihood
Varies over	the data, with the parameter fixed	the parameter, with the data fixed
Sums or integrates to 1	yes, over the data	no
Is a distribution over the parameter	no	no

For seven heads in ten flips the likelihood is p7 * (1-p)3. Read as a function of p it is not a probability distribution over p: it does not integrate to 1, and saying "the probability that p is 0.7" is a different claim, which needs a prior and belongs to Bayes Rule and Its Use. Maximum likelihood makes no statement about the probability of a parameter.

The recipe

Four steps, and every MLE derivation in the subject is these four.

Write the probability of the observed data under the model, as a function of the parameters. That is the likelihood L.
Take its logarithm, the log likelihood. The maximum is in the same place because the logarithm is increasing.
Differentiate with respect to each parameter and set the derivative to zero.
Solve, and check it is a maximum.

Worked for a coin. With h heads and t tails:

L(p) = ph * (1 - p)t

ln L(p) = h ln(p) + t ln(1 - p)

d/dp = h/p - t/(1 - p) = 0

=> h (1 - p) = t p

=> p = h / (h + t)

So the MLE of a coin's bias is the observed proportion of heads, which is the answer anyone would have guessed. That is the point: maximum likelihood derives the obvious estimate rather than assuming it, and the same four steps then give answers where nothing is obvious.

The computation

# Maximum likelihood estimation, computed: the likelihood of a coin searched by
# hand, why the LOG is used, the closed forms, and the bias the method carries.
import math

def lcg(seed):
    x = seed
    while True:
        x = (1664525 * x + 1013904223) % (2 ** 32)
        yield x / 2 ** 32

gen = lcg(9109)
rnd = lambda: next(gen)                                   # noqa: E731

# ---- 1. the likelihood of a coin, searched ---------------------------------
HEADS, TAILS = 7, 3
n = HEADS + TAILS
print("A COIN flipped %d times: %d heads, %d tails." % (n, HEADS, TAILS))
print("the parameter is p, the probability of heads. the LIKELIHOOD of the data")
print("for a given p is p**%d * (1-p)**%d. it is a function of p, NOT of the" % (HEADS, TAILS))
print("data, and it is not a probability distribution over p.")
print()
print("      p | likelihood        | log likelihood")
for i in range(1, 10):
    p = i / 10
    lik = p ** HEADS * (1 - p) ** TAILS
    print("   %.2f | %.12f | %14.6f" % (p, lik, math.log(lik)))
print()
best = max((i / 1000 for i in range(1, 1000)),
           key=lambda p: p ** HEADS * (1 - p) ** TAILS)
print("   searched on a grid of 1/1000, the maximum is at p = %.3f" % best)
print("   and the calculus says it is exactly heads/n = %d/%d = %.4f."
      % (HEADS, n, HEADS / n))
print("   differentiate the LOG likelihood, %d*ln(p) + %d*ln(1-p), set it to"
      % (HEADS, TAILS))
print("   zero: %d/p = %d/(1-p), so p = %d/%d. THAT is the MLE."
      % (HEADS, TAILS, HEADS, n))
print()

# ---- 2. why the log, measured ----------------------------------------------
print("WHY THE LOG, and it is not only convenience. the likelihood of a long")
print("sequence is a PRODUCT of many numbers below 1:")
for k in (10, 100, 1000, 1200):
    lik = 0.5 ** k
    print("   0.5 ** %-5d = %-24s log = %12.4f" % (k, repr(lik), k * math.log(0.5)))
print("   the fourth line is ZERO in floating point: the likelihood of 1200")
print("   coin flips UNDERFLOWS and the maximisation collapses. the log is a")
print("   sum, stays in range, and has its maximum in the same place because")
print("   the logarithm is increasing.")
print()

# ---- 3. a normal distribution: the closed forms ----------------------------
DATA = [62.0, 71.0, 58.0, 79.0, 66.0, 74.0, 69.0, 61.0]
print("A NORMAL DISTRIBUTION fitted to %d marks: %s"
      % (len(DATA), " ".join("%.0f" % v for v in DATA)))
m = len(DATA)
mu = sum(DATA) / m
var_n = sum((v - mu) ** 2 for v in DATA) / m
var_n1 = sum((v - mu) ** 2 for v in DATA) / (m - 1)
print("   the MLE of the mean is the sample mean          = %.4f" % mu)
print("   the MLE of the variance divides by n            = %.4f" % var_n)
print("   the UNBIASED estimate divides by n - 1          = %.4f" % var_n1)
print("   so the MLE standard deviation is %.4f and the unbiased %.4f."
      % (var_n ** 0.5, var_n1 ** 0.5))
print()

# ---- 4. the MLE variance is BIASED, and here is the measurement ------------
print("THE MLE OF A VARIANCE IS BIASED, and the bias is not a rounding matter.")
print("draw 20000 samples of 5 from a population whose true variance is known,")
print("estimate it both ways, and average the estimates.")
TRUE_MU, TRUE_SD = 70.0, 12.0

def normal():
    """Box-Muller, from the same deterministic stream."""
    u1, u2 = rnd(), rnd()
    if u1 < 1e-12:
        u1 = 1e-12
    return TRUE_MU + TRUE_SD * math.sqrt(-2 * math.log(u1)) * math.cos(2 * math.pi * u2)

trials, size = 20000, 5
sum_n = sum_n1 = 0.0
for _ in range(trials):
    s = [normal() for _ in range(size)]
    mean = sum(s) / size
    ss = sum((v - mean) ** 2 for v in s)
    sum_n += ss / size
    sum_n1 += ss / (size - 1)
print("   true variance                                   = %.4f" % (TRUE_SD ** 2))
print("   average of the MLE estimate (divide by n)       = %.4f" % (sum_n / trials))
print("   average of the estimate that divides by n - 1   = %.4f" % (sum_n1 / trials))
print("   the MLE is about %.1f per cent too small, and 1/n against 1/(n-1) at"
      % (100 * (1 - (sum_n / trials) / TRUE_SD ** 2)))
print("   n = %d is exactly %.1f per cent. the MLE reuses the sample mean, so"
      % (size, 100 * (1 - (size - 1) / size)))
print("   the deviations it squares are measured from the wrong centre and come")
print("   out too small. MAXIMUM LIKELIHOOD IS NOT THE SAME AS UNBIASED.")
print()

# ---- 5. the zero count -----------------------------------------------------
print("THE OTHER FAILURE: A COUNT OF ZERO. five flips, no heads.")
h, t = 0, 5
print("   MLE: p = %d/%d = %.4f" % (h, h + t, h / (h + t)))
print("   the estimate says heads is IMPOSSIBLE, on five flips. and any later")
print("   sequence containing a head then has likelihood exactly 0, so the")
print("   model cannot be used at all.")
print()
print("   LAPLACE SMOOTHING adds one imagined observation of each outcome:")
print("      p = (%d + 1)/(%d + 2) = %.4f" % (h, h + t, (h + 1) / (h + t + 2)))
print("   which is small but not impossible. more generally, add k:")
for k in (1, 2, 5):
    print("      k = %d:  p = (%d + %d)/(%d + %d) = %.4f"
          % (k, h, k, h + t, 2 * k, (h + k) / (h + t + 2 * k)))
print("   the effect of k fades as the data grows, which is what it should do.")
print("   the same repair is what the naive Bayes chapter needed for an unseen")
print("   word, and it is the point at which pure maximum likelihood is given up.")

Maximum Likelihood Estimation

A COIN flipped 10 times: 7 heads, 3 tails.
the parameter is p, the probability of heads. the LIKELIHOOD of the data
for a given p is p**7 * (1-p)**3. it is a function of p, NOT of the
data, and it is not a probability distribution over p.

      p | likelihood        | log likelihood
   0.10 | 0.000000072900 |     -16.434177
   0.20 | 0.000006553600 |     -11.935496
   0.30 | 0.000075014100 |      -9.497834
   0.40 | 0.000353894400 |      -7.946512
   0.50 | 0.000976562500 |      -6.931472
   0.60 | 0.001791590400 |      -6.324652
   0.70 | 0.002223566100 |      -6.108643
   0.80 | 0.001677721600 |      -6.390319
   0.90 | 0.000478296900 |      -7.645279

   searched on a grid of 1/1000, the maximum is at p = 0.700
   and the calculus says it is exactly heads/n = 7/10 = 0.7000.
   differentiate the LOG likelihood, 7*ln(p) + 3*ln(1-p), set it to
   zero: 7/p = 3/(1-p), so p = 7/10. THAT is the MLE.

WHY THE LOG, and it is not only convenience. the likelihood of a long
sequence is a PRODUCT of many numbers below 1:
   0.5 ** 10    = 0.0009765625             log =      -6.9315
   0.5 ** 100   = 7.888609052210118e-31    log =     -69.3147
   0.5 ** 1000  = 9.332636185032189e-302   log =    -693.1472
   0.5 ** 1200  = 0.0                      log =    -831.7766
   the fourth line is ZERO in floating point: the likelihood of 1200
   coin flips UNDERFLOWS and the maximisation collapses. the log is a
   sum, stays in range, and has its maximum in the same place because
   the logarithm is increasing.

A NORMAL DISTRIBUTION fitted to 8 marks: 62 71 58 79 66 74 69 61
   the MLE of the mean is the sample mean          = 67.5000
   the MLE of the variance divides by n            = 44.2500
   the UNBIASED estimate divides by n - 1          = 50.5714
   so the MLE standard deviation is 6.6521 and the unbiased 7.1114.

THE MLE OF A VARIANCE IS BIASED, and the bias is not a rounding matter.
draw 20000 samples of 5 from a population whose true variance is known,
estimate it both ways, and average the estimates.
   true variance                                   = 144.0000
   average of the MLE estimate (divide by n)       = 114.9999
   average of the estimate that divides by n - 1   = 143.7499
   the MLE is about 20.1 per cent too small, and 1/n against 1/(n-1) at
   n = 5 is exactly 20.0 per cent. the MLE reuses the sample mean, so
   the deviations it squares are measured from the wrong centre and come
   out too small. MAXIMUM LIKELIHOOD IS NOT THE SAME AS UNBIASED.

THE OTHER FAILURE: A COUNT OF ZERO. five flips, no heads.
   MLE: p = 0/5 = 0.0000
   the estimate says heads is IMPOSSIBLE, on five flips. and any later
   sequence containing a head then has likelihood exactly 0, so the
   model cannot be used at all.

   LAPLACE SMOOTHING adds one imagined observation of each outcome:
      p = (0 + 1)/(5 + 2) = 0.1429
   which is small but not impossible. more generally, add k:
      k = 1:  p = (0 + 1)/(5 + 2) = 0.1429
      k = 2:  p = (0 + 2)/(5 + 4) = 0.2222
      k = 5:  p = (0 + 5)/(5 + 10) = 0.3333
   the effect of k fades as the data grows, which is what it should do.
   the same repair is what the naive Bayes chapter needed for an unseen
   word, and it is the point at which pure maximum likelihood is given up.

Maximum Likelihood Estimation

Reading it

The grid search agrees with the calculus. The likelihood column rises to 0.002223566100 at p = 0.70 and falls either side, and a search on a grid of one thousandth puts the maximum at exactly 0.700, which is 7/10. The derivation and the arithmetic meet.

Maximum Likelihood Estimation

The log is not a convenience. Read the four underflow lines:

Maximum Likelihood Estimation

	Likelihood	Log likelihood
`0.5 ** 10`	`0.0009765625`	-6.9315
`0.5 ** 100`	`7.888609052210118e-31`	-69.3147
`0.5 ** 1000`	`9.332636185032189e-302`	-693.1472
`0.5 ** 1200`	`0.0`	-831.7766

The likelihood of 1200 coin flips is zero in floating point. Not small: zero. Every candidate parameter would score 0, the comparison between them is destroyed, and the maximisation returns whatever came first. The log likelihood of the same data is -831.7766, an ordinary number. That is why every real implementation works in logs, and it is a good answer to "why the log likelihood" that goes beyond "sums are easier than products".

The normal distribution's MLEs are worth memorising, since a paper may simply ask for them: the MLE of the mean is the sample mean, 67.5000 here, and the MLE of the variance is the average squared deviation from it, 44.2500, dividing by n.

The first failure: maximum likelihood is biased

The estimate that divides by n is the maximum likelihood one. The estimate that divides by n - 1 is the one every statistics course teaches. They are not the same, and the difference is measured here, not asserted.

Twenty thousand samples of five, drawn from a population whose true variance is 144.0000:

	Average estimate
MLE, dividing by `n`	114.9999
dividing by `n - 1`	143.7499

The MLE is 20.1 per cent too small, and 1 - (n-1)/n at n = 5 is exactly 20.0 per cent. The measurement lands on the theory.

The cause, which is the part worth marks. The MLE measures every deviation from the sample mean, and the sample mean is itself pulled towards the sample. So the deviations are taken from a centre that is already too close to the data, and the squares come out too small. Dividing by n - 1 compensates for the one degree of freedom spent on estimating the mean.

Maximum Likelihood Estimation

So: maximum likelihood is not the same as unbiased. It is the estimate that best explains the data in hand, not the estimate that is right on average over repeated samples. And the bias vanishes as n grows, which is why nobody minds much on large data and everybody minds on small.

The second failure: a count of zero

Five flips, no heads. The MLE is 0/5 = 0.0000, which asserts that heads is impossible on the evidence of five flips. Worse, it is not merely overconfident: any later sequence containing a head has likelihood exactly 0, so the model assigns probability zero to something that just happened and cannot be used at all.

The repair is Laplace smoothing, also called add-one: pretend to have seen one of each outcome before starting.

p = (h + k) / (h + t + 2k) k = 1 is Laplace, add-one

`k`	Estimate
0, plain MLE	0.0000
1	0.1429
2	0.2222
5	0.3333

Note that the choice of k matters a great deal on five observations and hardly at all on five thousand, which is exactly the right behaviour: the imagined observations are outvoted by real ones.

And be honest about what has happened: adding imagined counts is no longer maximum likelihood. It is the maximum a posteriori estimate under a prior that says extreme values of p are unlikely, and the constant k is that prior's strength. Naive Bayes needed the same repair for a word never seen in a class, and this is why.

MLE, MAP and Bayesian, in one table

A paper asking how maximum likelihood relates to Bayes wants these three separated.

	Maximum likelihood	Maximum a posteriori	Fully Bayesian
Maximises	`P(data given parameter)`	`P(data given parameter) * P(parameter)`	nothing; it keeps the whole distribution
Uses a prior	no	yes	yes
Returns	one value	one value	a distribution over parameters
On 0 heads in 5	`p = 0`	`p = 0.1429` with add-one	a distribution with little mass near 0.5
Cost	cheapest	cheap	expensive

MLE is MAP with a uniform prior. If every parameter value is equally likely beforehand, the prior is a constant, and maximising the product is the same as maximising the likelihood. That one line is worth stating.

The properties, honestly

What maximum likelihood is good for, and what it is not.

Property	Holds
Consistent: converges on the true parameter as `n` grows	yes
Asymptotically efficient: no estimator does better for large `n`	yes
Invariant: the MLE of a function of the parameter is that function of the MLE	yes
Unbiased	no, as measured above
Well behaved on small samples	no
Well behaved on zero counts	no

Maximum Likelihood Estimation

Distinctions

	Probability	Likelihood
Fixed	the parameter	the data
Varies	the data	the parameter
Normalised	yes	no

	MLE of a variance	Unbiased estimate
Divides by	`n`	`n - 1`
On the eight marks	44.2500	50.5714
Averaged over 20000 samples of 5	114.9999	143.7499
True value	144.0000	144.0000

	Plain MLE	With add-one
0 heads in 5	0.0000	0.1429
Can score a head afterwards	no, likelihood 0	yes
Is it still maximum likelihood	yes	no, it is MAP

What it does not mean

The likelihood is not the probability of the parameter. It is the probability of the data, read as a function of the parameter.

The log is not just for convenience. Without it the likelihood of 1200 flips is 0.0 and the maximisation is destroyed.

Maximum likelihood is not unbiased. Its variance estimate ran 20.1 per cent short of a true 144.

A zero count is not a small probability. It is zero, and it makes the model unusable.

Laplace smoothing is not maximum likelihood. It is MAP under a prior, and saying so is part of a correct answer.

A larger likelihood is not a better model. Adding parameters raises the likelihood always, which is why model comparison needs a penalty for complexity and not a likelihood alone.

Quick revision

Likelihood is P(data given parameter) read as a function of the parameter. It is not normalised and is not a distribution over the parameter.
The recipe: write L, take ln L, differentiate, set to zero, solve.
Coin: L = ph (1-p)t, ln L = h ln p + t ln(1-p), and the MLE is h/(h+t). Measured: 7 heads in 10 peaks at p = 0.700.
Why the log: 0.5 1200 is 0.0** in floating point; its log likelihood is -831.7766. Logs also turn products into sums and keep the maximum in the same place.
Normal: MLE mean = the sample mean (67.5000); MLE variance = average squared deviation, dividing by n (44.2500), against 50.5714 for n - 1.
The MLE variance is biased: over 20000 samples of 5 with a true variance of 144.0000, it averaged 114.9999, 20.1 per cent short, matching 1 - (n-1)/n = 20%. Cause: the deviations are taken from the sample mean, which is pulled towards the sample.
Zero counts: 0 heads in 5 gives p = 0, so heads becomes impossible and any head later has likelihood 0. Laplace smoothing gives (h+k)/(h+t+2k), which is 0.1429 at k = 1, and is MAP, not MLE.
MLE is MAP with a uniform prior. MAP returns one value with a prior; fully Bayesian keeps the whole distribution.
MLE is consistent, asymptotically efficient and invariant, and is not unbiased and not reliable on small samples or zero counts.

Maximum Likelihood Estimation

Test yourself

1. Distinguish likelihood from probability. Probability varies over the data with the parameter fixed and sums to one over the data. Likelihood is the same expression read as a function of the parameter with the data fixed; it does not sum or integrate to one and is not a distribution over the parameter.

2. Derive the maximum likelihood estimate of a coin's bias from h heads and t tails. The likelihood is ph * (1-p)t. Its logarithm is h ln p + t ln(1-p), whose derivative is h/p - t/(1-p). Setting that to zero gives h(1-p) = tp, so p = h/(h+t), the observed proportion of heads.

3. Give two reasons for maximising the log likelihood rather than the likelihood. Products of many probabilities underflow to zero in floating point, and 0.5 to the power 1200 is exactly 0.0 on this machine, which destroys the comparison between candidate parameters; the log of the same quantity is about -831.78 and computes normally. The log also turns the product into a sum, which differentiates term by term, and because the logarithm is increasing the maximum is in the same place.

4. State the MLEs for a normal distribution and say which one is biased. The MLE of the mean is the sample mean. The MLE of the variance is the average squared deviation from that mean, dividing by n. The variance estimate is biased low; the unbiased version divides by n - 1.

5. Explain why the MLE of a variance is too small, and quantify it at n = 5. The squared deviations are taken from the sample mean, which is itself pulled towards the sample, so they are smaller than deviations from the true mean. The estimate is short by a factor of (n-1)/n, which at n = 5 is 20 per cent. Measured over 20000 samples of five from a population with variance 144, the MLE averaged 115.0 against 143.7 for the n - 1 version.

6. Five coin flips give no heads. Give the MLE, say what is wrong with it, and repair it. The MLE is 0/5 = 0, which asserts heads is impossible on five flips, and any later sequence containing a head then has likelihood exactly zero, so the model cannot be used. Laplace smoothing adds one imagined observation of each outcome, giving (0+1)/(5+2) = 0.1429, which is small without being impossible, and the added counts are outvoted as real data accumulates.

Maximum Likelihood Estimation

7. How do MLE, MAP and a fully Bayesian treatment differ? MLE maximises the probability of the data given the parameter and uses no prior. MAP maximises the product of that likelihood and a prior over the parameter, returning a single value; MLE is the special case where the prior is uniform. A fully Bayesian treatment returns no single value at all but the whole posterior distribution over the parameter, at a greater computational cost.

The rest of this subject

These notes are cut from the University's printed syllabus. Open the syllabus itself, or the past papers, for the same subject.

Issue