{ "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 13: Bayesian Inference\n", "\n", "*Beginner to Ninja: Math for Causal Inference*\n", "\n", "Priors, likelihoods and posteriors; conjugate updating, shrinkage and credible intervals.\n", "\n", "**Learning outcomes**\n", "\n", "- Apply Bayes' theorem to parameters and compute Beta–binomial and normal–normal posteriors.\n", "- Show that the posterior mean is a precision-weighted average of the prior mean and the data.\n", "- Explain empirical-Bayes shrinkage and its link to ridge regression.\n", "- Contrast credible intervals with confidence intervals and compute a Bayesian A/B comparison." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Theory\n", "\n", "### Parameters as random variables\n", "\n", "Frequentists treat $\\theta$ as fixed and the data as random. Bayesians also put a distribution on $\\theta$ that describes what is known about it:\n", "\n", "$$\n", "\\underbrace{p(\\theta \\mid x)}_{\\text{posterior}} = \\frac{\\overbrace{p(x \\mid \\theta)}^{\\text{likelihood}}\\;\\overbrace{p(\\theta)}^{\\text{prior}}}{p(x)} \\;\\propto\\; p(x \\mid \\theta)\\,p(\\theta).\n", "$$\n", "\n", "The denominator $p(x) = \\int p(x \\mid \\theta)p(\\theta)\\,d\\theta$ only normalises, so we often work with \"∝\".\n", "\n", "### Two conjugate models\n", "\n", "**Beta–binomial.** Prior $p \\sim \\text{Beta}(\\alpha, \\beta)$, data $k$ successes and $f$ failures:\n", "\n", "$$\n", "p \\mid \\text{data} \\sim \\text{Beta}(\\alpha + k,\\; \\beta + f), \\qquad\n", "E[p \\mid \\text{data}] = \\frac{\\alpha + k}{\\alpha + \\beta + k + f}, \\qquad\n", "\\text{mode} = \\frac{\\alpha + k - 1}{\\alpha + \\beta + k + f - 2}.\n", "$$\n", "\n", "The prior acts like $\\alpha$ earlier successes and $\\beta$ earlier failures.\n", "\n", "**Normal–normal.** Prior $\\theta \\sim N(\\mu_0, \\tau_0^2)$ and an estimate $\\hat\\theta \\sim N(\\theta, s^2)$. With precisions (inverse variances) $\\pi_0 = 1/\\tau_0^2$ and $\\pi_d = 1/s^2$:\n", "\n", "$$\n", "\\theta \\mid \\hat\\theta \\sim N\\Big(\\frac{\\pi_0\\mu_0 + \\pi_d\\hat\\theta}{\\pi_0 + \\pi_d},\\; \\frac{1}{\\pi_0 + \\pi_d}\\Big).\n", "$$\n", "\n", "Precisions add, and the posterior mean **shrinks** the estimate toward the prior mean in proportion to how noisy the estimate is.\n", "\n", "### Summaries\n", "\n", "- **Posterior mean** (minimises expected squared loss), **MAP** (posterior mode).\n", "- A 95% **credible interval** contains $\\theta$ with posterior probability 0.95: the statement people often wrongly make about confidence intervals.\n", "- As $n$ grows, the likelihood dominates any reasonable prior and Bayesian and frequentist answers converge.\n", "\n", "> **Causal link.** Shrinkage is everywhere in applied causal work: teacher and hospital value-added, many-site experiments, and heterogeneous treatment effects all borrow strength across units this way.\n", "\n", "## Part 1 · Warm-up" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 1 · The three ingredients\n", "\n", "In the posterior $p(\\theta \\mid x) \\propto p(x \\mid \\theta)p(\\theta)$, name each term and say which one Day 11 maximised.\n", "\n", "
Show solution\n", "\n", "$p(x \\mid \\theta)$ is the likelihood (Day 11 maximised it), $p(\\theta)$ the prior, and $p(\\theta \\mid x)$ the posterior.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 2 · A flat prior\n", "\n", "With a Beta(1, 1) (uniform) prior and 18 adoptions in 30 households, find the posterior, its mean and its mode.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Posterior Beta(19, 13).\n", "\n", "**Step 2.** Mean $19/32 \\approx 0.594$; mode $18/30 = 0.6$.\n", "\n", "**Step 3.** With a flat prior, the MAP equals the MLE.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 3 · An informative prior\n", "\n", "Prior Beta(5, 5); data 8 successes and 12 failures. Find the posterior mean and mode.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Posterior Beta(13, 17).\n", "\n", "**Step 2.** Mean $13/30 \\approx 0.433$; mode $12/28 \\approx 0.429$.\n", "\n", "**Step 3.** The MLE is $8/20 = 0.40$, pulled toward the prior mean 0.5.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 4 · Posterior predictive\n", "\n", "Using Exercise 3, what is the probability that the next household adopts?\n", "\n", "
Show solution\n", "\n", "$P(\\text{next} = 1 \\mid \\text{data}) = E[p \\mid \\text{data}] = 13/30 \\approx 0.433$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 5 · Prior strength\n", "\n", "How many \"pseudo-observations\" does a Beta(5, 5) prior represent, compared with Beta(50, 50)?\n", "\n", "
Show solution\n", "\n", "$\\alpha + \\beta$: 10 versus 100. Both centre on 0.5, but the second takes ten times as much data to move.\n", "\n", "
\n", "\n", "## Part 2 · Proofs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 6 · Beta–binomial conjugacy\n", "\n", "Prove the posterior is Beta($\\alpha + k$, $\\beta + f$).\n", "\n", "
Show solution\n", "\n", "**Step 1.** Likelihood $\\propto p^k(1 - p)^f$; prior $\\propto p^{\\alpha - 1}(1 - p)^{\\beta - 1}$.\n", "\n", "**Step 2.** Their product is\n", "\n", "$$\n", "p^{\\alpha + k - 1}(1 - p)^{\\beta + f - 1},\n", "$$\n", "\n", "**Step 3.** the kernel of a Beta($\\alpha + k$, $\\beta + f$) density.\n", "\n", "**Step 4.** The normalising constant must be the Beta one, because a density integrates to 1.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 7 · The posterior mean as a weighted average\n", "\n", "Show $E[p \\mid \\text{data}] = w\\cdot\\frac{\\alpha}{\\alpha + \\beta} + (1 - w)\\cdot\\frac{k}{n}$ with $n = k + f$ and $w = \\frac{\\alpha + \\beta}{\\alpha + \\beta + n}$.\n", "\n", "
Show solution\n", "\n", "$$\n", "\\frac{\\alpha + k}{\\alpha + \\beta + n} = \\frac{\\alpha + \\beta}{\\alpha + \\beta + n}\\cdot\\frac{\\alpha}{\\alpha + \\beta} + \\frac{n}{\\alpha + \\beta + n}\\cdot\\frac kn.\n", "$$\n", "\n", "As $n \\to \\infty$, $w \\to 0$ and the posterior mean approaches the MLE: the data dominate.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 8 · Normal–normal updating\n", "\n", "Derive the normal–normal posterior by completing the square.\n", "\n", "
Show solution\n", "\n", "$$\n", "\\ln p(\\theta \\mid \\hat\\theta) = -\\frac{(\\hat\\theta - \\theta)^2}{2s^2} - \\frac{(\\theta - \\mu_0)^2}{2\\tau_0^2} + c = -\\frac{\\pi_d + \\pi_0}{2}\\theta^2 + (\\pi_d\\hat\\theta + \\pi_0\\mu_0)\\theta + c'.\n", "$$\n", "\n", "A quadratic $-\\frac a2\\theta^2 + b\\theta$ is the log-kernel of $N(b/a, 1/a)$. So the posterior is normal with precision $\\pi_0 + \\pi_d$ and mean $(\\pi_0\\mu_0 + \\pi_d\\hat\\theta)/(\\pi_0 + \\pi_d)$.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 9 · MAP with a normal prior is ridge\n", "\n", "In $y = X\\beta + \\varepsilon$ with $\\varepsilon \\sim N(0, \\sigma^2I)$ and prior $\\beta \\sim N(0, \\tau^2I)$, show the MAP minimises $\\|y - X\\beta\\|^2 + \\lambda\\|\\beta\\|^2$ with $\\lambda = \\sigma^2/\\tau^2$.\n", "\n", "
Show solution\n", "\n", "$$\n", "-\\ln p(\\beta \\mid y) = \\frac{1}{2\\sigma^2}\\|y - X\\beta\\|^2 + \\frac{1}{2\\tau^2}\\|\\beta\\|^2 + c.\n", "$$\n", "\n", "Multiply by $2\\sigma^2$: the MAP minimises $\\|y - X\\beta\\|^2 + (\\sigma^2/\\tau^2)\\|\\beta\\|^2$, the ridge objective. A tighter prior (small $\\tau$) means more penalisation.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 10 · Credible versus confidence\n", "\n", "With a flat prior and $\\hat\\theta \\sim N(\\theta, s^2)$, show that the 95% credible interval equals the 95% confidence interval numerically. Do they mean the same thing?\n", "\n", "
Show solution\n", "\n", "**Step 1.** A flat prior gives posterior $\\theta \\mid \\hat\\theta \\sim N(\\hat\\theta, s^2)$, so the credible interval is $\\hat\\theta \\pm 1.96s$, the same numbers as the CI.\n", "\n", "**Step 2.** The meanings differ: the credible interval says \"given the data, $\\theta$ is in here with probability 0.95\"; the CI says \"this procedure covers the true $\\theta$ in 95% of samples\".\n", "\n", "
\n", "\n", "## Part 3 · Applied" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 11 · A sceptical prior\n", "\n", "A policy evaluation reports $\\hat\\tau = 2.5$ (SE 1). Using a sceptical prior $\\tau \\sim N(0, 1)$, find the posterior and $P(\\tau > 0 \\mid \\text{data})$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Precisions: prior $1/1^2 = 1$; data $1/1^2 = 1$.\n", "\n", "**Step 2.** Posterior precision $= 2$, so posterior variance $= 0.5$ and sd $\\approx 0.707$.\n", "\n", "**Step 3.** Posterior mean $= (1 \\times 0 + 1 \\times 2.5)/2 = 1.25$.\n", "\n", "**Step 4.** $P(\\tau > 0 \\mid \\text{data}) = \\Phi(1.25/0.707) = \\Phi(1.77) \\approx 0.96$.\n", "\n", "**Answer:** the sceptical prior halves the estimate to 1.25, but there is still a 96% posterior probability that the policy helps.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 12 · Empirical-Bayes shrinkage\n", "\n", "School effects vary with $\\tau^2 = 4$ around a mean of 0. One school's estimate is 6 with sampling variance 12. Find its shrunken estimate.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Shrinkage factor: $\\tau^2/(\\tau^2 + s^2) = 4/(4 + 12) = 0.25$.\n", "\n", "**Step 2.** Posterior mean $= 0 + 0.25 \\times (6 - 0) = 1.5$.\n", "\n", "**Answer:** 1.5. The school's noisy estimate (sampling variance 12, three times the true spread of school effects) is pulled three-quarters of the way toward the average. Value-added rankings use exactly this adjustment.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 13 · Vaccination uptake: the data dominate\n", "\n", "A health agency's prior for the share of adults who will take a new vaccine is Beta(50, 50). A survey of 1,000 adults finds 600 intend to get it. Compare the posterior mean with the MLE, and give the prior's weight.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Posterior: Beta$(50 + 600,\\; 50 + 400)$ = Beta(650, 450).\n", "\n", "**Step 2.** Posterior mean $= 650/1100 \\approx 0.591$; MLE $= 600/1000 = 0.600$.\n", "\n", "**Step 3.** Prior weight $= (\\alpha + \\beta)/(\\alpha + \\beta + n) = 100/1100 \\approx 9\\%$.\n", "\n", "**Answer:** 59.1% versus 60.0%; the prior carries only 9% weight, and less with more data.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 14 · Peeking at an A/B test\n", "\n", "Why can a frequentist who checks a p-value every day and stops at the first $p < 0.05$ inflate the false-positive rate, and how does a Bayesian analysis treat stopping?\n", "\n", "
Show solution\n", "\n", "**Step 1 (frequentist).** Each daily look is another chance for noise to cross $p < 0.05$, so stopping at the first crossing makes the overall false-positive rate well above 5%. Sequential designs (e.g. alpha-spending) adjust the thresholds.\n", "\n", "**Step 2 (Bayesian).** The posterior depends on the data only through the likelihood, which is the same whatever the stopping rule, so the posterior is unchanged by when you stopped.\n", "\n", "**Step 3 (caveat).** Stopping on a lucky streak still yields optimistic estimates when the prior is weak.\n", "\n", "**Answer:** peeking invalidates fixed-level p-values; Bayesian summaries remain coherent, but a sensible prior is still needed.\n", "\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Your work:*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise 15 · Which vaccination reminder works better?\n", "\n", "A clinic randomises two SMS reminders. Message A: 120 of 1,000 patients book a vaccination; message B: 150 of 1,000. With Beta(1, 1) priors, estimate $P(p_B > p_A \\mid \\text{data})$.\n", "\n", "
Show solution\n", "\n", "**Step 1.** Posteriors: $p_A \\sim$ Beta(121, 881) and $p_B \\sim$ Beta(151, 851).\n", "\n", "**Step 2.** Draw many pairs $(p_A, p_B)$ from these posteriors and count how often $p_B > p_A$ (code below): about **0.975**.\n", "\n", "**Step 3.** Check with a normal approximation. Posterior means ≈ 0.121 and 0.151; sds ≈ 0.0103 and 0.0113. The difference has mean ≈ 0.030 and sd $\\approx \\sqrt{0.0103^2 + 0.0113^2} \\approx 0.0153$, so $P \\approx \\Phi(0.030/0.0153) = \\Phi(1.96) \\approx 0.975$.\n", "\n", "**Answer:** a 97.5% posterior probability that message B works better, with a 95% credible interval for the gain of roughly 0 to 6 points.\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 stats\n", "\n", "rng = np.random.default_rng(13)\n", "a = rng.beta(121, 881, 1_000_000); b = rng.beta(151, 851, 1_000_000)\n", "print(\"P(pB > pA | data) ~\", round((b > a).mean(), 4))\n", "print(\"95% credible interval for pB - pA:\", np.percentile(b - a, [2.5, 97.5]).round(4))\n", "\n", "# Normal-normal update (Exercise 11)\n", "prior_mu, prior_var, est, se2 = 0.0, 1.0, 2.5, 1.0\n", "post_var = 1 / (1 / prior_var + 1 / se2)\n", "post_mu = post_var * (prior_mu / prior_var + est / se2)\n", "print(f\"posterior N({post_mu:.2f}, {post_var:.2f}); P(tau>0) = {stats.norm.sf(0, post_mu, np.sqrt(post_var)):.3f}\")\n", "\n", "# Beta posterior summaries (Exercise 3)\n", "post = stats.beta(13, 17)\n", "print(\"mean\", round(post.mean(), 3), \" 95% credible interval\", np.round(post.ppf([0.025, 0.975]), 3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Key takeaways\n", "\n", "- Posterior ∝ likelihood × prior; conjugate priors give closed-form updates.\n", "- Posterior means are precision-weighted averages, so noisy estimates are shrunk toward the prior.\n", "- With enough data the prior matters little; with little data it matters a lot, so state it openly." ] } ] }