{ "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernelspec": { "name": "python3", "display_name": "Python 3" }, "language_info": { "name": "python" }, "colab": { "provenance": [] } }, "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Day 11: Maximum Likelihood Estimation\n", "\n", "*Beginner to Ninja: Math for Causal Inference*\n", "\n", "Choosing the parameter that makes the observed data most probable, and measuring its precision from the curvature.\n", "\n", "**Learning outcomes**\n", "\n", "- Write the likelihood and log-likelihood for Bernoulli, Poisson, normal and exponential data.\n", "- Derive MLEs from the first-order condition and check the second-order condition.\n", "- Use Fisher information to get standard errors, and the invariance property for transformations.\n", "- Run a likelihood-ratio test and connect OLS and logit to maximum likelihood." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory\n", "\n", "### Likelihood\n", "\n", "For iid data $x_1, \\dots, x_n$ from a density or pmf $f(x; \\theta)$, the **likelihood** is the joint probability of the data, read as a function of $\\theta$:\n", "\n", "$$\n", "L(\\theta) = \\prod_{i=1}^n f(x_i; \\theta), \\qquad \\ell(\\theta) = \\ln L(\\theta) = \\sum_{i=1}^n \\ln f(x_i; \\theta).\n", "$$\n", "\n", "Probability fixes $\\theta$ and varies the data; likelihood fixes the data and varies $\\theta$. Logs turn products into sums and do not move the maximiser.\n", "\n", "### The maximum likelihood estimator\n", "\n", "$\\hat\\theta_{MLE}$ maximises $\\ell(\\theta)$. Usually you solve the **score equation** $\\ell'(\\theta) = 0$ and check $\\ell''(\\hat\\theta) < 0$.\n", "\n", "### Precision from curvature\n", "\n", "The **Fisher information** per observation is $I(\\theta) = -E\\big[\\partial^2 \\ln f(X; \\theta)/\\partial\\theta^2\\big]$. Under regularity conditions,\n", "\n", "$$\n", "\\hat\\theta_{MLE} \\approx N\\Big(\\theta,\\; \\frac{1}{n\\,I(\\theta)}\\Big).\n", "$$\n", "\n", "A sharply peaked log-likelihood (large curvature) means a precise estimate. In practice, $\\operatorname{SE} \\approx 1/\\sqrt{-\\ell''(\\hat\\theta)}$.\n", "\n", "### Two more tools\n", "\n", "- **Invariance:** the MLE of $g(\\theta)$ is $g(\\hat\\theta)$.\n", "- **Likelihood-ratio test:** $LR = 2\\big[\\ell(\\hat\\theta) - \\ell(\\theta_0)\\big] \\approx \\chi^2_1$ under $H_0: \\theta = \\theta_0$; reject at 5% if $LR > 3.84$.\n", "\n", "> **Causal link.** Logit and probit propensity-score models, duration models and many structural models are estimated by maximum likelihood. On Day 13 you will see that the posterior mode with a flat prior is the MLE.\n", "\n", "## Part 1 · Warm-up" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1 · Probability or likelihood?\n", "\n", "Which is a probability and which a likelihood: (a) \"if $p = 0.5$, the chance of 3 heads in 4 tosses is 0.25\"; (b) \"given 3 heads in 4 tosses, $p = 0.75$ fits better than $p = 0.5$\"?\n", "\n", "
Show solution\n", "\n", "**Step 1.** (a) is a probability: the parameter is fixed and we ask about data.\n", "\n", "**Step 2.** (b) compares likelihoods: the data are fixed and we compare parameter values.\n", "\n", "**Step 3.** Likelihoods need not sum to 1 over $\\theta$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2 · A Bernoulli likelihood\n", "\n", "Data: 1, 1, 0, 1. Write $L(p)$ and evaluate it at $p = 0.5$ and $p = 0.75$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $L(p) = p^3(1 - p)$.\n", "\n", "**Step 2.** $L(0.5) = 0.0625$ and $L(0.75) = 0.421875 \\times 0.25 \\approx 0.105$.\n", "\n", "**Step 3.** The data are more probable under $p = 0.75$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3 · Log-likelihood\n", "\n", "Write $\\ell(p)$ for the same data.\n", "\n", "
Show solution\n", "\n", "$\\ell(p) = 3\\ln p + \\ln(1 - p)$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 4 · Maximise\n", "\n", "Find the MLE for these data.\n", "\n", "
Show solution\n", "\n", "$\\ell'(p) = 3/p - 1/(1 - p) = 0 \\Rightarrow 3(1 - p) = p \\Rightarrow \\hat p = 3/4$. And $\\ell''(p) = -3/p^2 - 1/(1 - p)^2 < 0$, so it is a maximum.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 5 · Exponential durations\n", "\n", "Five unemployment spells last 2, 3, 4, 5 and 6 months. With $f(t) = \\lambda e^{-\\lambda t}$, find $\\hat\\lambda$.\n", "\n", "
Show solution\n", "\n", "$\\ell(\\lambda) = n\\ln\\lambda - \\lambda\\sum t_i$, so $n/\\lambda - \\sum t_i = 0$ and $\\hat\\lambda = 1/\\bar t = 1/4 = 0.25$ per month (mean spell 4 months).\n", "\n", "
\n", "\n", "## Part 2 · Proofs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 6 · Bernoulli MLE in general\n", "\n", "With $k$ successes in $n$ trials, show $\\hat p = k/n$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\ell(p) = k\\ln p + (n - k)\\ln(1 - p)$.\n", "\n", "**Step 2.** Setting $k/p - (n - k)/(1 - p) = 0$ gives $k(1 - p) = (n - k)p$, so $\\hat p = k/n$.\n", "\n", "**Step 3.** The second derivative is negative everywhere, so this is the global maximum.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 7 · Poisson MLE\n", "\n", "Show that the MLE of $\\lambda$ for iid Poisson data is $\\bar x$.\n", "\n", "
Show solution\n", "\n", "$\\ell(\\lambda) = \\sum_i\\big(x_i\\ln\\lambda - \\lambda - \\ln x_i!\\big)$, so $\\ell'(\\lambda) = \\sum x_i/\\lambda - n = 0 \\Rightarrow \\hat\\lambda = \\bar x$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 8 · Normal mean\n", "\n", "For $X_i \\sim N(\\mu, \\sigma^2)$ iid, show $\\hat\\mu = \\bar x$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\ell(\\mu, \\sigma^2) = -\\frac n2\\ln(2\\pi\\sigma^2) - \\frac{1}{2\\sigma^2}\\sum(x_i - \\mu)^2$.\n", "\n", "**Step 2.** Maximising over $\\mu$ means minimising $\\sum(x_i - \\mu)^2$, a least-squares problem: $\\partial/\\partial\\mu = \\frac{1}{\\sigma^2}\\sum(x_i - \\mu) = 0 \\Rightarrow \\hat\\mu = \\bar x$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9 · Normal variance\n", "\n", "Show $\\hat\\sigma^2_{MLE} = \\frac1n\\sum(x_i - \\bar x)^2$. Is it unbiased?\n", "\n", "
Show solution\n", "\n", "**Step 1.** $\\partial\\ell/\\partial\\sigma^2 = -\\frac{n}{2\\sigma^2} + \\frac{1}{2\\sigma^4}\\sum(x_i - \\hat\\mu)^2 = 0 \\Rightarrow \\hat\\sigma^2 = \\frac1n\\sum(x_i - \\bar x)^2$.\n", "\n", "**Step 2.** This is the plug-in variance from Day 8: biased (by $-\\sigma^2/n$) but consistent.\n", "\n", "**Step 3.** MLEs are often biased in small samples.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 10 · Fisher information for a proportion\n", "\n", "Show $I(p) = 1/(p(1 - p))$ for one Bernoulli observation, so $\\operatorname{Var}(\\hat p) \\approx p(1 - p)/n$.\n", "\n", "
Show solution\n", "\n", "$\\ln f(x; p) = x\\ln p + (1 - x)\\ln(1 - p)$, so $\\partial^2/\\partial p^2 = -x/p^2 - (1 - x)/(1 - p)^2$. Taking expectations with $E[x] = p$:\n", "\n", "$$\n", "I(p) = \\frac{p}{p^2} + \\frac{1 - p}{(1 - p)^2} = \\frac1p + \\frac{1}{1 - p} = \\frac{1}{p(1 - p)}.\n", "$$\n", "\n", "So $1/(nI(p)) = p(1 - p)/n$, the familiar variance of a proportion.\n", "\n", "
\n", "\n", "## Part 3 · Applied" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 11 · A conversion rate\n", "\n", "18 of 30 households adopt a new technology. Find $\\hat p$ and its standard error.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Bernoulli MLE (Exercise 6): $\\hat p = 18/30 = 0.6$.\n", "\n", "**Step 2.** Fisher information gives $\\operatorname{Var}(\\hat p) \\approx \\hat p(1 - \\hat p)/n = 0.6 \\times 0.4/30 = 0.008$.\n", "\n", "**Step 3.** $\\operatorname{SE} = \\sqrt{0.008} \\approx 0.0894$.\n", "\n", "**Answer:** 60% adoption, SE ≈ 0.089. Slide $N$ and $S$ in the Lab's likelihood explorer to see the curve sharpen as $N$ grows.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 12 · Invariance\n", "\n", "Using Exercise 11, find the MLE of the odds $p/(1 - p)$ and of the log-odds.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Invariance: the MLE of a function of $p$ is that function of $\\hat p$.\n", "\n", "**Step 2.** Odds: $\\hat p/(1 - \\hat p) = 0.6/0.4 = 1.5$ (three adopters for every two non-adopters).\n", "\n", "**Step 3.** Log-odds: $\\ln 1.5 \\approx 0.405$, the scale on which a logit model works.\n", "\n", "**Answer:** odds 1.5 and log-odds 0.405. Their standard errors come from the delta method (Day 12).\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 13 · A likelihood-ratio test\n", "\n", "Test $H_0: p = 0.5$ for the data in Exercise 11.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Log-likelihood at the MLE:\n", "\n", "$$\n", "\\ell(0.6) = 18\\ln 0.6 + 12\\ln 0.4 = 18(-0.5108) + 12(-0.9163) \\approx -20.190.\n", "$$\n", "\n", "**Step 2.** Log-likelihood under $H_0$:\n", "\n", "$$\n", "\\ell(0.5) = 30\\ln 0.5 = 30(-0.6931) \\approx -20.794.\n", "$$\n", "\n", "**Step 3.** $LR = 2[\\ell(0.6) - \\ell(0.5)] = 2(0.604) \\approx 1.21$.\n", "\n", "**Step 4.** Compare with the $\\chi^2_1$ 5% critical value 3.84: $1.21 < 3.84$ ($p \\approx 0.27$).\n", "\n", "**Answer:** do not reject $p = 0.5$. The Wald test agrees: $z = 0.1/0.0894 \\approx 1.12$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 14 · OLS is MLE in a wage equation\n", "\n", "A labour economist models log wages as $y_i = \\beta_0 + \\beta_1\\,\\text{school}_i + \\varepsilon_i$ with $\\varepsilon_i \\sim N(0, \\sigma^2)$ iid. Show that the maximum-likelihood estimates of $(\\beta_0, \\beta_1)$ are the OLS estimates.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Each $y_i$ is $N(\\beta_0 + \\beta_1\\text{school}_i,\\; \\sigma^2)$, so the log-likelihood is\n", "\n", "$$\n", "\\ell(\\beta_0, \\beta_1, \\sigma^2) = -\\frac n2\\ln(2\\pi\\sigma^2) - \\frac{1}{2\\sigma^2}\\sum_i(y_i - \\beta_0 - \\beta_1\\text{school}_i)^2.\n", "$$\n", "\n", "**Step 2.** The $\\beta$s appear only in the sum of squares, with a negative sign.\n", "\n", "**Step 3.** For any $\\sigma^2 > 0$, maximising $\\ell$ over the $\\beta$s is therefore the same as *minimising* the sum of squared residuals: the definition of OLS.\n", "\n", "**Answer:** under normal errors, MLE and OLS coincide. (Without normality OLS is still a sensible estimator, Day 29.)\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 15 · The propensity-score score equation\n", "\n", "An evaluator models the probability that worker $i$ joins a training programme with a logit, $P(T_i = 1 \\mid x_i) = \\Lambda(x_i'\\beta)$, $\\Lambda(u) = 1/(1 + e^{-u})$, where $x_i$ holds age, education and past earnings. Show that the score equation is $\\sum_i (T_i - \\Lambda(x_i'\\beta))\\,x_i = 0$, and interpret the equation for the intercept.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Log-likelihood: $\\ell(\\beta) = \\sum_i\\big[T_i\\ln\\Lambda_i + (1 - T_i)\\ln(1 - \\Lambda_i)\\big]$ with $\\Lambda_i = \\Lambda(x_i'\\beta)$.\n", "\n", "**Step 2.** Useful fact: $\\Lambda'(u) = \\Lambda(u)(1 - \\Lambda(u))$, so $\\partial\\Lambda_i/\\partial\\beta = \\Lambda_i(1 - \\Lambda_i)x_i$.\n", "\n", "**Step 3.** Chain rule:\n", "\n", "$$\n", "\\frac{\\partial\\ell}{\\partial\\beta} = \\sum_i\\Big[\\frac{T_i}{\\Lambda_i} - \\frac{1 - T_i}{1 - \\Lambda_i}\\Big]\\Lambda_i(1 - \\Lambda_i)x_i = \\sum_i\\big[T_i(1 - \\Lambda_i) - (1 - T_i)\\Lambda_i\\big]x_i = \\sum_i(T_i - \\Lambda_i)x_i.\n", "$$\n", "\n", "**Step 4.** Set it to zero. The intercept's entry of $x_i$ is 1, giving $\\sum_iT_i = \\sum_i\\Lambda_i$.\n", "\n", "**Answer:** the fitted propensities average exactly to the share of workers trained. This logit is the standard propensity-score model of Phase 3.\n", "\n", "
\n", "\n", "## Verify in Python" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "code", "metadata": {}, "execution_count": null, "outputs": [], "source": [ "import numpy as np\n", "from scipy import optimize, stats\n", "\n", "k, n = 18, 30\n", "loglik = lambda p: k * np.log(p) + (n - k) * np.log(1 - p)\n", "res = optimize.minimize_scalar(lambda p: -loglik(p), bounds=(1e-6, 1 - 1e-6), method=\"bounded\")\n", "p_hat = res.x\n", "h = 1e-5\n", "curv = (loglik(p_hat + h) - 2 * loglik(p_hat) + loglik(p_hat - h)) / h**2\n", "print(f\"p_hat = {p_hat:.4f}, SE from curvature = {1/np.sqrt(-curv):.4f}, \"\n", " f\"formula = {np.sqrt(p_hat*(1-p_hat)/n):.4f}\")\n", "\n", "LR = 2 * (loglik(0.6) - loglik(0.5))\n", "print(f\"LR = {LR:.3f}, p-value = {stats.chi2.sf(LR, 1):.3f}\")\n", "\n", "# Logit by maximum likelihood on simulated data\n", "rng = np.random.default_rng(11)\n", "x = rng.normal(size=20_000); X = np.column_stack([np.ones_like(x), x])\n", "T = rng.random(20_000) < 1 / (1 + np.exp(-(-0.5 + 1.0 * x)))\n", "nll = lambda b: -np.sum(T * (X @ b) - np.log1p(np.exp(X @ b)))\n", "b_hat = optimize.minimize(nll, np.zeros(2)).x\n", "print(\"logit estimates:\", b_hat.round(3), \"(truth -0.5, 1.0)\")\n", "print(\"score at optimum:\", (X.T @ (T - 1 / (1 + np.exp(-X @ b_hat)))).round(3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Key takeaways\n", "\n", "- The MLE picks the parameter under which the observed data are most probable.\n", "- Curvature of the log-likelihood (Fisher information) gives the standard error.\n", "- OLS with normal errors and the logit propensity model are both maximum-likelihood estimators." ] } ] }