From dbc16354d1dfd0035b56ca15eb33e29f385f7010 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 29 Sep 2026 07:41:13 +1000 Subject: [PATCH 1/5] mle: make the lecture a sequel to fitting_distributions - Compare the lognormal MLE with the method-of-moments fit, using the KS statistic - Replace the by-eye tail comparison with log-scale Q-Q plots and KS - Show why the method of moments fails for a Pareto tail with b <= 1 - Add a Student's t MLE fit to Amazon returns, as promised in fitting_distributions - Fix the lognormal density, the sign of the derivative with respect to mu, and use \partial and aligned - Add figure captions, axis labels and lw=2 per the style guide Co-Authored-By: Claude Opus 5.5 --- lectures/mle.md | 395 +++++++++++++++++++++++++++++++++++++----------- 1 file changed, 310 insertions(+), 85 deletions(-) diff --git a/lectures/mle.md b/lectures/mle.md index a740b035..25601c33 100644 --- a/lectures/mle.md +++ b/lectures/mle.md @@ -13,13 +13,23 @@ kernelspec: # Maximum Likelihood Estimation +In addition to what's in Anaconda, this lecture will need the following libraries: + +```{code-cell} ipython3 +:tags: [hide-output] + +!pip install --upgrade yfinance +``` + ```{code-cell} ipython3 from scipy.stats import lognorm, pareto, expon +from scipy.stats import t as student_t +from scipy.stats import kurtosis import numpy as np from scipy.integrate import quad import matplotlib.pyplot as plt import pandas as pd -from math import exp +import yfinance as yf ``` ## Introduction @@ -94,7 +104,7 @@ url = 'https://raw.githubusercontent.com/QuantEcon/data-lectures/main/lectures/S df = pd.read_csv(url) df = df.dropna() df = df[df['year'] == 2016] -df = df.loc[df['n_wealth'] > 1 ] #restrcting data to net worth > 1 +df = df.loc[df['n_wealth'] > 1 ] # restricting data to net worth > 1 rv = df['n_wealth'].sample(n=n, random_state=1234) rv = rv.to_numpy() / 100_000 sample = rv @@ -103,14 +113,20 @@ sample = rv Let's histogram this sample. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Histogram of US household wealth + name: fig:mle-wealth-hist +--- fig, ax = plt.subplots() ax.set_xlim(-1, 20) density, edges = np.histogram(sample, bins=5000, density=True) prob = density * np.diff(edges) -plt.stairs(prob, edges, fill=True, alpha=0.8, label=r"unit: $\$100,000$") -plt.ylabel("prob") -plt.xlabel("net wealth") -plt.legend() +ax.stairs(prob, edges, fill=True, alpha=0.8, label=r"unit: $\$100,000$") +ax.set_ylabel("prob") +ax.set_xlabel("net wealth") +ax.legend() plt.show() ``` @@ -168,9 +184,17 @@ histogram log wealth instead of wealth, the picture starts to look something like a bell-shaped curve. ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Histogram of log wealth + name: fig:mle-log-wealth-hist +--- ln_sample = np.log(sample) fig, ax = plt.subplots() ax.hist(ln_sample, density=True, bins=200, histtype='stepfilled', alpha=0.8) +ax.set_xlabel("log of net wealth") +ax.set_ylabel("density") plt.show() ``` @@ -183,15 +207,15 @@ data. The pdf of a lognormally distributed random variable $X$ is given by: $$ - f(x, \mu, \sigma) - = \frac{1}{x}\frac{1}{\sigma \sqrt{2\pi}} - \exp\left(\frac{-1}{2}\left(\frac{\ln x-\mu}{\sigma}\right)\right)^2 + f(x; \mu, \sigma) + = \frac{1}{x \sigma \sqrt{2\pi}} + \exp \left( - \frac{(\ln x - \mu)^2}{2 \sigma^2} \right) $$ For our sample $w_1, w_2, \cdots, w_n$, the [likelihood function](https://en.wikipedia.org/wiki/Likelihood_function) is given by $$ - L(\mu, \sigma | w_i) = \prod_{i=1}^{n} f(w_i, \mu, \sigma) + L(\mu, \sigma) = \prod_{i=1}^{n} f(w_i; \mu, \sigma) $$ The likelihood function can be viewed as both @@ -203,35 +227,39 @@ Taking logs on both sides gives us the log likelihood function, which is $$ \begin{aligned} - \ell(\mu, \sigma | w_i) - & = \ln \left[ \prod_{i=1}^{n} f(w_i, \mu, \sigma) \right] \\ + \ell(\mu, \sigma) + & = \ln \left[ \prod_{i=1}^{n} f(w_i; \mu, \sigma) \right] \\ & = -\sum_{i=1}^{n} \ln w_i - \frac{n}{2} \ln(2\pi) - \frac{n}{2} \ln \sigma^2 - \frac{1}{2\sigma^2} \sum_{i=1}^n (\ln w_i - \mu)^2 \end{aligned} $$ -To find where this function is maximised we find its partial derivatives wrt $\mu$ and $\sigma ^2$ and equate them to $0$. +To find where this function is maximized we find its partial derivatives with respect to $\mu$ and $\sigma^2$ and set them to zero. Let's first find the maximum likelihood estimate (MLE) of $\mu$ $$ -\frac{\delta \ell}{\delta \mu} - = - \frac{1}{2\sigma^2} \times 2 \sum_{i=1}^n (\ln w_i - \mu) = 0 \\ -\implies \sum_{i=1}^n \ln w_i - n \mu = 0 \\ -\implies \hat{\mu} = \frac{\sum_{i=1}^n \ln w_i}{n} +\begin{aligned} +\frac{\partial \ell}{\partial \mu} + = \frac{1}{\sigma^2} \sum_{i=1}^n (\ln w_i - \mu) = 0 + & \implies \sum_{i=1}^n \ln w_i - n \mu = 0 \\ + & \implies \hat{\mu} = \frac{1}{n} \sum_{i=1}^n \ln w_i +\end{aligned} $$ Now let's find the MLE of $\sigma$ $$ -\frac{\delta \ell}{\delta \sigma^2} +\begin{aligned} +\frac{\partial \ell}{\partial \sigma^2} = - \frac{n}{2\sigma^2} + \frac{1}{2\sigma^4} - \sum_{i=1}^n (\ln w_i - \mu)^2 = 0 \\ - \implies \frac{n}{2\sigma^2} = + \sum_{i=1}^n (\ln w_i - \mu)^2 = 0 + & \implies \frac{n}{2\sigma^2} = \frac{1}{2\sigma^4} \sum_{i=1}^n (\ln w_i - \mu)^2 \\ - \implies \hat{\sigma} = - \left( \frac{\sum_{i=1}^{n}(\ln w_i - \hat{\mu})^2}{n} \right)^{1/2} + & \implies \hat{\sigma} = + \left( \frac{1}{n} \sum_{i=1}^{n}(\ln w_i - \hat{\mu})^2 \right)^{1/2} +\end{aligned} $$ Now that we have derived the expressions for $\hat{\mu}$ and $\hat{\sigma}$, @@ -251,21 +279,89 @@ num = (ln_sample - μ_hat)**2 Let's plot the lognormal pdf using the estimated parameters against our sample data. ```{code-cell} ipython3 -dist_lognorm = lognorm(σ_hat, scale = exp(μ_hat)) -x = np.linspace(0,50,10000) +--- +mystnb: + figure: + caption: Lognormal fit to wealth + name: fig:mle-lognormal-fit +--- +dist_lognorm = lognorm(σ_hat, scale=np.exp(μ_hat)) +x = np.linspace(0, 50, 10000) fig, ax = plt.subplots() -ax.set_xlim(-1,20) +ax.set_xlim(-1, 20) ax.hist(sample, density=True, bins=5_000, histtype='stepfilled', alpha=0.5) -ax.plot(x, dist_lognorm.pdf(x), 'k-', lw=0.5, label='lognormal pdf') +ax.plot(x, dist_lognorm.pdf(x), 'k-', lw=2, label='lognormal pdf') +ax.set_xlabel("net wealth") +ax.set_ylabel("density") ax.legend() plt.show() ``` Our estimated lognormal distribution appears to be a reasonable fit for the overall data. -We now use {eq}`eq:est_rev` to calculate total revenue. + +### Comparison with the method of moments + +Notice that $\hat{\mu}$ and $\hat{\sigma}$ are just the sample mean and sample +standard deviation of log wealth. + +So, for the lognormal class, maximum likelihood gives the same answer as +fitting a normal distribution to log wealth by the method of moments. + +This is *not* the same as the lognormal fit in {doc}`fitting_distributions`, +which matches the mean and variance of wealth itself, rather than of log +wealth. + +Let's compute that fit too. + +```{code-cell} ipython3 +def fit_lognormal_mom(sample): + "Lognormal fit matching the sample mean and variance of the data." + m, v = sample.mean(), sample.var() + σ_squared = np.log(1 + v / m**2) + μ = np.log(m) - σ_squared / 2 + return lognorm(s=np.sqrt(σ_squared), scale=np.exp(μ)) + +dist_lognorm_mom = fit_lognormal_mom(sample) +``` + +To compare the two fits we use the Kolmogorov-Smirnov statistic from +{doc}`fitting_distributions`, which is the largest vertical gap between the +ECDF of the data and the fitted CDF. + +```{code-cell} ipython3 +def ks_statistic(sample, u): + "Largest vertical distance between the ECDF of the sample and the CDF of u." + x_sorted = np.sort(sample) + n = len(x_sorted) + F = u.cdf(x_sorted) + above = np.arange(1, n+1) / n - F # gap just after each jump + below = F - np.arange(0, n) / n # gap just before each jump + return max(above.max(), below.max()) +``` + +```{code-cell} ipython3 +pd.DataFrame({ + 'σ': [σ_hat, dist_lognorm_mom.kwds['s']], + 'KS statistic': [ks_statistic(sample, dist_lognorm), + ks_statistic(sample, dist_lognorm_mom)]}, + index=['maximum likelihood', 'method of moments']) +``` + +The two methods give quite different parameters, and maximum likelihood gives +a much closer fit. + +The reason is that the sample variance of wealth is dominated by a handful of +very wealthy households, so matching it pulls the whole fitted distribution +toward them. + +Maximum likelihood works with log wealth, where these households are far less +extreme. + +We now use {eq}`eq:est_rev` and the maximum likelihood fit to calculate total +revenue. We will compute the integral using numerical integration via SciPy's [quad](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) @@ -301,7 +397,7 @@ with parameters $b$ and $x_m$. In this case, the maximum likelihood estimates are known to be $$ - \hat{b} = \frac{n}{\sum_{i=1}^{n} \ln (w_i/\hat{x_m})} + \hat{b} = \frac{n}{\sum_{i=1}^{n} \ln (w_i/\hat{x}_m)} \quad \text{and} \quad \hat{x}_m = \min_{i} w_i $$ @@ -322,7 +418,7 @@ b_hat Now let's recompute total revenue. ```{code-cell} ipython3 -dist_pareto = pareto(b = b_hat, scale = xm_hat) +dist_pareto = pareto(b=b_hat, scale=xm_hat) tr_pareto = total_revenue(dist_pareto) tr_pareto ``` @@ -340,31 +436,42 @@ We see that choosing the right distribution is extremely important. Let's compare the fitted Pareto distribution to the histogram: ```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Pareto fit to wealth + name: fig:mle-pareto-fit +--- fig, ax = plt.subplots() ax.set_xlim(-1, 20) -ax.set_ylim(0,1.75) +ax.set_ylim(0, 1.75) ax.hist(sample, density=True, bins=5_000, histtype='stepfilled', alpha=0.5) -ax.plot(x, dist_pareto.pdf(x), 'k-', lw=0.5, label='Pareto pdf') +ax.plot(x, dist_pareto.pdf(x), 'k-', lw=2, label='Pareto pdf') +ax.set_xlabel("net wealth") +ax.set_ylabel("density") ax.legend() plt.show() ``` -We observe that in this case the fit for the Pareto distribution is not very -good, so we can probably reject it. +The fit is poor, and the KS statistic confirms it. -## What is the best distribution? +```{code-cell} ipython3 +ks_statistic(sample, dist_pareto), ks_statistic(sample, dist_lognorm) +``` -There is no "best" distribution --- every choice we make is an assumption. +For the sample as a whole, the lognormal distribution is clearly preferred. -All we can do is try to pick a distribution that fits the data well. -The plots above suggested that the lognormal distribution is optimal. +## Fitting the right tail -However when we inspect the upper tail (the richest people), the Pareto distribution may be a better fit. +There is no "best" distribution --- every choice we make is an assumption. + +The lognormal distribution fits the sample as a whole, but for questions about +the richest households, what matters is the fit in the upper tail. -To see this, let's now set a minimum threshold of net worth in our dataset. +To study the tail, let's set a minimum threshold of net worth in our dataset. We set an arbitrary threshold of $500,000 and read the data into `sample_tail`. @@ -372,84 +479,202 @@ We set an arbitrary threshold of $500,000 and read the data into `sample_tail`. :tags: [hide-input] df_tail = df.loc[df['n_wealth'] > 500_000 ] -df_tail.head() rv_tail = df_tail['n_wealth'].sample(n=10_000, random_state=4321) rv_tail = rv_tail.to_numpy() sample_tail = rv_tail/500_000 ``` -Let's plot this data. +Let's fit a lognormal and a Pareto distribution to this data by maximum +likelihood, using the formulas derived above. ```{code-cell} ipython3 -fig, ax = plt.subplots() -ax.set_xlim(0,50) -ax.hist(sample_tail, density=True, bins=500, histtype='stepfilled', alpha=0.8) +ln_sample_tail = np.log(sample_tail) +μ_hat_tail = np.mean(ln_sample_tail) +σ_hat_tail = np.std(ln_sample_tail) +dist_lognorm_tail = lognorm(σ_hat_tail, scale=np.exp(μ_hat_tail)) + +xm_hat_tail = min(sample_tail) +b_hat_tail = 1/np.mean(np.log(sample_tail/xm_hat_tail)) +dist_pareto_tail = pareto(b=b_hat_tail, scale=xm_hat_tail) +``` + +Since we care about the tail, we compare the fits with Q-Q plots, which were +introduced in {doc}`fitting_distributions` and show clearly how a fit behaves +at the extremes. + +Tail wealth spans several orders of magnitude, so we draw the plots on log +scales. + +```{code-cell} ipython3 +def qq_plot(sample, u, ax, **kwargs): + "Plot sample quantiles against the quantiles of the distribution u." + x_sorted = np.sort(sample) + n = len(x_sorted) + p = (np.arange(1, n+1) - 0.5) / n + ax.plot(u.ppf(p), x_sorted, '.', ms=3, alpha=0.6, **kwargs) + lo, hi = u.ppf(p[0]), u.ppf(p[-1]) + ax.plot([lo, hi], [lo, hi], 'k--', lw=2) + ax.set_xlabel('fitted quantiles') + ax.set_ylabel('sample quantiles') +``` + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Q-Q plots for tail wealth + name: fig:mle-tail-qq +--- +fig, axes = plt.subplots(1, 2, figsize=(10, 4.5)) +for ax, u, label in zip(axes, + (dist_lognorm_tail, dist_pareto_tail), + ('lognormal', 'Pareto')): + qq_plot(sample_tail, u, ax) + ax.set_xscale('log') + ax.set_yscale('log') + ax.set_title(label) +plt.tight_layout() plt.show() ``` -Now let's try fitting some distributions to this data. +The lognormal fit bends away from the 45 degree line: it puts probability on +wealth below the threshold, which cannot occur in this sample, and it +underpredicts large fortunes. +The Pareto fit tracks the data closely until far into the tail, up to around 20 +times the threshold, or 10 million dollars. -### Lognormal distribution for the right hand tail +Beyond that point the data flatten out, because no household in our data set +is worth more than about 50 million dollars. -Let's start with the lognormal distribution +This reflects the survey rather than the population, since the SCF does not +capture the very wealthiest households. -We estimate the parameters again and plot the density against our data. +The KS statistic also favors the Pareto distribution. ```{code-cell} ipython3 -ln_sample_tail = np.log(sample_tail) -μ_hat_tail = np.mean(ln_sample_tail) -num_tail = (ln_sample_tail - μ_hat_tail)**2 -σ_hat_tail = (np.mean(num_tail))**(1/2) -dist_lognorm_tail = lognorm(σ_hat_tail, scale = exp(μ_hat_tail)) +ks_statistic(sample_tail, dist_lognorm_tail), ks_statistic(sample_tail, dist_pareto_tail) +``` -fig, ax = plt.subplots() -ax.set_xlim(0,50) -ax.hist(sample_tail, density=True, bins=500, histtype='stepfilled', alpha=0.5) -ax.plot(x, dist_lognorm_tail.pdf(x), 'k-', lw=0.5, label='lognormal pdf') -ax.legend() -plt.show() +So, while the lognormal distribution was a good fit for the entire dataset, the +Pareto distribution is a better fit for the right hand tail. + + +### Why not the method of moments? + +Here is the estimated Pareto parameter for the tail. + +```{code-cell} ipython3 +b_hat_tail ``` -While the lognormal distribution was a good fit for the entire dataset, -it is not a good fit for the right hand tail. +This value is below 1, which has a striking implication. + +A Pareto distribution with parameter $b$ has mean $b x_m / (b - 1)$ when +$b > 1$, and infinite mean when $b \leq 1$. + +So the fitted distribution has no finite mean, let alone a finite variance. + +(We discuss distributions of this kind in {doc}`heavy_tails`.) +This is a problem for the method of moments. -### Pareto distribution for the right hand tail +With $x_m$ fixed, matching the population mean to the sample mean $\bar w$ gives + +$$ +\hat{b} = \frac{\bar w}{\bar w - x_m} +$$ -Let's now assume the truncated dataset has a Pareto distribution. +which is greater than 1 for *every* data set, since $\bar w > x_m$. -We estimate the parameters again and plot the density against our data. +The method of moments therefore cannot recover a value of $b$ at or below 1, +however much data we have. + +Let's see what it gives here. ```{code-cell} ipython3 -xm_hat_tail = min(sample_tail) -den_tail = np.log(sample_tail/xm_hat_tail) -b_hat_tail = 1/np.mean(den_tail) -dist_pareto_tail = pareto(b = b_hat_tail, scale = xm_hat_tail) +w_bar_tail = np.mean(sample_tail) +w_bar_tail / (w_bar_tail - xm_hat_tail) +``` -fig, ax = plt.subplots() -ax.set_xlim(0, 50) -ax.set_ylim(0,0.65) -ax.hist(sample_tail, density=True, bins= 500, histtype='stepfilled', alpha=0.5) -ax.plot(x, dist_pareto_tail.pdf(x), 'k-', lw=0.5, label='pareto pdf') -plt.show() +The trouble is that the method of moments assumes the moments it matches +exist. + +Maximum likelihood makes no such assumption, since it works with the density +itself, which is well defined for every value of $b$. + +This is one reason why maximum likelihood is the standard method for fitting +heavy-tailed distributions. + + +## Heavy-tailed returns + +In {doc}`fitting_distributions` we fitted a Student's t distribution to monthly +returns on Amazon shares, using the method of moments. + +We noted there that the fitted value of the degrees of freedom $\nu$ was +unreliable, because it came from the sample kurtosis, a fourth moment. + +Let's fit the same class by maximum likelihood instead. + +```{code-cell} ipython3 +:tags: [hide-output] + +data = yf.download('AMZN', '2000-1-1', '2024-1-1', interval='1mo') +prices = data['Close']['AMZN'] +returns = prices.pct_change().dropna() * 100 ``` -The Pareto distribution is a better fit for the right hand tail of our dataset. +Unlike the lognormal and Pareto cases, there is no closed-form expression for +the maximum likelihood estimates of the Student's t distribution. -### So what is the best distribution? +Instead, the `fit` method in SciPy maximizes the log likelihood numerically. -As we said above, there is no "best" distribution --- each choice is an -assumption. +```{code-cell} ipython3 +ν_mle, loc_mle, scale_mle = student_t.fit(returns) +dist_t_mle = student_t(df=ν_mle, loc=loc_mle, scale=scale_mle) +ν_mle +``` -We just have to test what we think are reasonable distributions. +For comparison, here is the method of moments fit from {doc}`fitting_distributions`. -One test is to plot the data against the fitted distribution, as we did. +```{code-cell} ipython3 +def fit_t_mom(sample): + m, s = sample.mean(), sample.std() + ν = 4 + 6 / kurtosis(sample) + return student_t(df=ν, loc=m, scale=s * np.sqrt((ν - 2) / ν)) + +dist_t_mom = fit_t_mom(returns) +dist_t_mom.kwds['df'] +``` + +Maximum likelihood gives a noticeably smaller $\nu$, and hence heavier tails. + +It also gives a closer fit. + +```{code-cell} ipython3 +ks_statistic(returns, dist_t_mom), ks_statistic(returns, dist_t_mle) +``` + +The Q-Q plot confirms that the maximum likelihood fit tracks the extreme +returns well. + +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Amazon returns against a t fitted by MLE + name: fig:mle-qq-returns-t +--- +fig, ax = plt.subplots() +qq_plot(returns, dist_t_mle, ax) +plt.show() +``` -There are other more rigorous tests, such as the [Kolmogorov-Smirnov test](https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test). +The lesson is the same as for the wealth tail: when the tails are heavy, the +sample moments are dominated by a few extreme observations, and maximum +likelihood makes better use of the data. -We omit such advanced topics (but encourage readers to study them once -they have completed these lectures). ## Exercises @@ -481,7 +706,7 @@ $$ ``` ```{code-cell} ipython3 -dist_exp = expon(scale = 1/λ_hat) +dist_exp = expon(scale=1/λ_hat) tr_expo = total_revenue(dist_exp) tr_expo ``` @@ -507,7 +732,7 @@ fig, ax = plt.subplots() ax.set_xlim(-1, 20) ax.hist(sample, density=True, bins=5000, histtype='stepfilled', alpha=0.5) -ax.plot(x, dist_exp.pdf(x), 'k-', lw=0.5, label='exponential pdf') +ax.plot(x, dist_exp.pdf(x), 'k-', lw=2, label='exponential pdf') ax.legend() plt.show() From 90158313ab3bc33ccd8757ccdb7a83cda6af1fbf Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 29 Sep 2026 10:46:12 +1000 Subject: [PATCH 2/5] mle: add overview with role and history of maximum likelihood Co-Authored-By: Claude Opus 5.5 --- lectures/_static/quant-econ.bib | 29 +++++++++++++++++++ lectures/mle.md | 49 ++++++++++++++++++++++++++++++++- 2 files changed, 77 insertions(+), 1 deletion(-) diff --git a/lectures/_static/quant-econ.bib b/lectures/_static/quant-econ.bib index ece2dc72..397c4448 100644 --- a/lectures/_static/quant-econ.bib +++ b/lectures/_static/quant-econ.bib @@ -3119,3 +3119,32 @@ @article{jianakoplos1997wealth pages = {18--31}, year = {1997} } + +@article{fisher1922mathematical, + title = {On the mathematical foundations of theoretical statistics}, + author = {Fisher, Ronald A.}, + journal = {Philosophical Transactions of the Royal Society of London. Series A}, + volume = {222}, + pages = {309--368}, + year = {1922} +} + +@article{aldrich1997fisher, + title = {{R. A. Fisher} and the making of maximum likelihood 1912--1922}, + author = {Aldrich, John}, + journal = {Statistical Science}, + volume = {12}, + number = {3}, + pages = {162--176}, + year = {1997} +} + +@article{stigler2007epic, + title = {The epic story of maximum likelihood}, + author = {Stigler, Stephen M.}, + journal = {Statistical Science}, + volume = {22}, + number = {4}, + pages = {598--620}, + year = {2007} +} diff --git a/lectures/mle.md b/lectures/mle.md index 25601c33..af6ad3ae 100644 --- a/lectures/mle.md +++ b/lectures/mle.md @@ -32,7 +32,54 @@ import pandas as pd import yfinance as yf ``` -## Introduction +## Overview + +Maximum likelihood estimation is a method for choosing the parameters of a +statistical model. + +The idea is simple: among all candidate parameter values, choose the one under +which the data we actually observed are most probable. + +Put differently, we treat the probability that the model assigns to the data as +a function of the parameters, and then maximize it. + +This principle is one of the foundations of modern statistics. + +It is used to fit distributions, as in this lecture, and also regression models, +macroeconomic models, and models in physics, biology and machine learning. + +(For example, training a classifier by minimizing cross-entropy loss is maximum +likelihood estimation.) + +One reason for its wide use is generality: once we have a model that assigns +probabilities to data, maximum likelihood tells us how to fit it, with no need +to invent a new method for each model. + +Another is that, under fairly general conditions, maximum likelihood estimates +converge to the true parameter values as the sample grows, and in large samples +no other method gives more precise estimates. + +The idea has a long history. + +Daniel Bernoulli and Carl Friedrich Gauss used versions of it in the late 18th +and early 19th centuries, and Gauss justified the method of least squares by +showing that, when errors are normally distributed, it picks out the most +probable parameter values. + +The method was developed systematically, and given its name, by Ronald Fisher +between 1912 and 1922 {cite}`fisher1922mathematical`. + +{cite:t}`aldrich1997fisher` and {cite:t}`stigler2007epic` tell the story. + +In this lecture we apply maximum likelihood to a policy problem: estimating the +revenue raised by a wealth tax. + +We also compare it with the method of moments, which we studied in +{doc}`fitting_distributions`, and see why maximum likelihood is the standard +choice when data have heavy tails. + + +## A wealth tax Consider a situation where a policymaker is trying to estimate how much revenue a proposed wealth tax will raise. From 3fb835cff5763fb18a661bd59ce0b518ae9f71e5 Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 29 Sep 2026 10:49:34 +1000 Subject: [PATCH 3/5] mle: rewrite as an elementary lecture on maximum likelihood Build up the method through Bernoulli (loan defaults), Poisson, normal and Pareto examples; add consistency and asymptotic normality simulations, a comparison with the method of moments for Pareto tails, numerical MLE for the gamma distribution, and keep the Student's t fit to stock returns. The wealth tax application moves to a separate lecture. Co-Authored-By: Claude Opus 5.5 --- lectures/mle.md | 1094 +++++++++++++++++++++++++++++------------------ 1 file changed, 669 insertions(+), 425 deletions(-) diff --git a/lectures/mle.md b/lectures/mle.md index af6ad3ae..c3eb4b65 100644 --- a/lectures/mle.md +++ b/lectures/mle.md @@ -13,25 +13,6 @@ kernelspec: # Maximum Likelihood Estimation -In addition to what's in Anaconda, this lecture will need the following libraries: - -```{code-cell} ipython3 -:tags: [hide-output] - -!pip install --upgrade yfinance -``` - -```{code-cell} ipython3 -from scipy.stats import lognorm, pareto, expon -from scipy.stats import t as student_t -from scipy.stats import kurtosis -import numpy as np -from scipy.integrate import quad -import matplotlib.pyplot as plt -import pandas as pd -import yfinance as yf -``` - ## Overview Maximum likelihood estimation is a method for choosing the parameters of a @@ -71,587 +52,722 @@ between 1912 and 1922 {cite}`fisher1922mathematical`. {cite:t}`aldrich1997fisher` and {cite:t}`stigler2007epic` tell the story. -In this lecture we apply maximum likelihood to a policy problem: estimating the -revenue raised by a wealth tax. +In this lecture we build up the idea through a sequence of examples, check by +simulation that maximum likelihood estimates converge to the truth, and compare +maximum likelihood with the method of moments from {doc}`fitting_distributions`. -We also compare it with the method of moments, which we studied in -{doc}`fitting_distributions`, and see why maximum likelihood is the standard -choice when data have heavy tails. +We also show how to compute maximum likelihood estimates numerically when no +formula is available, and apply this to heavy-tailed stock returns. +In addition to what's in Anaconda, this lecture will need the following libraries: -## A wealth tax - -Consider a situation where a policymaker is trying to estimate how much revenue -a proposed wealth tax will raise. +```{code-cell} ipython3 +:tags: [hide-output] -The proposed tax is +!pip install --upgrade yfinance +``` -$$ - h(w) = - \begin{cases} - a w & \text{if } w \leq \bar w \\ - a \bar{w} + b (w-\bar{w}) & \text{if } w > \bar w - \end{cases} -$$ +We will use the following imports. -where $w$ is wealth. +```{code-cell} ipython3 +import numpy as np +import matplotlib.pyplot as plt +import pandas as pd +import scipy.stats +from scipy.optimize import minimize +import yfinance as yf +``` -```{prf:example} -:label: mle_ex_wt -For example, if $a = 0.05$, $b = 0.1$, and $\bar w = 2.5$, this means +## The likelihood function -* a 5% tax on wealth up to 2.5 and -* a 10% tax on wealth in excess of 2.5. +We start with an example that we will meet again in {doc}`bayes_intro`. -The unit is 100,000, so $w= 2.5$ means 250,000 dollars. -``` -Let's go ahead and define $h$: +A bank makes a series of small loans, each of which either defaults or is +repaid. -```{code-cell} ipython3 -def h(w, a=0.05, b=0.1, w_bar=2.5): - if w <= w_bar: - return a * w - else: - return a * w_bar + b * (w - w_bar) -``` +We encode the outcome of loan $i$ as $Y_i = 1$ if it defaults and $Y_i = 0$ if +it is repaid. -For a population of size $N$, where individual $i$ has wealth $w_i$, total revenue raised by -the tax will be +We assume that the outcomes are IID Bernoulli draws with $$ - T = \sum_{i=1}^{N} h(w_i) +\mathbb{P}\{Y_i = 1\} = \theta +\quad \text{and} \quad +\mathbb{P}\{Y_i = 0\} = 1 - \theta $$ -We wish to calculate this quantity. +The default probability θ is unknown, and our job is to estimate it from the +observed outcomes $y_1, \ldots, y_n$. -The problem we face is that, in most countries, wealth is not observed for all individuals. +Both cases can be written as a single expression, since +$\mathbb{P}\{Y_i = y\} = \theta^y (1 - \theta)^{1 - y}$ for $y \in \{0, 1\}$. -Collecting and maintaining accurate wealth data for all individuals or households in a country -is just too hard. +By independence, the probability of observing the whole sample is therefore -So let's suppose instead that we obtain a sample $w_1, w_2, \cdots, w_n$ telling us the wealth of $n$ randomly selected individuals. +$$ +L(\theta) += \prod_{i=1}^n \theta^{y_i} (1 - \theta)^{1 - y_i} += \theta^k (1 - \theta)^{n - k} +$$ (eq:bern_lik) -For our exercise we are going to use a sample of $n = 10,000$ observations from wealth data in the US in 2016. +where $k = \sum_i y_i$ is the number of defaults. -```{code-cell} ipython3 -n = 10_000 -``` - -The data is derived from the -[Survey of Consumer Finances](https://en.wikipedia.org/wiki/Survey_of_Consumer_Finances) (SCF). +Viewed as a function of θ, with the data held fixed, $L$ is called the +**likelihood function**. +The **maximum likelihood estimate** $\hat \theta$ is the value of θ that +maximizes $L(\theta)$. -The following code imports this data and reads it into an array called `sample`. +Let's generate some data, using a true default probability of 0.15. ```{code-cell} ipython3 -:tags: [hide-input] - -url = 'https://raw.githubusercontent.com/QuantEcon/data-lectures/main/lectures/SCF_plus_mini_no_weights.csv' -df = pd.read_csv(url) -df = df.dropna() -df = df[df['year'] == 2016] -df = df.loc[df['n_wealth'] > 1 ] # restricting data to net worth > 1 -rv = df['n_wealth'].sample(n=n, random_state=1234) -rv = rv.to_numpy() / 100_000 -sample = rv +rng = np.random.default_rng(1234) +θ_true = 0.15 +n = 20 +y = (rng.random(n) < θ_true).astype(int) +k = y.sum() +y ``` -Let's histogram this sample. +Now we plot the likelihood function. ```{code-cell} ipython3 --- mystnb: figure: - caption: Histogram of US household wealth - name: fig:mle-wealth-hist + caption: Likelihood function for the default probability + name: fig:mle-bernoulli-lik --- +θ_grid = np.linspace(0.001, 0.999, 500) + +def bernoulli_likelihood(θ, k, n): + return θ**k * (1 - θ)**(n - k) + fig, ax = plt.subplots() -ax.set_xlim(-1, 20) -density, edges = np.histogram(sample, bins=5000, density=True) -prob = density * np.diff(edges) -ax.stairs(prob, edges, fill=True, alpha=0.8, label=r"unit: $\$100,000$") -ax.set_ylabel("prob") -ax.set_xlabel("net wealth") +ax.plot(θ_grid, bernoulli_likelihood(θ_grid, k, n), lw=2) +ax.axvline(k / n, color='k', ls='--', lw=1, label='sample mean') +ax.set_xlabel(r'$\theta$') +ax.set_ylabel('likelihood') ax.legend() plt.show() ``` -The histogram shows that many people have very low wealth and a few people have -very high wealth. +The likelihood peaks at the sample mean $k/n$, the fraction of loans that +defaulted. +This is the value of θ under which the data we observed are most probable. -We will take the full population size to be +It is important to be clear about what $L$ is. -```{code-cell} ipython3 -N = 100_000_000 -``` +For a fixed θ, the expression in {eq}`eq:bern_lik` is a probability, and +summing it over all possible data sets gives one. -How can we estimate total revenue from the full population using only the sample data? +As a function of θ with the data fixed, it is not a probability density over θ +and need not integrate to one. -Our plan is to assume that wealth of each individual is a draw from a distribution with density $f$. +That is why we give it a different name. -If we obtain an estimate of $f$ we can then approximate $T$ as follows: +Let's see what happens to the likelihood as we observe more loans. -$$ - T = \sum_{i=1}^{N} h(w_i) - = N \frac{1}{N} \sum_{i=1}^{N} h(w_i) - \approx N \int_{0}^{\infty} h(w)f(w) dw -$$ (eq:est_rev) +Since the likelihood gets very small when $n$ is large, we divide each curve by +its maximum value, so that all the curves peak at one. -(The sample mean should be close to the mean by the law of large numbers.) +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Likelihood functions for increasing sample sizes + name: fig:mle-bernoulli-lik-n +--- +fig, ax = plt.subplots() +for n in (10, 100, 1000): + y = (rng.random(n) < θ_true).astype(int) + k = y.sum() + L = bernoulli_likelihood(θ_grid, k, n) + ax.plot(θ_grid, L / L.max(), lw=2, label=f'$n = {n}$') +ax.axvline(θ_true, color='k', ls=':', lw=1, label=r'true $\theta$') +ax.set_xlim(0, 0.5) +ax.set_xlabel(r'$\theta$') +ax.set_ylabel('likelihood (normalized)') +ax.legend() +plt.show() +``` -The problem now is: how do we estimate $f$? +As $n$ grows, the likelihood concentrates around the true value. +With more data, fewer values of θ are consistent with what we observe. -## Maximum likelihood estimation +This is a first sign of the consistency of maximum likelihood, which we study +{ref}`below `. -[Maximum likelihood estimation](https://en.wikipedia.org/wiki/Maximum_likelihood_estimation) -is a method of estimating an unknown distribution. -Maximum likelihood estimation has two steps: +## The log likelihood -1. Guess what the underlying distribution is (e.g., normal with mean $\mu$ and - standard deviation $\sigma$). -2. Estimate the parameter values (e.g., estimate $\mu$ and $\sigma$ for the - normal distribution) +In practice we almost always work with the logarithm of the likelihood, +called the **log likelihood**, rather than the likelihood itself. -These are the same two steps discussed in {doc}`fitting_distributions`, which -chooses parameters by the method of moments instead, and which also shows how -to judge the resulting fit. +There are two reasons. -One possible assumption for the wealth is that each -$w_i$ is [log-normally distributed](https://en.wikipedia.org/wiki/Log-normal_distribution), -with parameters $\mu \in (-\infty,\infty)$ and $\sigma \in (0,\infty)$. +The first is that the logarithm turns products into sums, which are easier to +differentiate. -(This means that $\ln w_i$ is normally distributed with mean $\mu$ and standard deviation $\sigma$.) +The second is numerical. -You can see that this assumption is not completely unreasonable because, if we -histogram log wealth instead of wealth, the picture starts to look something -like a bell-shaped curve. +The likelihood is a product of many numbers less than one, and for large +samples it becomes too small to represent on a computer. ```{code-cell} ipython3 ---- -mystnb: - figure: - caption: Histogram of log wealth - name: fig:mle-log-wealth-hist ---- -ln_sample = np.log(sample) -fig, ax = plt.subplots() -ax.hist(ln_sample, density=True, bins=200, histtype='stepfilled', alpha=0.8) -ax.set_xlabel("log of net wealth") -ax.set_ylabel("density") -plt.show() +n = 5_000 +y = (rng.random(n) < θ_true).astype(int) +k = y.sum() +bernoulli_likelihood(θ_true, k, n) ``` -Now our job is to obtain the maximum likelihood estimates of $\mu$ and $\sigma$, which -we denote by $\hat{\mu}$ and $\hat{\sigma}$. +The true value is positive, but it is smaller than the smallest positive +floating point number, so the computer rounds it to zero. + +The log likelihood has no such problem. -These estimates can be found by maximizing the likelihood function given the -data. +```{code-cell} ipython3 +def bernoulli_log_likelihood(θ, k, n): + return k * np.log(θ) + (n - k) * np.log(1 - θ) + +bernoulli_log_likelihood(θ_true, k, n) +``` + +Since the logarithm is strictly increasing, maximizing the log likelihood gives +the same answer as maximizing the likelihood. -The pdf of a lognormally distributed random variable $X$ is given by: +For the loan data, the log likelihood is $$ - f(x; \mu, \sigma) - = \frac{1}{x \sigma \sqrt{2\pi}} - \exp \left( - \frac{(\ln x - \mu)^2}{2 \sigma^2} \right) +\ell(\theta) = \ln L(\theta) = k \ln \theta + (n - k) \ln (1 - \theta) $$ -For our sample $w_1, w_2, \cdots, w_n$, the [likelihood function](https://en.wikipedia.org/wiki/Likelihood_function) is given by +Setting the derivative to zero gives $$ - L(\mu, \sigma) = \prod_{i=1}^{n} f(w_i; \mu, \sigma) +\frac{k}{\theta} - \frac{n - k}{1 - \theta} = 0 +\quad \implies \quad +\hat \theta = \frac{k}{n} $$ -The likelihood function can be viewed as both +So the maximum likelihood estimate is the sample mean, confirming what we saw +in the figure. -* the joint distribution of the sample (which is assumed to be IID) and -* the "likelihood" of parameters $(\mu, \sigma)$ given the data. -Taking logs on both sides gives us the log likelihood function, which is +## More examples + +In this section we derive maximum likelihood estimates for three more +distributions. + +In each case the data $x_1, \ldots, x_n$ are assumed to be IID draws from a +distribution with density or probability mass function $f(x; \theta)$. + +The log likelihood is then $$ -\begin{aligned} - \ell(\mu, \sigma) - & = \ln \left[ \prod_{i=1}^{n} f(w_i; \mu, \sigma) \right] \\ - & = -\sum_{i=1}^{n} \ln w_i - - \frac{n}{2} \ln(2\pi) - \frac{n}{2} \ln \sigma^2 - \frac{1}{2\sigma^2} - \sum_{i=1}^n (\ln w_i - \mu)^2 -\end{aligned} +\ell(\theta) = \sum_{i=1}^n \ln f(x_i; \theta) $$ -To find where this function is maximized we find its partial derivatives with respect to $\mu$ and $\sigma^2$ and set them to zero. +### Poisson + +In {doc}`fitting_distributions` we modeled the number of goals scored in a +football match with a Poisson distribution. -Let's first find the maximum likelihood estimate (MLE) of $\mu$ +The Poisson distribution with parameter λ has probability mass function +$f(x; \lambda) = e^{-\lambda} \lambda^x / x!$, so $$ -\begin{aligned} -\frac{\partial \ell}{\partial \mu} - = \frac{1}{\sigma^2} \sum_{i=1}^n (\ln w_i - \mu) = 0 - & \implies \sum_{i=1}^n \ln w_i - n \mu = 0 \\ - & \implies \hat{\mu} = \frac{1}{n} \sum_{i=1}^n \ln w_i -\end{aligned} +\ell(\lambda) = -n \lambda + \ln \lambda \sum_{i=1}^n x_i - \sum_{i=1}^n \ln x_i! $$ -Now let's find the MLE of $\sigma$ +Setting the derivative to zero gives $$ -\begin{aligned} -\frac{\partial \ell}{\partial \sigma^2} - = - \frac{n}{2\sigma^2} + \frac{1}{2\sigma^4} - \sum_{i=1}^n (\ln w_i - \mu)^2 = 0 - & \implies \frac{n}{2\sigma^2} = - \frac{1}{2\sigma^4} \sum_{i=1}^n (\ln w_i - \mu)^2 \\ - & \implies \hat{\sigma} = - \left( \frac{1}{n} \sum_{i=1}^{n}(\ln w_i - \hat{\mu})^2 \right)^{1/2} -\end{aligned} +-n + \frac{1}{\lambda} \sum_{i=1}^n x_i = 0 +\quad \implies \quad +\hat \lambda = \frac{1}{n} \sum_{i=1}^n x_i $$ -Now that we have derived the expressions for $\hat{\mu}$ and $\hat{\sigma}$, -let's compute them for our wealth sample. +Let's check this on the goals data. ```{code-cell} ipython3 -μ_hat = np.mean(ln_sample) -μ_hat +url = ('https://github.com/QuantEcon/data-lectures/raw/main/' + 'lectures/epl_match_goals.csv') +matches = pd.read_csv(url) +goals = matches['home_goals'] + matches['away_goals'] ``` -```{code-cell} ipython3 -num = (ln_sample - μ_hat)**2 -σ_hat = (np.mean(num))**(1/2) -σ_hat -``` - -Let's plot the lognormal pdf using the estimated parameters against our sample data. +We compute the log likelihood on a grid of values for λ, using the log +probability mass function provided by SciPy. ```{code-cell} ipython3 --- mystnb: figure: - caption: Lognormal fit to wealth - name: fig:mle-lognormal-fit + caption: Log likelihood for the Poisson parameter + name: fig:mle-poisson-loglik --- -dist_lognorm = lognorm(σ_hat, scale=np.exp(μ_hat)) -x = np.linspace(0, 50, 10000) +λ_grid = np.linspace(2, 3.5, 200) +ℓ = [scipy.stats.poisson.logpmf(goals, λ).sum() for λ in λ_grid] fig, ax = plt.subplots() -ax.set_xlim(-1, 20) - -ax.hist(sample, density=True, bins=5_000, histtype='stepfilled', alpha=0.5) -ax.plot(x, dist_lognorm.pdf(x), 'k-', lw=2, label='lognormal pdf') -ax.set_xlabel("net wealth") -ax.set_ylabel("density") +ax.plot(λ_grid, ℓ, lw=2) +ax.axvline(goals.mean(), color='k', ls='--', lw=1, label='sample mean') +ax.set_xlabel(r'$\lambda$') +ax.set_ylabel('log likelihood') ax.legend() plt.show() ``` -Our estimated lognormal distribution appears to be a reasonable fit for the overall data. +The maximum is at the sample mean, as the formula predicts. +```{code-cell} ipython3 +λ_grid[np.argmax(ℓ)], goals.mean() +``` -### Comparison with the method of moments +(The two numbers differ slightly only because the grid is finite.) -Notice that $\hat{\mu}$ and $\hat{\sigma}$ are just the sample mean and sample -standard deviation of log wealth. +### Normal -So, for the lognormal class, maximum likelihood gives the same answer as -fitting a normal distribution to log wealth by the method of moments. +The normal distribution has two parameters, μ and σ, and density -This is *not* the same as the lognormal fit in {doc}`fitting_distributions`, -which matches the mean and variance of wealth itself, rather than of log -wealth. +$$ +f(x; \mu, \sigma) += \frac{1}{\sigma \sqrt{2 \pi}} + \exp \left( - \frac{(x - \mu)^2}{2 \sigma^2} \right) +$$ -Let's compute that fit too. +The log likelihood is -```{code-cell} ipython3 -def fit_lognormal_mom(sample): - "Lognormal fit matching the sample mean and variance of the data." - m, v = sample.mean(), sample.var() - σ_squared = np.log(1 + v / m**2) - μ = np.log(m) - σ_squared / 2 - return lognorm(s=np.sqrt(σ_squared), scale=np.exp(μ)) +$$ +\ell(\mu, \sigma) += - \frac{n}{2} \ln (2 \pi) - n \ln \sigma + - \frac{1}{2 \sigma^2} \sum_{i=1}^n (x_i - \mu)^2 +$$ (eq:normal_loglik) -dist_lognorm_mom = fit_lognormal_mom(sample) -``` +Setting the partial derivatives with respect to μ and σ to zero gives -To compare the two fits we use the Kolmogorov-Smirnov statistic from -{doc}`fitting_distributions`, which is the largest vertical gap between the -ECDF of the data and the fitted CDF. +$$ +\hat \mu = \frac{1}{n} \sum_{i=1}^n x_i +\quad \text{and} \quad +\hat \sigma = \left( \frac{1}{n} \sum_{i=1}^n (x_i - \hat \mu)^2 \right)^{1/2} +$$ -```{code-cell} ipython3 -def ks_statistic(sample, u): - "Largest vertical distance between the ECDF of the sample and the CDF of u." - x_sorted = np.sort(sample) - n = len(x_sorted) - F = u.cdf(x_sorted) - above = np.arange(1, n+1) / n - F # gap just after each jump - below = F - np.arange(0, n) / n # gap just before each jump - return max(above.max(), below.max()) -``` +Let's look at the log likelihood for the heights of US adult women, which we +fitted with a normal distribution in {doc}`fitting_distributions`. ```{code-cell} ipython3 -pd.DataFrame({ - 'σ': [σ_hat, dist_lognorm_mom.kwds['s']], - 'KS statistic': [ks_statistic(sample, dist_lognorm), - ks_statistic(sample, dist_lognorm_mom)]}, - index=['maximum likelihood', 'method of moments']) +url = ('https://github.com/QuantEcon/data-lectures/raw/main/' + 'lectures/us_adult_heights.csv') +heights = pd.read_csv(url) +female = heights[heights['sex'] == 'female']['height_cm'].to_numpy() ``` -The two methods give quite different parameters, and maximum likelihood gives -a much closer fit. +With two parameters, the log likelihood is a surface, which we display with a +contour plot. -The reason is that the sample variance of wealth is dominated by a handful of -very wealthy households, so matching it pulls the whole fitted distribution -toward them. +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Log likelihood for the normal distribution + name: fig:mle-normal-contour +--- +μ_hat, σ_hat = female.mean(), female.std() -Maximum likelihood works with log wealth, where these households are far less -extreme. +μ_grid = np.linspace(μ_hat - 1.5, μ_hat + 1.5, 200) +σ_grid = np.linspace(σ_hat - 1.2, σ_hat + 1.2, 200) +M, S = np.meshgrid(μ_grid, σ_grid) -We now use {eq}`eq:est_rev` and the maximum likelihood fit to calculate total -revenue. +def normal_log_likelihood(μ, σ, x): + n = len(x) + return (- n / 2 * np.log(2 * np.pi) - n * np.log(σ) + - ((x - μ)**2).sum() / (2 * σ**2)) -We will compute the integral using numerical integration via SciPy's -[quad](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) -function +Z = np.vectorize(lambda μ, σ: normal_log_likelihood(μ, σ, female))(M, S) -```{code-cell} ipython3 -def total_revenue(dist): - integral, _ = quad(lambda x: h(x) * dist.pdf(x), 0, 100_000) - T = N * integral - return T +fig, ax = plt.subplots() +cs = ax.contour(M, S, Z, levels=20, cmap='viridis') +ax.plot(μ_hat, σ_hat, 'ko', label='maximum') +ax.set_xlabel(r'$\mu$') +ax.set_ylabel(r'$\sigma$') +ax.legend() +plt.show() ``` +The maximum sits at the sample mean and standard deviation. + ```{code-cell} ipython3 -tr_lognorm = total_revenue(dist_lognorm) -tr_lognorm +μ_hat, σ_hat ``` -(Our unit was 100,000 dollars, so this means that actual revenue is 100,000 -times as large.) +Notice that $\hat \sigma$ divides by $n$ rather than $n - 1$. +This is also what NumPy computes by default: `female.std()` uses $n$, while +`female.std(ddof=1)` uses $n - 1$. -## Pareto distribution +For large samples the difference is negligible, but for small samples it is +not, as you are asked to check in {ref}`an exercise `. -We mentioned above that using maximum likelihood estimation requires us to make -a prior assumption of the underlying distribution. +The normal log likelihood also sheds light on {doc}`simple_linear_regression`. -Previously we assumed that the distribution is lognormal. +Suppose that $y_i = \alpha + \beta x_i + \epsilon_i$, where the errors +$\epsilon_i$ are IID and normal with mean zero and standard deviation σ. -Suppose instead we assume that $w_i$ are drawn from the -[Pareto Distribution](https://en.wikipedia.org/wiki/Pareto_distribution) -with parameters $b$ and $x_m$. +Then, replacing $x_i - \mu$ by $y_i - \alpha - \beta x_i$ in +{eq}`eq:normal_loglik`, the log likelihood is -In this case, the maximum likelihood estimates are known to be +$$ +\ell(\alpha, \beta, \sigma) += - \frac{n}{2} \ln (2 \pi) - n \ln \sigma + - \frac{1}{2 \sigma^2} \sum_{i=1}^n (y_i - \alpha - \beta x_i)^2 +$$ + +For any σ, maximizing this over α and β is the same as minimizing the sum of +squared residuals. + +So, with normal errors, ordinary least squares *is* maximum likelihood. + +### Pareto + +Our last example is the Pareto distribution from {doc}`heavy_tails`. + +We assume that the minimum value $\bar x$ is known, so that the only unknown +parameter is the tail index α. + +The density is $f(x; \alpha) = \alpha \bar x^\alpha / x^{\alpha + 1}$ for +$x \geq \bar x$, so $$ - \hat{b} = \frac{n}{\sum_{i=1}^{n} \ln (w_i/\hat{x}_m)} - \quad \text{and} \quad - \hat{x}_m = \min_{i} w_i +\ell(\alpha) = n \ln \alpha + n \alpha \ln \bar x + - (\alpha + 1) \sum_{i=1}^n \ln x_i $$ -Let's calculate them. +Setting the derivative to zero gives -```{code-cell} ipython3 -xm_hat = min(sample) -xm_hat -``` +$$ +\frac{n}{\alpha} - \sum_{i=1}^n \ln \frac{x_i}{\bar x} = 0 +\quad \implies \quad +\hat \alpha = \frac{n}{\sum_{i=1}^n \ln (x_i / \bar x)} +$$ (eq:pareto_mle) + +Here's a function that computes it. ```{code-cell} ipython3 -den = np.log(sample/xm_hat) -b_hat = 1/np.mean(den) -b_hat +def pareto_mle(x, x_bar=1.0): + return len(x) / np.sum(np.log(x / x_bar)) ``` -Now let's recompute total revenue. +For the Bernoulli, Poisson and normal distributions, the maximum likelihood +estimates coincide with the method of moments estimates from +{doc}`fitting_distributions`, since both are built from the sample mean and +variance. + +For the Pareto distribution this is not the case. + +We will see {ref}`below ` that the difference matters. + + +(mle_consistency)= +## Consistency + +An estimator is called **consistent** if it converges to the true parameter +value as the sample size grows. + +Let's check this by simulation for the Pareto estimate {eq}`eq:pareto_mle`. + +We set $\bar x = 1$ and α = 1.5, and draw many samples of size $n$. + +To draw from the Pareto distribution we use the inverse transform method: if +$U$ is uniform on $(0, 1)$, then $\bar x U^{-1/\alpha}$ is Pareto with minimum +$\bar x$ and tail index α. ```{code-cell} ipython3 -dist_pareto = pareto(b=b_hat, scale=xm_hat) -tr_pareto = total_revenue(dist_pareto) -tr_pareto +def pareto_draws(α, size, rng, x_bar=1.0): + return x_bar * rng.random(size)**(-1 / α) ``` -The number is very different! +For each sample size, we compute the estimate on 2,000 independent samples and +plot a kernel density estimate (see {doc}`observed_distributions`) of the +results. ```{code-cell} ipython3 -tr_pareto / tr_lognorm +--- +mystnb: + figure: + caption: Maximum likelihood estimates of the tail index + name: fig:mle-consistency +--- +α_true = 1.5 +num_reps = 2_000 +sample_sizes = (50, 200, 1_000, 5_000) + +α_grid = np.linspace(0.9, 2.3, 400) +fig, ax = plt.subplots() +for n in sample_sizes: + x = pareto_draws(α_true, (num_reps, n), rng) + α_hats = n / np.sum(np.log(x), axis=1) + kde = scipy.stats.gaussian_kde(α_hats) + ax.plot(α_grid, kde(α_grid), lw=2, label=f'$n = {n}$') +ax.axvline(α_true, color='k', ls=':', lw=1, label=r'true $\alpha$') +ax.set_xlabel(r'estimate of $\alpha$') +ax.set_ylabel('density') +ax.legend() +plt.show() ``` -We see that choosing the right distribution is extremely important. +As $n$ grows, the estimates concentrate on the true value. +How fast do they concentrate? +The {doc}`central limit theorem ` suggests that the spread should +shrink like $1/\sqrt{n}$. -Let's compare the fitted Pareto distribution to the histogram: +For this estimator it is known that, for large $n$, $\hat \alpha$ is +approximately normal with mean α and standard deviation $\alpha / \sqrt{n}$. + +Let's check by standardizing the estimates and comparing them with a standard +normal density. ```{code-cell} ipython3 --- mystnb: figure: - caption: Pareto fit to wealth - name: fig:mle-pareto-fit + caption: Standardized estimates and the standard normal + name: fig:mle-asymptotic-normal --- -fig, ax = plt.subplots() -ax.set_xlim(-1, 20) -ax.set_ylim(0, 1.75) +n = 1_000 +x = pareto_draws(α_true, (num_reps, n), rng) +α_hats = n / np.sum(np.log(x), axis=1) +z = (α_hats - α_true) / (α_true / np.sqrt(n)) -ax.hist(sample, density=True, bins=5_000, histtype='stepfilled', alpha=0.5) -ax.plot(x, dist_pareto.pdf(x), 'k-', lw=2, label='Pareto pdf') -ax.set_xlabel("net wealth") -ax.set_ylabel("density") +z_grid = np.linspace(-4, 4, 200) +fig, ax = plt.subplots() +ax.hist(z, bins=50, density=True, alpha=0.4, label='standardized estimates') +ax.plot(z_grid, scipy.stats.norm.pdf(z_grid), lw=2, label='standard normal') +ax.set_xlabel('standardized estimate') +ax.set_ylabel('density') ax.legend() - plt.show() ``` -The fit is poor, and the KS statistic confirms it. +The fit is close. -```{code-cell} ipython3 -ks_statistic(sample, dist_pareto), ks_statistic(sample, dist_lognorm) -``` +These two properties, consistency and approximate normality in large samples, +hold for maximum likelihood estimates in general, under regularity conditions +that are satisfied in most applications. -For the sample as a whole, the lognormal distribution is clearly preferred. +Moreover, the spread of the approximating normal distribution is as small as +that of any reasonable estimator, which is the sense in which maximum +likelihood makes the best use of the data. +Proofs can be found in any graduate text on statistics. -## Fitting the right tail -There is no "best" distribution --- every choice we make is an assumption. +(mle_vs_mom)= +## Maximum likelihood and the method of moments -The lognormal distribution fits the sample as a whole, but for questions about -the richest households, what matters is the fit in the upper tail. +Let's now compare maximum likelihood with the method of moments for the Pareto +distribution. -To study the tail, let's set a minimum threshold of net worth in our dataset. +When α > 1, the Pareto distribution has mean $\alpha \bar x / (\alpha - 1)$. -We set an arbitrary threshold of $500,000 and read the data into `sample_tail`. +Setting this equal to the sample mean $m_n = \frac{1}{n} \sum_i x_i$ and +solving for α gives the method of moments estimate -```{code-cell} ipython3 -:tags: [hide-input] +$$ +\tilde \alpha = \frac{m_n}{m_n - \bar x} +$$ -df_tail = df.loc[df['n_wealth'] > 500_000 ] -rv_tail = df_tail['n_wealth'].sample(n=10_000, random_state=4321) -rv_tail = rv_tail.to_numpy() -sample_tail = rv_tail/500_000 +```{code-cell} ipython3 +def pareto_mom(x, x_bar=1.0): + m = np.mean(x, axis=-1) + return m / (m - x_bar) ``` -Let's fit a lognormal and a Pareto distribution to this data by maximum -likelihood, using the formulas derived above. +Since every observation is at least $\bar x$, we have $m_n > \bar x$ and hence +$\tilde \alpha > 1$ for *every* data set. -```{code-cell} ipython3 -ln_sample_tail = np.log(sample_tail) -μ_hat_tail = np.mean(ln_sample_tail) -σ_hat_tail = np.std(ln_sample_tail) -dist_lognorm_tail = lognorm(σ_hat_tail, scale=np.exp(μ_hat_tail)) +So if the true α is less than or equal to one, the method of moments can never +recover it, however much data we have. -xm_hat_tail = min(sample_tail) -b_hat_tail = 1/np.mean(np.log(sample_tail/xm_hat_tail)) -dist_pareto_tail = pareto(b=b_hat_tail, scale=xm_hat_tail) -``` +The underlying problem is that when α ≤ 1 the mean is infinite, so there is no +population moment for the sample mean to match. -Since we care about the tail, we compare the fits with Q-Q plots, which were -introduced in {doc}`fitting_distributions` and show clearly how a fit behaves -at the extremes. +Even when α is a little above one, so that the mean exists, the variance is +infinite whenever α ≤ 2. -Tail wealth spans several orders of magnitude, so we draw the plots on log -scales. +Then, as we saw in {doc}`heavy_tails`, the sample mean settles down only +slowly, and the method of moments estimate inherits this instability. -```{code-cell} ipython3 -def qq_plot(sample, u, ax, **kwargs): - "Plot sample quantiles against the quantiles of the distribution u." - x_sorted = np.sort(sample) - n = len(x_sorted) - p = (np.arange(1, n+1) - 0.5) / n - ax.plot(u.ppf(p), x_sorted, '.', ms=3, alpha=0.6, **kwargs) - lo, hi = u.ppf(p[0]), u.ppf(p[-1]) - ax.plot([lo, hi], [lo, hi], 'k--', lw=2) - ax.set_xlabel('fitted quantiles') - ax.set_ylabel('sample quantiles') -``` +Maximum likelihood suffers from neither problem, since it works with the +density, which is well defined for every α > 0. + +Let's compare the two estimators by simulation, for α = 0.8 and α = 1.5, with +samples of size 1,000. ```{code-cell} ipython3 --- mystnb: figure: - caption: Q-Q plots for tail wealth - name: fig:mle-tail-qq + caption: Maximum likelihood and method of moments estimates + name: fig:mle-vs-mom --- -fig, axes = plt.subplots(1, 2, figsize=(10, 4.5)) -for ax, u, label in zip(axes, - (dist_lognorm_tail, dist_pareto_tail), - ('lognormal', 'Pareto')): - qq_plot(sample_tail, u, ax) - ax.set_xscale('log') - ax.set_yscale('log') - ax.set_title(label) +n = 1_000 +fig, axes = plt.subplots(1, 2, figsize=(10, 4)) +for ax, α, lims in zip(axes, (0.8, 1.5), ((0.65, 1.2), (1.0, 2.0))): + x = pareto_draws(α, (num_reps, n), rng) + α_mle = n / np.sum(np.log(x), axis=1) + α_mom = pareto_mom(x) + bins = np.linspace(*lims, 80) + ax.hist(α_mle, bins=bins, density=True, alpha=0.5, + label='maximum likelihood') + ax.hist(α_mom, bins=bins, density=True, alpha=0.5, + label='method of moments') + ax.axvline(α, color='k', ls=':', lw=2, label=r'true $\alpha$') + ax.set_title(rf'$\alpha = {α}$') + ax.set_xlabel(r'estimate of $\alpha$') + ax.set_ylabel('density') + ax.legend() plt.tight_layout() plt.show() ``` -The lognormal fit bends away from the 45 degree line: it puts probability on -wealth below the threshold, which cannot occur in this sample, and it -underpredicts large fortunes. +When α = 0.8 (left), the method of moments estimates all lie above one, far +from the truth, while the maximum likelihood estimates are centered on it. + +When α = 1.5 (right), both estimators are concentrated near the truth, but +the method of moments estimates are much more dispersed, with a long left tail +caused by samples whose mean is inflated by a few very large draws. + +Here are the standard deviations of the two sets of estimates for α = 1.5. + +```{code-cell} ipython3 +np.std(α_mle), np.std(α_mom) +``` + +Heavy tails are common in economic data, as we saw in {doc}`heavy_tails`, and +this is one reason why maximum likelihood is the standard method for fitting +heavy-tailed distributions. + -The Pareto fit tracks the data closely until far into the tail, up to around 20 -times the threshold, or 10 million dollars. +## Numerical maximum likelihood -Beyond that point the data flatten out, because no household in our data set -is worth more than about 50 million dollars. +In all the examples so far, we could find the maximizer with pencil and paper. -This reflects the survey rather than the population, since the SCF does not -capture the very wealthiest households. +For many distributions this is not possible, and we have to maximize the log +likelihood numerically. -The KS statistic also favors the Pareto distribution. +Let's see how this works for the gamma distribution, which we fitted to house +prices by the method of moments in {doc}`fitting_distributions`. ```{code-cell} ipython3 -ks_statistic(sample_tail, dist_lognorm_tail), ks_statistic(sample_tail, dist_pareto_tail) +url = ('https://github.com/QuantEcon/data-lectures/raw/main/' + 'lectures/ames_house_prices.csv') +houses = pd.read_csv(url) +price = houses['price'].to_numpy() / 1_000 # price in thousands of dollars ``` -So, while the lognormal distribution was a good fit for the entire dataset, the -Pareto distribution is a better fit for the right hand tail. +The gamma distribution has a shape parameter $a$ and a scale parameter $s$, +both positive. +There is no closed-form expression for the maximum likelihood estimate of $a$, +so we will compute it numerically. -### Why not the method of moments? +Most numerical optimization routines minimize rather than maximize, so we +minimize the *negative* of the log likelihood. -Here is the estimated Pareto parameter for the tail. +To keep both parameters positive, we optimize over their logarithms, which can +take any real value. ```{code-cell} ipython3 -b_hat_tail +def gamma_neg_log_likelihood(params, x): + a, s = np.exp(params) + return -np.sum(scipy.stats.gamma.logpdf(x, a, scale=s)) ``` -This value is below 1, which has a striking implication. +We need an initial guess, and the method of moments estimates from +{doc}`fitting_distributions` are a natural choice. -A Pareto distribution with parameter $b$ has mean $b x_m / (b - 1)$ when -$b > 1$, and infinite mean when $b \leq 1$. +```{code-cell} ipython3 +m, v = price.mean(), price.var() +a_mom, s_mom = m**2 / v, v / m + +result = minimize(gamma_neg_log_likelihood, + x0=np.log([a_mom, s_mom]), + args=(price,)) +a_mle, s_mle = np.exp(result.x) +a_mle, s_mle +``` -So the fitted distribution has no finite mean, let alone a finite variance. +SciPy's distributions have a `fit` method that does the same computation, which +we can use as a check. -(We discuss distributions of this kind in {doc}`heavy_tails`.) +(Setting `floc=0` tells `fit` not to add a location parameter.) -This is a problem for the method of moments. +```{code-cell} ipython3 +a_fit, _, s_fit = scipy.stats.gamma.fit(price, floc=0) +a_fit, s_fit +``` -With $x_m$ fixed, matching the population mean to the sample mean $\bar w$ gives +The answers agree. -$$ -\hat{b} = \frac{\bar w}{\bar w - x_m} -$$ +Here are the two fitted densities. -which is greater than 1 for *every* data set, since $\bar w > x_m$. +```{code-cell} ipython3 +--- +mystnb: + figure: + caption: Gamma densities fitted to house prices + name: fig:mle-gamma-prices +--- +u_mom = scipy.stats.gamma(a_mom, scale=s_mom) +u_mle = scipy.stats.gamma(a_mle, scale=s_mle) -The method of moments therefore cannot recover a value of $b$ at or below 1, -however much data we have. +x_grid = np.linspace(0, price.max(), 400) +fig, ax = plt.subplots() +ax.hist(price, bins=50, density=True, alpha=0.25) +ax.plot(x_grid, u_mom.pdf(x_grid), lw=2, label='method of moments') +ax.plot(x_grid, u_mle.pdf(x_grid), lw=2, label='maximum likelihood') +ax.set_xlabel('sale price (thousands of US$)') +ax.set_ylabel('density') +ax.legend() +plt.show() +``` -Let's see what it gives here. +To compare the fits we use the Kolmogorov-Smirnov statistic from +{doc}`fitting_distributions`, the largest vertical gap between the ECDF of the +data and the fitted CDF. ```{code-cell} ipython3 -w_bar_tail = np.mean(sample_tail) -w_bar_tail / (w_bar_tail - xm_hat_tail) +def ks_statistic(sample, u): + "Largest vertical distance between the ECDF of the sample and the CDF of u." + x_sorted = np.sort(sample) + n = len(x_sorted) + F = u.cdf(x_sorted) + above = np.arange(1, n+1) / n - F # gap just after each jump + below = F - np.arange(0, n) / n # gap just before each jump + return max(above.max(), below.max()) + +ks_statistic(price, u_mom), ks_statistic(price, u_mle) ``` -The trouble is that the method of moments assumes the moments it matches -exist. +Maximum likelihood gives a slightly closer fit. -Maximum likelihood makes no such assumption, since it works with the density -itself, which is well defined for every value of $b$. +Here the two methods give similar answers, because the house price data do not +have especially heavy tails. -This is one reason why maximum likelihood is the standard method for fitting -heavy-tailed distributions. +Remember, though, that maximum likelihood finds the best member of the family +we chose, not the best family. + +In {doc}`fitting_distributions` the lognormal family fitted these data better +than the gamma family, and no choice of gamma parameters can change that. + +Choosing the family is a separate step, for which the Q-Q plots and KS +statistics of {doc}`fitting_distributions` remain the right tools. ## Heavy-tailed returns @@ -659,10 +775,10 @@ heavy-tailed distributions. In {doc}`fitting_distributions` we fitted a Student's t distribution to monthly returns on Amazon shares, using the method of moments. -We noted there that the fitted value of the degrees of freedom $\nu$ was -unreliable, because it came from the sample kurtosis, a fourth moment. +We noted there that the fitted degrees of freedom ν were unreliable, because +they came from the sample kurtosis, a fourth moment. -Let's fit the same class by maximum likelihood instead. +Let's fit the same family by maximum likelihood instead. ```{code-cell} ipython3 :tags: [hide-output] @@ -672,39 +788,52 @@ prices = data['Close']['AMZN'] returns = prices.pct_change().dropna() * 100 ``` -Unlike the lognormal and Pareto cases, there is no closed-form expression for -the maximum likelihood estimates of the Student's t distribution. - -Instead, the `fit` method in SciPy maximizes the log likelihood numerically. +As for the gamma distribution, there is no closed-form expression for the +estimates, so we use the `fit` method, which maximizes the log likelihood +numerically. ```{code-cell} ipython3 -ν_mle, loc_mle, scale_mle = student_t.fit(returns) -dist_t_mle = student_t(df=ν_mle, loc=loc_mle, scale=scale_mle) +ν_mle, loc_mle, scale_mle = scipy.stats.t.fit(returns) +u_t_mle = scipy.stats.t(df=ν_mle, loc=loc_mle, scale=scale_mle) ν_mle ``` -For comparison, here is the method of moments fit from {doc}`fitting_distributions`. +For comparison, here is the method of moments fit from +{doc}`fitting_distributions`. ```{code-cell} ipython3 def fit_t_mom(sample): m, s = sample.mean(), sample.std() - ν = 4 + 6 / kurtosis(sample) - return student_t(df=ν, loc=m, scale=s * np.sqrt((ν - 2) / ν)) + ν = 4 + 6 / scipy.stats.kurtosis(sample) + return scipy.stats.t(df=ν, loc=m, scale=s * np.sqrt((ν - 2) / ν)) -dist_t_mom = fit_t_mom(returns) -dist_t_mom.kwds['df'] +u_t_mom = fit_t_mom(returns) +u_t_mom.kwds['df'] ``` -Maximum likelihood gives a noticeably smaller $\nu$, and hence heavier tails. +Maximum likelihood gives a noticeably smaller ν, and hence heavier tails. It also gives a closer fit. ```{code-cell} ipython3 -ks_statistic(returns, dist_t_mom), ks_statistic(returns, dist_t_mle) +ks_statistic(returns, u_t_mom), ks_statistic(returns, u_t_mle) ``` -The Q-Q plot confirms that the maximum likelihood fit tracks the extreme -returns well. +The Q-Q plot from {doc}`fitting_distributions` confirms that the maximum +likelihood fit tracks the extreme returns well. + +```{code-cell} ipython3 +def qq_plot(sample, u, ax): + "Plot sample quantiles against the quantiles of the distribution u." + x_sorted = np.sort(sample) + n = len(x_sorted) + p = (np.arange(1, n+1) - 0.5) / n + ax.plot(u.ppf(p), x_sorted, '.', ms=3, alpha=0.6) + lo, hi = u.ppf(p[0]), u.ppf(p[-1]) + ax.plot([lo, hi], [lo, hi], 'k--', lw=2) + ax.set_xlabel('fitted quantiles') + ax.set_ylabel('sample quantiles') +``` ```{code-cell} ipython3 --- @@ -714,11 +843,11 @@ mystnb: name: fig:mle-qq-returns-t --- fig, ax = plt.subplots() -qq_plot(returns, dist_t_mle, ax) +qq_plot(returns, u_t_mle, ax) plt.show() ``` -The lesson is the same as for the wealth tail: when the tails are heavy, the +The lesson is the same as for the Pareto simulations: when tails are heavy, sample moments are dominated by a few extreme observations, and maximum likelihood makes better use of the data. @@ -726,66 +855,181 @@ likelihood makes better use of the data. ## Exercises ```{exercise-start} -:label: mle_ex1 +:label: mle_ex_exp +``` + +The exponential distribution with rate λ > 0 has density +$f(x; \lambda) = \lambda e^{-\lambda x}$ for $x \geq 0$. + +1. Derive the maximum likelihood estimate of λ from IID observations + $x_1, \ldots, x_n$. +2. Draw a sample of size 500 from the exponential distribution with λ = 2 and + check that your formula agrees with `scipy.stats.expon.fit(x, floc=0)`. + (SciPy parameterizes the exponential distribution by the scale $1/\lambda$.) + +```{exercise-end} +``` + +```{solution-start} mle_ex_exp +:class: dropdown ``` -Suppose we assume wealth is [exponentially](https://en.wikipedia.org/wiki/Exponential_distribution) -distributed with parameter $\lambda > 0$. -The maximum likelihood estimate of $\lambda$ is given by +The log likelihood is + +$$ +\ell(\lambda) = n \ln \lambda - \lambda \sum_{i=1}^n x_i +$$ + +Setting the derivative to zero gives $n / \lambda - \sum_i x_i = 0$, so $$ -\hat{\lambda} = \frac{n}{\sum_{i=1}^n w_i} +\hat \lambda = \frac{n}{\sum_{i=1}^n x_i} = \frac{1}{\bar x} $$ -1. Compute $\hat{\lambda}$ for our initial sample. -2. Use $\hat{\lambda}$ to find the total revenue +Here is the numerical check. + +```{code-cell} ipython3 +rng = np.random.default_rng(1) +x = rng.exponential(scale=1/2, size=500) + +_, scale_fit = scipy.stats.expon.fit(x, floc=0) +1 / x.mean(), 1 / scale_fit +``` + +```{solution-end} +``` + +```{exercise-start} +:label: mle_ex_poisson +``` + +Repeat the consistency simulation of {ref}`mle_consistency` for the Poisson +distribution with λ = 2.7, using sample sizes 10, 100 and 1,000. + +Plot histograms of the estimates $\hat \lambda$ for each sample size. ```{exercise-end} ``` -```{solution-start} mle_ex1 +```{solution-start} mle_ex_poisson :class: dropdown ``` ```{code-cell} ipython3 -λ_hat = 1/np.mean(sample) -λ_hat -``` +rng = np.random.default_rng(2) +λ_true = 2.7 +num_reps = 2_000 -```{code-cell} ipython3 -dist_exp = expon(scale=1/λ_hat) -tr_expo = total_revenue(dist_exp) -tr_expo +fig, ax = plt.subplots() +for n in (10, 100, 1_000): + x = rng.poisson(λ_true, size=(num_reps, n)) + λ_hats = x.mean(axis=1) + ax.hist(λ_hats, bins=40, density=True, histtype='step', lw=2, + label=f'$n = {n}$') +ax.axvline(λ_true, color='k', ls=':', lw=1) +ax.set_xlabel(r'estimate of $\lambda$') +ax.set_ylabel('density') +ax.legend() +plt.show() ``` +The estimates concentrate on the true value as $n$ grows. + ```{solution-end} ``` ```{exercise-start} -:label: mle_ex2 +:label: mle_ex_uniform ``` -Plot the exponential distribution against the sample and check if it is a good fit or not. +Suppose that $x_1, \ldots, x_n$ are IID draws from the uniform distribution on +$[0, \theta]$, where θ > 0 is unknown. + +1. Show that the likelihood is $\theta^{-n}$ when $\theta \geq \max_i x_i$ and + zero otherwise. +2. Conclude that the maximum likelihood estimate is + $\hat \theta = \max_i x_i$. +3. Plot the likelihood for a sample of size 10 with θ = 1. + +Notice that the maximizer is not found by setting a derivative to zero. ```{exercise-end} ``` -```{solution-start} mle_ex2 +```{solution-start} mle_ex_uniform :class: dropdown ``` +The density of each observation is $1/\theta$ on $[0, \theta]$ and zero +elsewhere. + +If any observation exceeds θ, that observation has zero density, so the +likelihood is zero. + +Otherwise the likelihood is $(1/\theta)^n = \theta^{-n}$. + +Since $\theta^{-n}$ is decreasing in θ, the likelihood is maximized by the +smallest θ that is compatible with the data, which is $\max_i x_i$. + ```{code-cell} ipython3 -fig, ax = plt.subplots() -ax.set_xlim(-1, 20) +rng = np.random.default_rng(3) +x = rng.uniform(0, 1, size=10) -ax.hist(sample, density=True, bins=5000, histtype='stepfilled', alpha=0.5) -ax.plot(x, dist_exp.pdf(x), 'k-', lw=2, label='exponential pdf') -ax.legend() +θ_grid = np.linspace(0.5, 1.5, 500) +L = np.where(θ_grid >= x.max(), θ_grid**(-len(x)), 0.0) +fig, ax = plt.subplots() +ax.plot(θ_grid, L, lw=2) +ax.axvline(x.max(), color='k', ls='--', lw=1, label=r'$\max_i x_i$') +ax.set_xlabel(r'$\theta$') +ax.set_ylabel('likelihood') +ax.legend() plt.show() ``` -Clearly, this distribution is not a good fit for our data. +The likelihood jumps up at $\max_i x_i$ and declines afterwards, so the maximum +is at the jump. + +```{solution-end} +``` + +```{exercise-start} +:label: mle_ex_bias +``` + +The maximum likelihood estimate of $\sigma^2$ for the normal distribution is +$\hat \sigma^2 = \frac{1}{n} \sum_i (x_i - \bar x)^2$. + +Using simulation with samples of size $n = 5$ from the standard normal +distribution, estimate the mean of $\hat \sigma^2$ and compare it with the +true value $\sigma^2 = 1$. + +Repeat the exercise using the divisor $n - 1$ instead of $n$. + +What do you conclude? + +```{exercise-end} +``` + +```{solution-start} mle_ex_bias +:class: dropdown +``` + +```{code-cell} ipython3 +rng = np.random.default_rng(4) +n, num_reps = 5, 100_000 +x = rng.standard_normal((num_reps, n)) + +x.var(axis=1).mean(), x.var(axis=1, ddof=1).mean() +``` + +With divisor $n$, the mean of the estimates is close to $(n - 1)/n = 0.8$, +so the maximum likelihood estimate of $\sigma^2$ is biased downward. + +With divisor $n - 1$, the mean is close to one, so this estimate is unbiased. + +The bias of the maximum likelihood estimate vanishes as $n$ grows, which is +consistent with what we saw in {ref}`mle_consistency`. ```{solution-end} ``` From e367333d9bd4e3b00a26814c7e6fce4139bd502e Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Tue, 29 Sep 2026 10:58:35 +1000 Subject: [PATCH 4/5] mle: use a seed that gives typical first samples Co-Authored-By: Claude Opus 5.5 --- lectures/mle.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lectures/mle.md b/lectures/mle.md index c3eb4b65..276c8c21 100644 --- a/lectures/mle.md +++ b/lectures/mle.md @@ -122,7 +122,7 @@ maximizes $L(\theta)$. Let's generate some data, using a true default probability of 0.15. ```{code-cell} ipython3 -rng = np.random.default_rng(1234) +rng = np.random.default_rng(9) θ_true = 0.15 n = 20 y = (rng.random(n) < θ_true).astype(int) From b9715a8f44b3a26306502911a7b2462c82b60add Mon Sep 17 00:00:00 2001 From: John Stachurski Date: Wed, 30 Sep 2026 07:07:50 +1000 Subject: [PATCH 5/5] mle: revise first example, normal case, and method of moments comparison - Go slowly through the Bernoulli example: likelihood as a product, plot, then the log likelihood derivation; use a seed where the estimate differs from the truth - Move the underflow point to the numerical MLE section, using house prices - Identify the normal MLEs as sample mean and standard deviation - Drop the asymptotic normality discussion - Reframe the method of moments comparison: invalid for alpha <= 1, unstable for heavy tails Co-Authored-By: John Stachurski Co-Authored-By: Claude Opus 5.5 --- lectures/mle.md | 436 +++++++++++++++++++++++++++++------------------- 1 file changed, 265 insertions(+), 171 deletions(-) diff --git a/lectures/mle.md b/lectures/mle.md index 276c8c21..da1713d7 100644 --- a/lectures/mle.md +++ b/lectures/mle.md @@ -15,49 +15,34 @@ kernelspec: ## Overview -Maximum likelihood estimation is a method for choosing the parameters of a -statistical model. +Maximum likelihood estimation is one of the most important methods for choosing +the parameters of a statistical model. -The idea is simple: among all candidate parameter values, choose the one under -which the data we actually observed are most probable. +The basic idea is to choose the parameter values that make the data we actually +observed most probable. Put differently, we treat the probability that the model assigns to the data as a function of the parameters, and then maximize it. -This principle is one of the foundations of modern statistics. +This principle is one of the foundations of modern statistics, with applications +across all scientific domains, including machine learning. -It is used to fit distributions, as in this lecture, and also regression models, -macroeconomic models, and models in physics, biology and machine learning. - -(For example, training a classifier by minimizing cross-entropy loss is maximum -likelihood estimation.) - -One reason for its wide use is generality: once we have a model that assigns -probabilities to data, maximum likelihood tells us how to fit it, with no need -to invent a new method for each model. - -Another is that, under fairly general conditions, maximum likelihood estimates -converge to the true parameter values as the sample grows, and in large samples -no other method gives more precise estimates. - -The idea has a long history. - -Daniel Bernoulli and Carl Friedrich Gauss used versions of it in the late 18th +Daniel Bernoulli and Carl Friedrich Gauss used versions of maximum in the late 18th and early 19th centuries, and Gauss justified the method of least squares by showing that, when errors are normally distributed, it picks out the most probable parameter values. -The method was developed systematically, and given its name, by Ronald Fisher -between 1912 and 1922 {cite}`fisher1922mathematical`. +The method was developed systematically by Ronald Fisher over the decade from 1912 {cite}`fisher1922mathematical`. -{cite:t}`aldrich1997fisher` and {cite:t}`stigler2007epic` tell the story. +(See {cite:t}`aldrich1997fisher` and {cite:t}`stigler2007epic` for more historical details.) -In this lecture we build up the idea through a sequence of examples, check by -simulation that maximum likelihood estimates converge to the truth, and compare -maximum likelihood with the method of moments from {doc}`fitting_distributions`. +This lecture provides a quick introduction to maximum likelihood, focusing on examples. -We also show how to compute maximum likelihood estimates numerically when no -formula is available, and apply this to heavy-tailed stock returns. +We will also + +* check by simulation that maximum likelihood estimates converge to the truth as more information becomes available, +* compare maximum likelihood with the method of moments from {doc}`fitting_distributions`, and +* how to compute maximum likelihood estimates numerically when no formula is available. In addition to what's in Anaconda, this lecture will need the following libraries: @@ -83,8 +68,11 @@ import yfinance as yf We start with an example that we will meet again in {doc}`bayes_intro`. -A bank makes a series of small loans, each of which either defaults or is -repaid. +Since this is our first example, we will go slowly. + +### The model + +A bank makes a series of small loans, each of which either defaults or is repaid. We encode the outcome of loan $i$ as $Y_i = 1$ if it defaults and $Y_i = 0$ if it is repaid. @@ -100,37 +88,60 @@ $$ The default probability θ is unknown, and our job is to estimate it from the observed outcomes $y_1, \ldots, y_n$. -Both cases can be written as a single expression, since -$\mathbb{P}\{Y_i = y\} = \theta^y (1 - \theta)^{1 - y}$ for $y \in \{0, 1\}$. +Note that both cases above can be written via the single expression + +$$ + \mathbb{P}\{Y_i = y\} = \theta^y (1 - \theta)^{1 - y} + \quad \text{for} \quad y \in \{0, 1\}. +$$ + +(To check this, set $y = 1$ and then $y = 0$.) + +### The likelihood -By independence, the probability of observing the whole sample is therefore +By independence, the probability distribution for the random vector +$(Y_1,\ldots, Y_n)$ is the product of the marginals. + +Thus, the probability of observing the whole sample $y_1, \ldots, y_n$ is $$ -L(\theta) -= \prod_{i=1}^n \theta^{y_i} (1 - \theta)^{1 - y_i} -= \theta^k (1 - \theta)^{n - k} +L(\theta) := \prod_{i=1}^n \theta^{y_i} (1 - \theta)^{1 - y_i} $$ (eq:bern_lik) -where $k = \sum_i y_i$ is the number of defaults. +The function $L$ is called the **likelihood function**. -Viewed as a function of θ, with the data held fixed, $L$ is called the -**likelihood function**. +It takes the observed data as fixed and writes the probability of observing the +data as a function of the parameter. The **maximum likelihood estimate** $\hat \theta$ is the value of θ that maximizes $L(\theta)$. +In other words, it is the value of θ under which the data we actually observed +are most probable. + +### Plotting the likelihood + Let's generate some data, using a true default probability of 0.15. ```{code-cell} ipython3 -rng = np.random.default_rng(9) -θ_true = 0.15 -n = 20 -y = (rng.random(n) < θ_true).astype(int) -k = y.sum() +seed = 21 # Arbitrarily chosen seed for the random number generator +rng = np.random.default_rng(seed) +θ_true = 0.15 # Set the "unknown" true parameter to 0.15 +n = 20 # Fix a relatively small sample size +y = (rng.random(n) < θ_true).astype(int) # Generate n independent draws y ``` -Now we plot the likelihood function. +Here is a function that computes the likelihood {eq}`eq:bern_lik`, directly +from its definition as a product. + +```{code-cell} ipython3 +def likelihood(θ, y): + return np.prod(θ**y * (1 - θ)**(1 - y)) +``` + +Now we plot the likelihood function, by evaluating it on a grid of values +for θ. ```{code-cell} ipython3 --- @@ -140,23 +151,113 @@ mystnb: name: fig:mle-bernoulli-lik --- θ_grid = np.linspace(0.001, 0.999, 500) - -def bernoulli_likelihood(θ, k, n): - return θ**k * (1 - θ)**(n - k) +L_values = [likelihood(θ, y) for θ in θ_grid] fig, ax = plt.subplots() -ax.plot(θ_grid, bernoulli_likelihood(θ_grid, k, n), lw=2) -ax.axvline(k / n, color='k', ls='--', lw=1, label='sample mean') +ax.plot(θ_grid, L_values, lw=2) ax.set_xlabel(r'$\theta$') ax.set_ylabel('likelihood') -ax.legend() plt.show() ``` -The likelihood peaks at the sample mean $k/n$, the fraction of loans that -defaulted. +The likelihood is maximized at a value of θ somewhere around 0.1. + +We can locate the maximizer more precisely by finding the grid point with the +largest likelihood. + +```{code-cell} ipython3 +θ_grid[np.argmax(L_values)] +``` + +### Maximizing the likelihood with calculus -This is the value of θ under which the data we observed are most probable. +Let's now find the maximizer exactly, using calculus. + +We could differentiate $L$ directly, but differentiating a product of $n$ +terms is awkward. + +It is easier to work with the logarithm of the likelihood, called the +**log likelihood**: + +$$ +\ell(\theta) := \ln L(\theta) +$$ + +Since the logarithm is strictly increasing, $L$ and $\ell$ have the same +maximizer. + +So we can maximize $\ell$ instead of $L$. + +The advantage is that the logarithm turns products into sums. + +Using $\ln (ab) = \ln a + \ln b$ and $\ln (a^b) = b \ln a$, we get + +$$ +\ell(\theta) += \sum_{i=1}^n \ln \left[ \theta^{y_i} (1 - \theta)^{1 - y_i} \right] += \sum_{i=1}^n \left[ y_i \ln \theta + (1 - y_i) \ln (1 - \theta) \right] +$$ + +Let $k := \sum_{i=1}^n y_i$ be the number of defaults. + +Then $\sum_{i=1}^n (1 - y_i) = n - k$ is the number of loans repaid, and +collecting terms gives + +$$ +\ell(\theta) = k \ln \theta + (n - k) \ln (1 - \theta) +$$ + +Its derivative is + +$$ +\ell'(\theta) = \frac{k}{\theta} - \frac{n - k}{1 - \theta} +$$ + +Assuming that $0 < k < n$, we set the derivative to zero and multiply through by +$\theta (1 - \theta)$, which gives + +$$ +k (1 - \theta) - (n - k) \theta = 0 +\quad \iff \quad +k - n \theta = 0 +\quad \iff \quad +\theta = \frac{k}{n} +$$ + +This point is a maximum, since the second derivative + +$$ +\ell''(\theta) = - \frac{k}{\theta^2} - \frac{n - k}{(1 - \theta)^2} +$$ + +is negative everywhere on $(0, 1)$. + +Hence the maximum likelihood estimate is + +$$ +\hat \theta = \frac{k}{n} = \frac{1}{n} \sum_{i=1}^n y_i +$$ + +This is the sample mean of the data, or, in other words, the fraction of loans +that defaulted. + +Let's check that it agrees with the figure. + +```{code-cell} ipython3 +k = y.sum() +k / n +``` + +This matches the grid search above, up to the spacing of the grid. + +Notice that the estimate is not equal to the true value 0.15. + +With only 20 loans, the fraction that defaulted is a noisy guide to θ. + +We will see {ref}`below ` that the estimate improves as we +observe more loans. + +### Likelihood is not probability It is important to be clear about what $L$ is. @@ -168,11 +269,25 @@ and need not integrate to one. That is why we give it a different name. +(mle_more_data)= +### More data + Let's see what happens to the likelihood as we observe more loans. +Exponentiating our expression for $\ell$ shows that +$L(\theta) = \theta^k (1 - \theta)^{n - k}$, which is quicker to compute +than the product when $n$ is large. + +```{code-cell} ipython3 +def bernoulli_likelihood(θ, k, n): + return θ**k * (1 - θ)**(n - k) +``` + Since the likelihood gets very small when $n$ is large, we divide each curve by its maximum value, so that all the curves peak at one. +(We return to the problem of very small likelihoods {ref}`below `.) + ```{code-cell} ipython3 --- mystnb: @@ -202,61 +317,6 @@ This is a first sign of the consistency of maximum likelihood, which we study {ref}`below `. -## The log likelihood - -In practice we almost always work with the logarithm of the likelihood, -called the **log likelihood**, rather than the likelihood itself. - -There are two reasons. - -The first is that the logarithm turns products into sums, which are easier to -differentiate. - -The second is numerical. - -The likelihood is a product of many numbers less than one, and for large -samples it becomes too small to represent on a computer. - -```{code-cell} ipython3 -n = 5_000 -y = (rng.random(n) < θ_true).astype(int) -k = y.sum() -bernoulli_likelihood(θ_true, k, n) -``` - -The true value is positive, but it is smaller than the smallest positive -floating point number, so the computer rounds it to zero. - -The log likelihood has no such problem. - -```{code-cell} ipython3 -def bernoulli_log_likelihood(θ, k, n): - return k * np.log(θ) + (n - k) * np.log(1 - θ) - -bernoulli_log_likelihood(θ_true, k, n) -``` - -Since the logarithm is strictly increasing, maximizing the log likelihood gives -the same answer as maximizing the likelihood. - -For the loan data, the log likelihood is - -$$ -\ell(\theta) = \ln L(\theta) = k \ln \theta + (n - k) \ln (1 - \theta) -$$ - -Setting the derivative to zero gives - -$$ -\frac{k}{\theta} - \frac{n - k}{1 - \theta} = 0 -\quad \implies \quad -\hat \theta = \frac{k}{n} -$$ - -So the maximum likelihood estimate is the sample mean, confirming what we saw -in the figure. - - ## More examples In this section we derive maximum likelihood estimates for three more @@ -356,6 +416,12 @@ $$ \hat \sigma = \left( \frac{1}{n} \sum_{i=1}^n (x_i - \hat \mu)^2 \right)^{1/2} $$ +In other words, $\hat \mu$ is the sample mean of the data, and $\hat \sigma$ +is the sample standard deviation. + +(Here the sample standard deviation is computed with divisor $n$, a point we +return to below.) + Let's look at the log likelihood for the heights of US adult women, which we fitted with a normal distribution in {doc}`fitting_distributions`. @@ -369,6 +435,10 @@ female = heights[heights['sex'] == 'female']['height_cm'].to_numpy() With two parameters, the log likelihood is a surface, which we display with a contour plot. +We compute $\hat \mu$ and $\hat \sigma$ using NumPy's `mean` and `std` +methods, which implement the formulas above, and mark the point +$(\hat \mu, \hat \sigma)$ on the plot. + ```{code-cell} ipython3 --- mystnb: @@ -398,7 +468,9 @@ ax.legend() plt.show() ``` -The maximum sits at the sample mean and standard deviation. +As expected, the log likelihood is maximized at $(\hat \mu, \hat \sigma)$. + +Here are the values. ```{code-cell} ipython3 μ_hat, σ_hat @@ -521,59 +593,34 @@ plt.show() As $n$ grows, the estimates concentrate on the true value. -How fast do they concentrate? - -The {doc}`central limit theorem ` suggests that the spread should -shrink like $1/\sqrt{n}$. +This property holds for maximum likelihood estimates in general, under +conditions that are satisfied in most applications. -For this estimator it is known that, for large $n$, $\hat \alpha$ is -approximately normal with mean α and standard deviation $\alpha / \sqrt{n}$. +(Proofs can be found in any graduate text on statistics.) -Let's check by standardizing the estimates and comparing them with a standard -normal density. -```{code-cell} ipython3 ---- -mystnb: - figure: - caption: Standardized estimates and the standard normal - name: fig:mle-asymptotic-normal ---- -n = 1_000 -x = pareto_draws(α_true, (num_reps, n), rng) -α_hats = n / np.sum(np.log(x), axis=1) -z = (α_hats - α_true) / (α_true / np.sqrt(n)) - -z_grid = np.linspace(-4, 4, 200) -fig, ax = plt.subplots() -ax.hist(z, bins=50, density=True, alpha=0.4, label='standardized estimates') -ax.plot(z_grid, scipy.stats.norm.pdf(z_grid), lw=2, label='standard normal') -ax.set_xlabel('standardized estimate') -ax.set_ylabel('density') -ax.legend() -plt.show() -``` - -The fit is close. +(mle_vs_mom)= +## Maximum likelihood and the method of moments -These two properties, consistency and approximate normality in large samples, -hold for maximum likelihood estimates in general, under regularity conditions -that are satisfied in most applications. +In {doc}`fitting_distributions` we fitted distributions using the **method of +moments**, which chooses parameters so that the moments of the fitted +distribution (its mean, variance and so on) match the corresponding moments of +the data. -Moreover, the spread of the approximating normal distribution is as small as -that of any reasonable estimator, which is the sense in which maximum -likelihood makes the best use of the data. +Let's compare this method with maximum likelihood for the Pareto distribution. -Proofs can be found in any graduate text on statistics. +The Pareto distribution has one unknown parameter α, so the method of moments +uses one moment, the mean. +This immediately raises a problem. -(mle_vs_mom)= -## Maximum likelihood and the method of moments +As we saw in {doc}`heavy_tails`, the mean of the Pareto distribution is finite +only when α > 1. -Let's now compare maximum likelihood with the method of moments for the Pareto -distribution. +When α ≤ 1 the mean is infinite, so there is no population mean for the sample +mean to match, and the method of moments cannot be used. -When α > 1, the Pareto distribution has mean $\alpha \bar x / (\alpha - 1)$. +When α > 1, the mean is $\alpha \bar x / (\alpha - 1)$. Setting this equal to the sample mean $m_n = \frac{1}{n} \sum_i x_i$ and solving for α gives the method of moments estimate @@ -588,26 +635,41 @@ def pareto_mom(x, x_bar=1.0): return m / (m - x_bar) ``` -Since every observation is at least $\bar x$, we have $m_n > \bar x$ and hence -$\tilde \alpha > 1$ for *every* data set. +In practice, of course, we do not know α in advance, so we do not know whether +the method is valid. + +The danger is that the formula can be computed for any sample, since every +sample has a finite mean. + +Moreover, since every observation is at least $\bar x$, we have +$m_n > \bar x$, so the formula always returns a value greater than one. + +So if the true α is at or below one, the method of moments quietly gives an +answer that is on the wrong side of one, suggesting that the mean is finite +when in fact it is not. -So if the true α is less than or equal to one, the method of moments can never -recover it, however much data we have. +Even when α > 1, so that the method is valid, heavy tails cause trouble. -The underlying problem is that when α ≤ 1 the mean is infinite, so there is no -population moment for the sample mean to match. +When α ≤ 2 the variance is infinite, and, as we saw in {doc}`heavy_tails`, the +sample mean then settles down only slowly, because it is dominated by a few very +large draws. -Even when α is a little above one, so that the mean exists, the variance is -infinite whenever α ≤ 2. +The method of moments estimate inherits this instability. -Then, as we saw in {doc}`heavy_tails`, the sample mean settles down only -slowly, and the method of moments estimate inherits this instability. +Maximum likelihood has neither problem. -Maximum likelihood suffers from neither problem, since it works with the -density, which is well defined for every α > 0. +It works with the density, which is well defined for every α > 0. -Let's compare the two estimators by simulation, for α = 0.8 and α = 1.5, with -samples of size 1,000. +Also, the estimate {eq}`eq:pareto_mle` depends on the data only through the +average of $\ln (x_i / \bar x)$, and taking logs tames the very large draws. + +(In fact, if $X$ is Pareto, then $\ln (X / \bar x)$ is exponentially +distributed, which is a light-tailed distribution.) + +Let's compare the two estimators by simulation, with samples of size 1,000. + +We take α = 0.8, where the method of moments is not valid, and α = 1.5, where +it is valid but the variance is infinite. ```{code-cell} ipython3 --- @@ -636,8 +698,9 @@ plt.tight_layout() plt.show() ``` -When α = 0.8 (left), the method of moments estimates all lie above one, far -from the truth, while the maximum likelihood estimates are centered on it. +When α = 0.8 (left), the method of moments estimates all lie above one, as +predicted, far from the truth, while the maximum likelihood estimates are +centered on it. When α = 1.5 (right), both estimators are concentrated near the truth, but the method of moments estimates are much more dispersed, with a long left tail @@ -654,6 +717,7 @@ this is one reason why maximum likelihood is the standard method for fitting heavy-tailed distributions. +(mle_numerical)= ## Numerical maximum likelihood In all the examples so far, we could find the maximizer with pencil and paper. @@ -680,6 +744,36 @@ so we will compute it numerically. Most numerical optimization routines minimize rather than maximize, so we minimize the *negative* of the log likelihood. +Why the log likelihood, rather than the likelihood itself? + +When we maximized by hand, logs made the calculus easier. + +For numerical work there is a second, more important reason. + +The likelihood is a product of $n$ densities, one for each observation, and +with thousands of observations this product is too small to represent on a +computer. + +Here is the likelihood of the house price data at a typical parameter value, +$a = 6$ and $s = 30$. + +```{code-cell} ipython3 +u = scipy.stats.gamma(6, scale=30) +np.prod(u.pdf(price)) +``` + +The true value is positive, but it is smaller than the smallest positive +floating point number, so the computer rounds it to zero. + +The log likelihood, which is a sum rather than a product, has no such problem. + +```{code-cell} ipython3 +np.sum(u.logpdf(price)) +``` + +For this reason, numerical routines for maximum likelihood always work with the +log likelihood. + To keep both parameters positive, we optimize over their logarithms, which can take any real value.