{
 "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 7: The Law of Large Numbers & Central Limit Theorem\n",
    "\n",
    "*Beginner to Ninja: Math for Causal Inference*\n",
    "\n",
    "Why sample means settle down, and why their errors are approximately normal.\n",
    "\n",
    "**Learning outcomes**\n",
    "\n",
    "- Compute the mean, variance and standard error of a sample mean.\n",
    "- Prove Markov's and Chebyshev's inequalities and the weak law of large numbers.\n",
    "- Apply the central limit theorem to sample means, proportions and totals.\n",
    "- Explain how clustering and fat tails change these results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory\n",
    "\n",
    "### The sample mean\n",
    "\n",
    "Let $X_1, \\dots, X_n$ be independent and identically distributed (**iid**) with mean $\\mu$ and variance $\\sigma^2$. The sample mean $\\bar X_n = \\frac1n\\sum_i X_i$ has\n",
    "\n",
    "$$\n",
    "E[\\bar X_n] = \\mu, \\qquad \\operatorname{Var}(\\bar X_n) = \\frac{\\sigma^2}{n}, \\qquad \\operatorname{SE}(\\bar X_n) = \\frac{\\sigma}{\\sqrt n}.\n",
    "$$\n",
    "\n",
    "The **standard error** shrinks like $1/\\sqrt n$: to halve it, you need four times the data.\n",
    "\n",
    "### Two inequalities\n",
    "\n",
    "- **Markov:** if $X \\ge 0$, then $P(X \\ge a) \\le E[X]/a$ for $a > 0$.\n",
    "- **Chebyshev:** $P(|X - \\mu| \\ge k) \\le \\sigma^2/k^2$.\n",
    "\n",
    "### Law of large numbers (LLN)\n",
    "\n",
    "For every $\\varepsilon > 0$, $P(|\\bar X_n - \\mu| > \\varepsilon) \\to 0$ as $n \\to \\infty$. We write $\\bar X_n \\xrightarrow{p} \\mu$ (\"converges in probability\"). This is why a large experiment's group means are close to the population means.\n",
    "\n",
    "### Central limit theorem (CLT)\n",
    "\n",
    "Whatever the shape of the distribution of $X$ (as long as $\\sigma^2$ is finite),\n",
    "\n",
    "$$\n",
    "\\frac{\\bar X_n - \\mu}{\\sigma/\\sqrt n} \\xrightarrow{d} N(0, 1).\n",
    "$$\n",
    "\n",
    "For moderate $n$, $\\bar X_n \\approx N(\\mu, \\sigma^2/n)$. Almost every confidence interval and p-value in applied work rests on this approximation, including those for regression coefficients, which are weighted averages.\n",
    "\n",
    "### When things go wrong\n",
    "\n",
    "- **Dependence.** With $m$ units per cluster and intra-cluster correlation $\\rho$, the variance of the mean is inflated by the **design effect** $1 + (m - 1)\\rho$.\n",
    "- **Infinite variance.** For very fat-tailed data (e.g. Cauchy), neither the LLN nor the CLT applies.\n",
    "\n",
    "## Part 1 · Warm-up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1 · Standard error\n",
    "\n",
    "Test scores have $\\sigma = 20$. What is the standard error of a mean from $n = 100$ students?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$\\operatorname{SE} = 20/\\sqrt{100} = 2$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2 · Diminishing returns\n",
    "\n",
    "How many students are needed to cut the standard error to 1?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$20/\\sqrt n = 1 \\Rightarrow n = 400$: four times as many for half the error.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3 · A probability for the mean\n",
    "\n",
    "$X$ has $\\mu = 50$ and $\\sigma = 10$. For $n = 25$, approximate $P(\\bar X > 53)$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$\\operatorname{SE} = 10/5 = 2$, $z = (53 - 50)/2 = 1.5$, so $P \\approx 1 - \\Phi(1.5) \\approx 0.067$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4 · A proportion\n",
    "\n",
    "The true employment rate is 0.3. For a sample of 400, find the standard error of the sample proportion and an interval containing it with about 95% probability.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** A proportion is the mean of 0/1 variables with $\\sigma^2 = p(1 - p) = 0.3 \\times 0.7 = 0.21$.\n",
    "\n",
    "**Step 2.** Standard error: $\\sqrt{0.21/400} = \\sqrt{0.000525} \\approx 0.0229$.\n",
    "\n",
    "**Step 3.** By the CLT, about 95% of samples fall within $1.96 \\times 0.0229 \\approx 0.045$ of 0.3.\n",
    "\n",
    "**Answer:** SE ≈ 0.023; the sample employment rate lands in $[0.255,\\; 0.345]$ about 95% of the time.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5 · A total\n",
    "\n",
    "An insurer has 100 independent claims with mean \\$500 and sd \\$300. Approximate the probability that total claims exceed \\$55,000.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Total $S = \\sum_{i=1}^{100}X_i$ has mean $100 \\times 500 = 50{,}000$.\n",
    "\n",
    "**Step 2.** Independent claims, so variances add: $\\operatorname{Var}(S) = 100 \\times 300^2$ and $\\operatorname{sd}(S) = 300\\sqrt{100} = 3{,}000$.\n",
    "\n",
    "**Step 3.** By the CLT, $S \\approx N(50{,}000,\\; 3{,}000^2)$. Standardise: $z = (55{,}000 - 50{,}000)/3{,}000 \\approx 1.67$.\n",
    "\n",
    "**Step 4.** $P(S > 55{,}000) \\approx 1 - \\Phi(1.67) \\approx 0.048$.\n",
    "\n",
    "**Answer:** about a 5% chance that total claims exceed \\$55,000.\n",
    "\n",
    "</details>\n",
    "\n",
    "## Part 2 · Proofs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 6 · Unbiasedness of the sample mean\n",
    "\n",
    "Prove $E[\\bar X_n] = \\mu$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "By linearity, $E[\\bar X_n] = \\frac1n\\sum_i E[X_i] = \\frac1n(n\\mu) = \\mu$. (Identical distribution is needed; independence is not.)\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 7 · Variance of the sample mean\n",
    "\n",
    "Prove $\\operatorname{Var}(\\bar X_n) = \\sigma^2/n$ for iid data.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "$$\n",
    "\\operatorname{Var}(\\bar X_n) = \\frac{1}{n^2}\\operatorname{Var}\\Big(\\sum_i X_i\\Big) = \\frac{1}{n^2}\\Big(\\sum_i \\operatorname{Var}(X_i) + \\sum_{i \\ne j}\\operatorname{Cov}(X_i, X_j)\\Big) = \\frac{1}{n^2}\\,n\\sigma^2 = \\frac{\\sigma^2}{n},\n",
    "$$\n",
    "\n",
    "because independence makes every covariance zero. Exercise 13 shows what happens when it does not.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 8 · Markov's inequality\n",
    "\n",
    "Prove $P(X \\ge a) \\le E[X]/a$ for $X \\ge 0$ and $a > 0$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Since $X \\ge 0$, $X \\ge a\\,\\mathbf{1}\\{X \\ge a\\}$ always holds.\n",
    "\n",
    "**Step 2.** Take expectations: $E[X] \\ge a\\,P(X \\ge a)$.\n",
    "\n",
    "**Step 3.** Divide by $a$.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 9 · Chebyshev's inequality\n",
    "\n",
    "Derive Chebyshev's inequality from Markov's.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "Apply Markov to the non-negative variable $(X - \\mu)^2$ with $a = k^2$:\n",
    "\n",
    "$$\n",
    "P(|X - \\mu| \\ge k) = P\\big((X - \\mu)^2 \\ge k^2\\big) \\le \\frac{E[(X - \\mu)^2]}{k^2} = \\frac{\\sigma^2}{k^2}.\n",
    "$$\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 10 · The weak law of large numbers\n",
    "\n",
    "Prove $\\bar X_n \\xrightarrow{p} \\mu$ when $\\sigma^2 < \\infty$.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "Apply Chebyshev to $\\bar X_n$, whose variance is $\\sigma^2/n$:\n",
    "\n",
    "$$\n",
    "P(|\\bar X_n - \\mu| \\ge \\varepsilon) \\le \\frac{\\sigma^2}{n\\varepsilon^2} \\to 0 \\quad\\text{as } n \\to \\infty.\n",
    "$$\n",
    "\n",
    "</details>\n",
    "\n",
    "## Part 3 · Applied"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 11 · Budgeting an RCT\n",
    "\n",
    "An outcome has $\\sigma = 30$ in both arms. (a) With 225 per arm, what is the standard error of the difference in means? (b) How many per arm are needed for a standard error of 2?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**(a) SE with 225 per arm**\n",
    "\n",
    "**Step 1.** Each arm mean has variance $\\sigma^2/n = 900/225 = 4$.\n",
    "\n",
    "**Step 2.** The arms are independent, so the variance of the difference is $4 + 4 = 8$ and $\\operatorname{SE} = \\sqrt 8 \\approx 2.83$.\n",
    "\n",
    "**(b) Sample size for SE = 2**\n",
    "\n",
    "**Step 3.** Solve $\\sqrt{2 \\times 900/n} = 2$: square both sides, $1800/n = 4$, so $n = 450$ per arm.\n",
    "\n",
    "**Answer:** SE ≈ 2.83 with 225 per arm; 450 per arm are needed for SE = 2 (twice the sample for a 29% smaller SE).\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 12 · How sure is a survey mean?\n",
    "\n",
    "A statistics office estimates mean household income from a random sample of $n$ households. Bound the probability that the sample mean misses the true mean by at least two standard errors, $P(|\\bar X_n - \\mu| \\ge 2\\sigma/\\sqrt n)$, using Chebyshev's inequality, and compare with the CLT approximation. Why might the office prefer the Chebyshev figure for very skewed incomes in a small sample?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Chebyshev with $\\operatorname{Var}(\\bar X_n) = \\sigma^2/n$ and $k = 2\\sigma/\\sqrt n$:\n",
    "\n",
    "$$\n",
    "P(|\\bar X_n - \\mu| \\ge k) \\le \\frac{\\sigma^2/n}{4\\sigma^2/n} = 0.25.\n",
    "$$\n",
    "\n",
    "**Step 2.** CLT: $(\\bar X_n - \\mu)/(\\sigma/\\sqrt n) \\approx N(0, 1)$, so the probability is about $2(1 - \\Phi(2)) \\approx 0.046$.\n",
    "\n",
    "**Answer:** at most 25% (Chebyshev) versus about 4.6% (CLT). Chebyshev holds for any distribution and any $n$, so it is a safe worst case when incomes are very skewed and the sample small; the CLT figure is sharper but only an approximation.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 13 · Clustered samples\n",
    "\n",
    "A school survey samples 100 schools with 20 pupils each. The intra-class correlation is 0.05. Find the design effect and the effective sample size.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** Design effect: $\\text{DEFF} = 1 + (m - 1)\\rho = 1 + 19 \\times 0.05 = 1.95$.\n",
    "\n",
    "**Step 2.** Effective sample size: $2{,}000/1.95 \\approx 1{,}026$ independent pupils.\n",
    "\n",
    "**Step 3.** Standard errors that ignore clustering are too small by $\\sqrt{1.95} \\approx 1.40$.\n",
    "\n",
    "**Answer:** DEFF = 1.95; the 2,000 pupils carry the information of about 1,026 independent ones.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 14 · When the LLN fails: very fat tails\n",
    "\n",
    "Suppose the profits of start-ups in a venture-capital portfolio follow a Cauchy distribution (no finite mean or variance). Explain why the average profit across more and more start-ups does not settle down, and what this suggests for heavy-tailed economic outcomes such as wealth or firm revenue.\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** The LLN and CLT both need a finite mean (and the CLT a finite variance). The Cauchy has neither.\n",
    "\n",
    "**Step 2.** A known property: the average of $n$ independent Cauchy draws is itself Cauchy with the *same* spread, so averaging never reduces uncertainty.\n",
    "\n",
    "**Step 3.** Intuition: every so often one start-up has a gigantic gain or loss that dominates the whole average (see the running means in the code).\n",
    "\n",
    "**Answer:** the portfolio average does not converge. Real outcomes such as wealth have finite but huge variances, so convergence is slow and a few extreme units drive the mean; analysts use logs, winsorising, medians or quantiles.\n",
    "\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Your work:*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 15 · Seeing the CLT in hospital stays\n",
    "\n",
    "Hospital length of stay is strongly right-skewed: suppose it is exponential with mean 1 week. A health authority audits random samples of 50 patients and records each sample's mean stay. If you simulate 10,000 such audits, what mean, variance and shape should the sample means have?\n",
    "\n",
    "<details><summary>Show solution</summary>\n",
    "\n",
    "**Step 1.** For the exponential with mean 1, $\\mu = 1$ and $\\sigma^2 = 1$.\n",
    "\n",
    "**Step 2.** $E[\\bar X] = \\mu = 1$ week.\n",
    "\n",
    "**Step 3.** $\\operatorname{Var}(\\bar X) = \\sigma^2/n = 1/50 = 0.02$, so sd ≈ 0.141 weeks.\n",
    "\n",
    "**Step 4.** Shape: by the CLT, approximately normal. The skewness shrinks from 2 for one patient to about $2/\\sqrt{50} \\approx 0.28$ for the mean.\n",
    "\n",
    "**Answer:** mean 1, variance 0.02, nearly bell-shaped despite the skewed stays (see the code and the CLT simulator).\n",
    "\n",
    "</details>\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(2026)\n",
    "means = rng.exponential(1.0, size=(10_000, 50)).mean(axis=1)\n",
    "print(f\"mean of means {means.mean():.3f} (theory 1), var {means.var():.4f} (theory 0.02)\")\n",
    "print(f\"skewness of draws ~2, of means {stats.skew(means):.2f} (theory 2/sqrt(50) = 0.28)\")\n",
    "\n",
    "print(\"P(Xbar > 53):\", round(stats.norm.sf(53, 50, 2), 4))\n",
    "print(\"Total claims tail:\", round(stats.norm.sf(55_000, 50_000, 3_000), 4))\n",
    "print(\"Chebyshev bound 0.25 vs CLT:\", round(2 * stats.norm.sf(2), 4))\n",
    "\n",
    "# Cauchy: running mean never settles\n",
    "c = rng.standard_cauchy(1_000_000)\n",
    "print(\"Cauchy running means:\", [round(float(c[:k].mean()), 2) for k in (10**3, 10**4, 10**5, 10**6)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Key takeaways\n",
    "\n",
    "- Sample means are unbiased with standard error $\\sigma/\\sqrt n$.\n",
    "- The LLN says they converge to $\\mu$; the CLT says their errors are approximately normal.\n",
    "- Dependence (clusters) and fat tails weaken both results."
   ]
  }
 ]
}