{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "---\n",
    "execute:\n",
    "  enabled: false\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-1 fw=all -->\n",
    "\n",
    "# Statistics\n",
    "\n",
    "A trained model is only ever fit to a *finite* sample, so every quantity we read off it (an accuracy, a learned weight, an estimated mean) is a guess computed from random data and would come out differently on a fresh draw. Statistics is the discipline that quantifies that randomness: it tells us how far a guess typically sits from the truth, when an apparent improvement is real rather than noise, and how confident we are entitled to be. This section develops the three ideas a deep-learning practitioner reaches for most often. We define an *estimator* and the two ways it can be wrong, *bias* and *variance*, and prove the decomposition that ties them together; this single identity is behind the U-curve that governs under- and over-fitting in that section. We then turn to *hypothesis testing*, the framework behind A/B tests and benchmark comparisons, and close with *confidence intervals*, which attach a notion of uncertainty to a point estimate. Throughout we take the true parameter $\\theta$ to be a scalar; the vector case is identical with sums of squares replaced by squared norms.\n",
    "\n",
    "We use the following imports throughout the section, plus plain NumPy as `onp` for the label-shuffling in the permutation test and the resampling in the bootstrap."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "statistics-imports"
   },
   "source": [
    "%matplotlib inline\n",
    "from d2l import torch as d2l\n",
    "import torch\n",
    "import numpy as onp  # plain NumPy for the permutation test and bootstrap below"
   ],
   "execution_count": 1,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-2 fw=all -->\n",
    "\n",
    "## Estimators and Their Quality\n",
    "\n",
    "### Estimators\n",
    "\n",
    "An *estimator* is a rule that turns data into a guess for an unknown parameter. Formally, given samples $x_1,\\ldots,x_n$ drawn from a distribution governed by a parameter $\\theta$, an estimator is a function\n",
    "\n",
    "$$\n",
    "\\hat\\theta_n = \\hat f(x_1,\\ldots,x_n)\n",
    "$$\n",
    "\n",
    "that we hope lands near $\\theta$. We have met estimators already: in that section the maximum-likelihood estimate of a Bernoulli probability was the fraction of observed ones, and the maximum-likelihood estimate of a Gaussian mean was the sample average. The key fact is that $\\hat\\theta_n$ is itself a *random variable*: it depends on the random sample, so it would come out differently on a fresh dataset. Asking whether an estimator is *good* is therefore asking about the distribution of $\\hat\\theta_n$ over repeated datasets (its *sampling distribution*), and that distribution has two features that matter, its center and its spread.\n",
    "\n",
    "### Bias and Variance\n",
    "\n",
    "The first feature is the *center*. The **bias** of $\\hat\\theta_n$ measures the systematic gap between where the estimator centers and the truth,\n",
    "\n",
    "$$\n",
    "\\operatorname{Bias}(\\hat\\theta_n) = \\mathbb{E}[\\hat\\theta_n] - \\theta ,\n",
    "$$\n",
    "\n",
    "the expectation taken over the random sample. When $\\operatorname{Bias}(\\hat\\theta_n)=0$ for every $\\theta$ we call $\\hat\\theta_n$ *unbiased*: it is right *on average*, even though any single estimate misses. Bias is systematic error at a fixed sample size. It may persist as data grows,\n",
    "or it may vanish for an asymptotically unbiased estimator; that distinction is\n",
    "made below.\n",
    "\n",
    "The second feature is the *spread*. The **variance** measures how much the estimator fluctuates around its own center, with the **standard error** its square root,\n",
    "\n",
    "$$\n",
    "\\operatorname{Var}(\\hat\\theta_n) = \\mathbb{E}\\!\\left[(\\hat\\theta_n - \\mathbb{E}[\\hat\\theta_n])^2\\right],\n",
    "\\qquad\n",
    "\\operatorname{se}(\\hat\\theta_n) = \\sqrt{\\operatorname{Var}(\\hat\\theta_n)} .\n",
    "$$\n",
    "\n",
    "Note carefully that variance is measured against $\\mathbb{E}[\\hat\\theta_n]$, *not* against the true $\\theta$: it captures the noise in the estimator, not its accuracy. the figure makes the two features visible by drawing the sampling distribution for two estimators of the same $\\theta$. Bias is the offset of the distribution's center from $\\theta$; variance is its width.\n",
    "\n",
    "![The sampling distribution of an estimator, drawn over many datasets. Bias is the offset of the center $\\mathbb{E}\\left(\\hat\\theta\\right)$ from the true $\\theta$; variance is the spread. Left: low bias, high variance. Right: high bias, low variance. A good estimator wants both small, but usually one is traded for the other.](../img/mdl-prob-sampling-distribution.svg)\n",
    "\n",
    "### Consistency and Efficiency\n",
    "\n",
    "Bias and variance describe an estimator at a *fixed* sample size. Two further notions describe how it behaves as data accumulates. An estimator is *asymptotically unbiased* if its bias vanishes in the limit, $\\lim_{n\\to\\infty}\\operatorname{Bias}(\\hat\\theta_n)=0$; many estimators used in practice are biased at finite $n$ but asymptotically unbiased, which is usually good enough. A stronger and more useful guarantee is *consistency*: $\\hat\\theta_n$ is consistent if it *converges in probability* to $\\theta$,\n",
    "\n",
    "$$\n",
    "\\hat\\theta_n \\xrightarrow{P} \\theta,\n",
    "\\qquad\\textrm{i.e.}\\qquad\n",
    "P\\bigl(|\\hat\\theta_n-\\theta|>\\varepsilon\\bigr)\\to 0 \\quad\\textrm{for every } \\varepsilon>0 .\n",
    "$$\n",
    "\n",
    "Consistency is the formal content of the slogan \"more data gets us arbitrarily close to the truth.\" Its prototype is the *weak law of large numbers*: the sample mean is a consistent estimator of the population mean, a proposition we state and prove below once the bias-variance decomposition supplies the tools. We have also already proved a headline instance: that section established consistency for the maximum-likelihood estimator under its regularity hypotheses. A sufficient condition for consistency in general is that *both* the bias and the variance tend to zero, since then the whole sampling distribution collapses onto $\\theta$. (The two limits are independent: an estimator can be asymptotically unbiased yet inconsistent because its variance does not vanish, and its variance can vanish while a persistent bias keeps it inconsistent.)\n",
    "\n",
    "Finally, among *unbiased* estimators we prefer the one that fluctuates least, and we call it *efficient*: efficiency ranks unbiased estimators by their variance, the smaller the better. There is a hard floor here, and we have met it too: the Cramér--Rao bound of that section places a lower bound on the variance of any unbiased estimator, $1/(n\\,I(\\theta))$, the inverse Fisher information, and the coin-flip simulation there watched the maximum-likelihood estimator land right on that floor. An estimator attaining the bound is as good as unbiased estimation can be (the precise sense in which the MLE is asymptotically efficient). What this section takes from the bound: once unbiasedness is secured, what remains to minimize is variance, which is exactly the second half of the decomposition we turn to next.\n",
    "\n",
    "## The Bias-Variance Decomposition\n",
    "\n",
    "We now have two distinct ways an estimator can be wrong, a systematic offset (bias) and random fluctuation (variance), and a single number that ought to combine them: the *mean squared error*. The MSE is *exactly* the sum of these two contributions, with no cross term.\n",
    "\n",
    "### Mean Squared Error and the Decomposition\n",
    "\n",
    "The simplest summary of how far an estimator lands from the truth is the **mean squared error**,\n",
    "\n",
    "$$\n",
    "\\operatorname{MSE}(\\hat\\theta_n) = \\mathbb{E}\\!\\left[(\\hat\\theta_n-\\theta)^2\\right] .\n",
    "$$\n",
    "\n",
    "It is always non-negative, and the smaller it is the closer $\\hat\\theta_n$ sits to $\\theta$ on average. If you have read that section you will recognize it as the squared-error loss, now applied to an estimator rather than a prediction.\n",
    "\n",
    "**Proposition (bias-variance decomposition).** *For any estimator $\\hat\\theta_n$ of a fixed parameter $\\theta$,*\n",
    "\n",
    "$$\n",
    "\\operatorname{MSE}(\\hat\\theta_n) = \\operatorname{Bias}(\\hat\\theta_n)^2 + \\operatorname{Var}(\\hat\\theta_n) .\n",
    "$$\n",
    "\n",
    "**Proof.** Abbreviate $\\mu = \\mathbb{E}[\\hat\\theta_n]$, the center of the estimator, and add and subtract it inside the square:\n",
    "\n",
    "$$\n",
    "\\operatorname{MSE}(\\hat\\theta_n)\n",
    " = \\mathbb{E}\\!\\left[(\\hat\\theta_n - \\theta)^2\\right]\n",
    " = \\mathbb{E}\\!\\left[\\bigl((\\hat\\theta_n - \\mu) + (\\mu - \\theta)\\bigr)^2\\right].\n",
    "$$\n",
    "\n",
    "Expanding the square gives three terms. The first is $\\mathbb{E}[(\\hat\\theta_n-\\mu)^2]=\\operatorname{Var}(\\hat\\theta_n)$; the last is $(\\mu-\\theta)^2=\\operatorname{Bias}(\\hat\\theta_n)^2$, a constant. The middle, cross term *vanishes*, because $\\mu-\\theta$ is a constant and $\\hat\\theta_n-\\mu$ has mean zero by the definition of $\\mu$:\n",
    "\n",
    "$$\n",
    "2\\,(\\mu-\\theta)\\,\\mathbb{E}[\\hat\\theta_n - \\mu]\n",
    " = 2\\,(\\mu-\\theta)\\,(\\mu - \\mu) = 0 .\n",
    "$$\n",
    "\n",
    "What remains is $\\operatorname{Var}(\\hat\\theta_n)+\\operatorname{Bias}(\\hat\\theta_n)^2$. $\\blacksquare$\n",
    "\n",
    "Because the deviation from the center has mean zero, the systematic part and the fluctuating part of the error never interfere, and the squared error splits into the two pieces of the figure. One immediate payoff is the consistency criterion promised above: if both $\\operatorname{Bias}(\\hat\\theta_n)\\to0$ and $\\operatorname{Var}(\\hat\\theta_n)\\to0$, then the equation forces $\\operatorname{MSE}(\\hat\\theta_n)\\to0$, which implies $\\hat\\theta_n\\xrightarrow{P}\\theta$: applying Markov's inequality the equation from that section to the nonnegative random variable $(\\hat\\theta_n-\\theta)^2$ gives $P(|\\hat\\theta_n-\\theta|>\\varepsilon) \\le \\operatorname{MSE}(\\hat\\theta_n)/\\varepsilon^2 \\to 0$. (For an unbiased estimator this is Chebyshev's inequality the equation; aiming Markov at $\\theta$ directly covers the biased case too. Sharper, exponentially decaying rates come from the concentration inequalities of that section when the data are bounded or sub-Gaussian.)\n",
    "\n",
    "### The Law of Large Numbers\n",
    "\n",
    "Assembling the pieces just derived yields the guarantee promised above, the same one that that section invoked to justify Monte Carlo integration.\n",
    "\n",
    "**Proposition (weak law of large numbers).** *If $x_1,\\ldots,x_n$ are i.i.d. with mean $\\mu$ and finite variance $\\sigma^2$, the sample mean $\\bar x_n = \\frac1n\\sum_{i=1}^n x_i$ converges in probability to $\\mu$:*\n",
    "\n",
    "$$\n",
    "\\bar x_n \\xrightarrow{P} \\mu ,\n",
    "\\qquad\\textrm{indeed}\\qquad\n",
    "P\\bigl(|\\bar x_n - \\mu| > \\varepsilon\\bigr) \\le \\frac{\\sigma^2}{n\\,\\varepsilon^2}\n",
    "\\quad\\textrm{for every } \\varepsilon > 0 .\n",
    "$$\n",
    "\n",
    "**Proof.** The sample mean is unbiased, $\\mathbb{E}[\\bar x_n]=\\mu$, by linearity of expectation, and its variance is $\\operatorname{Var}(\\bar x_n)=\\sigma^2/n$ by the variance of a sum of independent variables (that section). By the equation its MSE is therefore $0 + \\sigma^2/n$, and the Markov step above turns that into the displayed bound, which vanishes as $n\\to\\infty$. $\\blacksquare$\n",
    "\n",
    "Finite variance is more than the theorem needs: a finite mean $\\mathbb{E}|x_1|<\\infty$ already suffices, at the cost of a longer argument [@Wasserman.2013]. In whichever form, this is the law that makes averaging work: it is why a Monte Carlo estimate of an integral converges, why a test-set accuracy estimates the true error rate, and why the empirical loss a model is trained on tracks the expected loss it is meant to minimize.\n",
    "\n",
    "### The Trade-off and Generalization\n",
    "\n",
    "Identity the equation also explains the central tension of model fitting, told in full in that section. Read $\\hat\\theta_n$ as a *fitted model* and $\\theta$ as the function we wish it had learned: then *underfitting* is bias (a too-simple model misses systematically, however much data it sees) and *overfitting* is variance (a too-flexible model chases the noise of its particular training set) [@Geman.Bienenstock.Doursat.1992]. As model complexity grows, the squared bias falls while the variance rises, and their sum traces the U of the figure; for prediction the expected test error acquires one more term, $\\operatorname{Bias}^2 + \\operatorname{Var} + \\sigma^2$, where $\\sigma^2$ is the *irreducible noise* in the labels, a constant floor that shifts the U upward without moving its minimum. Regularization such as weight decay (that section) trades a little added bias for a large reduction in variance, sliding leftward on the curve to a lower total error. (Heavily overparameterized deep networks can defy this textbook U; the decomposition remains exactly true, and that section reproduces the resulting double-descent curve from scratch.)\n",
    "\n",
    "![As model complexity grows, squared bias falls and variance rises; their sum, the MSE (test error), is a U-curve with a minimum at the sweet spot.](../img/mdl-prob-bias-variance-u-curve.svg)\n",
    "\n",
    "### The Decomposition in Code\n",
    "\n",
    "The decomposition is an exact algebraic identity, so it should hold to numerical precision on a concrete example. We first define bias and MSE as the formulas the equation and the equation say: averages over a collection of estimates, against the true parameter."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "statistics-estimator-metrics"
   },
   "source": [
    "import numpy as onp\n",
    "def stat_bias(true_theta, est_theta):  # E[theta_hat] - theta\n",
    "    return onp.mean(est_theta) - true_theta\n",
    "\n",
    "def mse(est_theta, true_theta):        # E[(theta_hat - theta)^2]\n",
    "    return onp.mean(onp.square(est_theta - true_theta))"
   ],
   "execution_count": 2,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-3 fw=all -->\n",
    "\n",
    "To exercise these we need the *sampling distribution* itself, not a single dataset: we draw many independent datasets from $\\mathcal{N}(\\theta,\\sigma^2)$, compute the sample mean on each, and collect the resulting estimates. Their spread is the variance and their center the bias."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "statistics-sampling-distribution"
   },
   "source": [
    "import numpy as onp\n",
    "theta_true, sigma = 1.0, 4.0\n",
    "num_datasets, n = 10000, 30  # 10k datasets, each of n=30 points\n",
    "samples = onp.random.normal(theta_true, sigma, (num_datasets, n))\n",
    "theta_hats = samples.mean(axis=1)  # one sample-mean estimate per dataset"
   ],
   "execution_count": 3,
   "outputs": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-4 fw=all -->\n",
    "\n",
    "Now we read the decomposition off the empirical sampling distribution: the MSE of the estimates around the true $\\theta$ should match the squared bias plus the variance of the estimates around their own mean, the two sides of the equation. One detail matters for exactness. The identity is a statement about expectations under a *single* distribution, here the empirical distribution of our $10{,}000$ estimates, whose expectations are plain averages; that dictates the *plug-in* variance that divides by the number of estimates (`ddof=0`, the default in most libraries), not the unbiased $n-1$ variant we meet in the next subsection. With the plug-in choice the proof's algebra goes through verbatim for the empirical averages."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "statistics-verify-decomposition"
   },
   "source": [
    "import numpy as onp\n",
    "bias = stat_bias(theta_true, theta_hats)\n",
    "# Default ddof=0: the plug-in variance of the estimates, which makes the\n",
    "# identity exact for empirical averages (the n-1 variant is the next subsection)\n",
    "var = onp.var(theta_hats)\n",
    "mse(theta_hats, theta_true), var + onp.square(bias)"
   ],
   "execution_count": 4,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(np.float64(0.5284503820141454), np.float64(0.5284503820141454))"
     },
     "metadata": {},
     "execution_count": 4
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-5 fw=all -->\n",
    "\n",
    "The two numbers agree to floating-point round-off (the identity is exact), and both are close to the theoretical value. For the sample mean of $\\mathcal{N}(\\theta,\\sigma^2)$ the bias is exactly zero (the average of unbiased draws is unbiased) and the variance is $\\sigma^2/n$ (the variance-of-a-sum result from that section), so $\\operatorname{MSE}=\\sigma^2/n = 16/30 \\approx 0.53$. Because *both* the bias ($0$) and the variance ($\\sigma^2/n\\to0$) vanish as $n\\to\\infty$, the sample mean is consistent, exactly as the weak law of large numbers the equation asserts: the simulation is that proposition run in the small.\n",
    "\n",
    "### Why the Unbiased Variance Divides by $n-1$\n",
    "\n",
    "The sample mean was unbiased for free. The sample *variance* is more delicate, and it exposes a subtlety that the library's `std` function encodes in a `ddof` flag. Given samples $x_1,\\ldots,x_n$ with sample mean $\\bar x=\\frac1n\\sum_i x_i$, the natural estimator of the population variance $\\sigma^2$ would average the squared deviations,\n",
    "\n",
    "$$\n",
    "s_0^2 = \\frac1n\\sum_{i=1}^n (x_i-\\bar x)^2 .\n",
    "$$\n",
    "\n",
    "This is *biased*: it systematically underestimates $\\sigma^2$, because the deviations are measured from $\\bar x$ (the point that *minimizes* the sum of squared deviations for this particular sample) rather than from the unknown true mean $\\mu$. The $n$-divided estimator matters in its own right: $s_0^2$ *is* the maximum-likelihood estimator of a Gaussian's variance, so the bias computed here is precisely the finite-sample bias of the MLE that that section promised when it called maximum likelihood only *asymptotically* unbiased. The fix is to divide by $n-1$ instead of $n$, and the factor is exactly what unbiasedness requires.\n",
    "\n",
    "**Proposition (unbiased sample variance).** *For i.i.d. samples with variance $\\sigma^2$,*\n",
    "\n",
    "$$\n",
    "s^2 = \\frac{1}{n-1}\\sum_{i=1}^n (x_i-\\bar x)^2\n",
    "\\qquad\\textrm{satisfies}\\qquad\n",
    "\\mathbb{E}[s^2] = \\sigma^2 .\n",
    "$$\n",
    "\n",
    "**Proof.** Center the data at the true mean $\\mu$ by writing $x_i-\\bar x = (x_i-\\mu)-(\\bar x-\\mu)$, and expand the sum of squared deviations:\n",
    "\n",
    "$$\n",
    "\\sum_{i=1}^n (x_i-\\bar x)^2\n",
    " = \\sum_{i=1}^n (x_i-\\mu)^2 - n\\,(\\bar x-\\mu)^2 ,\n",
    "$$\n",
    "\n",
    "where the cross term collapsed because $\\sum_i (x_i-\\mu) = n(\\bar x-\\mu)$. Now take expectations. Each $\\mathbb{E}[(x_i-\\mu)^2]=\\sigma^2$, so the first sum has expectation $n\\sigma^2$. The second uses the variance of the sample mean, $\\mathbb{E}[(\\bar x-\\mu)^2]=\\operatorname{Var}(\\bar x)=\\sigma^2/n$, so that term has expectation $n\\cdot\\sigma^2/n=\\sigma^2$. Hence\n",
    "\n",
    "$$\n",
    "\\mathbb{E}\\!\\left[\\sum_{i=1}^n (x_i-\\bar x)^2\\right] = n\\sigma^2 - \\sigma^2 = (n-1)\\,\\sigma^2 .\n",
    "$$\n",
    "\n",
    "Dividing by $n-1$ gives $\\mathbb{E}[s^2]=\\sigma^2$. $\\blacksquare$\n",
    "\n",
    "The intuition is *degrees of freedom*: estimating $\\bar x$ from the same data consumes one degree of freedom, so only $n-1$ of the deviations are free to vary, and dividing by $n-1$ rather than $n$ corrects for it exactly. (As $n\\to\\infty$ the two estimators agree, so $s_0^2$ is biased but asymptotically unbiased and consistent.) We can watch the bias appear and the correction remove it by estimating both variances over many datasets and averaging."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "statistics-unbiased-variance"
   },
   "source": [
    "import numpy as onp\n",
    "n = 3  # small n makes the 1/n bias glaring; the gap shrinks like 1/n\n",
    "data = onp.random.normal(0, 2, (100000, n))  # sigma^2 = 4\n",
    "dev2 = onp.square(data - data.mean(axis=1, keepdims=True)).sum(axis=1)\n",
    "print('true variance    = 4')\n",
    "print(f'E[divide by n]   = {float((dev2 / n).mean()):.3f}  (biased)')\n",
    "print(f'E[divide by n-1] = {float((dev2 / (n - 1)).mean()):.3f}  (unbiased)')"
   ],
   "execution_count": 5,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "true variance    = 4\nE[divide by n]   = 2.666  (biased)\nE[divide by n-1] = 3.998  (unbiased)\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-6 fw=all -->\n",
    "\n",
    "With $n=3$ the biased estimator averages near $\\tfrac{n-1}{n}\\sigma^2 = \\tfrac23\\cdot 4 \\approx 2.67$, while dividing by $n-1$ recovers $4$, confirming the equation. The `ddof` argument (\"delta degrees of freedom\") selects the denominator: `ddof=1` divides by $n-1$, `ddof=0` by $n$.\n",
    "\n",
    "## Hypothesis Testing\n",
    "\n",
    "The bias-variance picture asks *how good* a single estimate is. Hypothesis testing asks a different question that dominates experimental practice: given two estimates (a baseline and a new model, a control group and a treatment), is the observed difference *real*, or could it be a fluke of the particular sample? This is the framework behind A/B testing and behind claims that one architecture beats another on a benchmark.\n",
    "\n",
    "### The Setup: Null, Alternative, and Two Kinds of Error\n",
    "\n",
    "A *hypothesis test* weighs evidence against a default claim. The **null hypothesis** $H_0$ is that default, typically \"there is no effect\" (e.g. the new model is no better than the baseline), and the **alternative** $H_A$ is the effect we hope to detect: sometimes the null's outright negation, but often one-sided or otherwise composite, e.g. \"the new model is *better*.\" The asymmetry is deliberate: we never *prove* $H_0$; we either gather enough evidence to *reject* it in favor of $H_A$, or we fail to, much as a court returns \"guilty\" or \"not guilty\" rather than \"innocent.\"\n",
    "\n",
    "Because the data are random, the decision can go wrong in two ways. A **type I error** (false positive) is rejecting $H_0$ when it is in fact true: declaring an effect that is not there. A **type II error** (false negative) is failing to reject $H_0$ when it is in fact false: missing a real effect. Their rates have standard names, the significance level $\\alpha$ and $\\beta$,\n",
    "\n",
    "$$\n",
    "\\alpha = P(\\textrm{reject } H_0 \\mid H_0 \\textrm{ true}),\n",
    "\\qquad\n",
    "\\beta = P(\\textrm{fail to reject } H_0 \\mid H_0 \\textrm{ false}),\n",
    "$$\n",
    "\n",
    "and the four possible outcomes arrange into the $2\\times2$ decision matrix of the figure: rows are whether $H_0$ is true or false, columns are our decision. The diagonal cells are correct; the off-diagonal cells are the two errors. One caveat on the notation: when $H_A$ is composite, \"$H_0$ false\" does not pin down a single distribution, so $\\beta$ (and the power below) is a *function* of the true parameter, evaluated at a given alternative; that is exactly how the figure will draw it.\n",
    "\n",
    "![The $2\\times 2$ hypothesis-test decision matrix: correct decisions on the diagonal, the type I error (rate $\\alpha$) and type II error (rate $\\beta$) off the diagonal, and the power $1-\\beta$ in the correct-rejection cell.](../img/mdl-prob-type-i-ii-matrix.svg)\n",
    "\n",
    "### Significance and Power\n",
    "\n",
    "We *choose* the type I error rate up front: the **significance level** $\\alpha$ is the risk of a false positive we are willing to tolerate, conventionally $\\alpha=0.05$. The complement $1-\\alpha$ is the *confidence level*, the name that returns with the confidence intervals below. The bottom-right cell of the matrix is the quantity we want to be large: the **statistical power**\n",
    "\n",
    "$$\n",
    "1 - \\beta = P(\\textrm{reject } H_0 \\mid H_0 \\textrm{ false})\n",
    "$$\n",
    "\n",
    "is the probability the test *detects* a real effect, again evaluated at a given alternative. A test with $\\alpha=0.05$ but power $0.2$ rejects a true null only $5\\%$ of the time yet still misses $80\\%$ of genuine effects of that size: underpowered, and useless for confirming improvements. A common target is $1-\\beta=0.8$.\n",
    "\n",
    "Power is what determines how much data we need, and the simplest test of all makes the accounting explicit. In the one-sample two-sided **$z$-test** the data are Gaussian with *known* variance $\\sigma^2$, the null is $H_0\\!:\\mu=0$, the statistic is $T=\\sqrt n\\,\\bar x/\\sigma$, which under $H_0$ is exactly standard normal, and we reject when $|T|\\ge z_{1-\\alpha/2}$, where $z_q$ denotes the standard-normal quantile ($z_{0.975}\\approx1.960$ at $\\alpha=0.05$). Now suppose the true mean is some $\\mu\\neq0$ and write $\\delta=|\\mu|/\\sigma$ for the **standardized effect**, the size of the effect measured in units of the noise. Then $T$ is Gaussian with mean $\\pm\\delta\\sqrt n$ and unit variance, so the test detects the effect with probability $\\Phi(\\delta\\sqrt n - z_{1-\\alpha/2})$, where $\\Phi$ is the standard-normal c.d.f. (the far tail on the opposite side contributes negligibly). Setting this power equal to $1-\\beta$, so that $\\delta\\sqrt n - z_{1-\\alpha/2} = z_{1-\\beta}$, and solving for $n$ gives the sample size needed to detect a standardized effect $\\delta$ at level $\\alpha$ with power $1-\\beta$:\n",
    "\n",
    "$$\n",
    "n \\;=\\; \\frac{\\bigl(z_{1-\\alpha/2} + z_{1-\\beta}\\bigr)^{2}}{\\delta^{2}} .\n",
    "$$\n",
    "\n",
    "At $\\alpha=0.05$ and power $0.8$ we have $z_{0.975}\\approx1.960$ and $z_{0.8}\\approx0.842$, so $n\\approx(2.80)^2/\\delta^2 = 7.85/\\delta^2$: about $8$ samples suffice for a large effect $\\delta=1$, while a tiny effect $\\delta=0.01$ needs $78{,}489$. This is the $1/\\delta^2$ law with its constant filled in, the formula behind exercise 6, and the reason marginal benchmark gains demand enormous test sets to confirm. the figure traces the whole family of *power curves*: for each effect size $\\delta$, the probability of detection climbs from $\\alpha$ (a false-positive rate is all a test delivers at $\\delta=0$) toward $1$ as $n$ grows, crossing the conventional target $0.8$ at a sample size proportional to $1/\\delta^2$.\n",
    "\n",
    "![Power of the one-sample two-sided $z$-test at $\\alpha=0.05$ as a function of the sample size $n$, one curve per effect size $\\delta$. Every curve starts near $\\alpha$ and climbs toward $1$; the dashed line marks the conventional target $0.8$, reached at a sample size that scales like $1/\\delta^2$: about $8$ samples for $\\delta=1$ but nearly $80{,}000$ for $\\delta=0.01$.](../img/mdl-prob-power.svg)\n",
    "\n",
    "### Test Statistics, $p$-values, and Significance\n",
    "\n",
    "To run a test we compress the data into a single **test statistic** $T(x)$, a scalar chosen so that extreme values are unlikely under $H_0$; the mean difference between two groups is a natural choice. Under $H_0$ the statistic has a known (often approximately Gaussian) *null distribution*, which is what lets us judge whether an observed value is surprising.\n",
    "\n",
    "The verdict is delivered by the **$p$-value**: the probability, *computed under $H_0$*, of seeing a statistic at least as extreme as the one we observed. For a two-sided test (the common case, where a deviation in either direction counts as evidence),\n",
    "\n",
    "$$\n",
    "p\\textrm{-value} = P_{H_0}\\bigl(|T(X)| \\ge |T(x)|\\bigr),\n",
    "$$\n",
    "\n",
    "valid when the null distribution is symmetric about $0$; in general the two-sided $p$-value is $2\\,\\min\\{P_{H_0}(T\\ge t),\\, P_{H_0}(T\\le t)\\}$, which doubles the smaller tail. The one-sided version uses a single tail. We reject $H_0$ when $p \\le \\alpha$. Geometrically, the rejection region is the set of statistic values whose $p$-value falls below $\\alpha$; the figure shows it for a Gaussian null at $\\alpha=0.05$ as the two tails beyond the critical values $\\pm 1.96$ (the $z$-test's rejection region $|T|\\ge1.96$ from above), together holding $5\\%$ of the probability. A statistic landing in those tails would be very unlikely if $H_0$ held, so we reject.\n",
    "\n",
    "![Statistical significance: under the null distribution, the central region holds probability $1-\\alpha$ and the two tails together hold $\\alpha$. A test statistic in the tails is unlikely under $H_0$, so we reject it.](../img/mdl-prob-significance.svg)\n",
    "\n",
    "A persistent warning is in order, because the $p$-value is among the most misread numbers in science [@Wasserstein.Lazar.2016]. It is $P(\\textrm{data this extreme}\\mid H_0)$, a statement about the data *given* the null; it is *not* $P(H_0\\mid\\textrm{data})$, the probability the null is true, which would require a prior and Bayes' rule. A large $p$-value does *not* confirm $H_0$; it means only that we failed to detect an effect, possibly because the test was underpowered.\n",
    "\n",
    "A subtler trap is *multiple testing*. The $\\alpha=0.05$ guarantee holds for a *single* pre-specified test; run $m$ of them under a true null (sweeping hyperparameters, comparing across benchmarks, retrying until something \"works\") and the chance of at least one spurious win is $1-(1-\\alpha)^m$, which already exceeds $0.4$ at $m=10$. Reporting only the test that cleared $p\\le\\alpha$ is *$p$-hacking*, and it is how noise gets published as a result. The simplest guard is the *Bonferroni correction* [@Bonferroni.1936]: to hold the family-wide false-positive rate at $\\alpha$, test each of the $m$ hypotheses at the stricter level $\\alpha/m$. When $m$ runs into the hundreds or thousands, as in a hyperparameter sweep or a screen of model variants, Bonferroni grows hopelessly conservative, and large-scale practice instead controls the *false discovery rate*, the expected fraction of rejections that are false, via the Benjamini--Hochberg procedure [@Benjamini.Hochberg.1995], which sorts the $m$ $p$-values in increasing order and rejects the $k$ smallest, with $k$ the largest index satisfying $p_{(k)}\\le k\\alpha/m$.\n",
    "\n",
    "To summarize, a hypothesis test proceeds in five steps:\n",
    "\n",
    "1. State $H_0$ and $H_A$.\n",
    "2. Fix the significance level $\\alpha$ and a target power $1-\\beta$ (which, with the expected effect size, sets the sample size).\n",
    "3. Collect the data.\n",
    "4. Compute the test statistic and its $p$-value under $H_0$.\n",
    "5. Reject $H_0$ if $p \\le \\alpha$; otherwise fail to reject.\n",
    "\n",
    "### A Worked Test: Comparing Two Models\n",
    "\n",
    "Let us walk through the five steps once, on the comparison practitioners face most often: is model B really better than model A, or did it just draw lucky seeds? We simulate per-seed test accuracies for the two models: twenty training runs each, with a true gap of $0.8\\%$ buried in seed-to-seed noise of comparable size. Step 1: $H_0$ is that the two models are equally good, i.e. the two accuracy samples come from the same distribution, and $H_A$ is that they differ; step 2: $\\alpha=0.05$. The test statistic is the gap between the mean accuracies. Rather than assume a Gaussian null distribution, we use a **permutation test** [@Fisher.1935], which manufactures the null distribution from the data itself: if $H_0$ holds, the labels \"A\" and \"B\" carry no information (the $40$ numbers are *exchangeable*), so shuffling the labels and recomputing the gap, many times over, shows exactly how large a gap arises by pure chance. The two-sided $p$-value is the fraction of shuffles producing a gap at least as extreme as the observed one (counting the observed labeling itself among them, which keeps the estimate valid and never exactly zero)."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "statistics-permutation-test"
   },
   "source": [
    "rng = onp.random.default_rng(1)\n",
    "num_seeds = 20                               # 20 training runs per model\n",
    "acc_a = rng.normal(0.850, 0.010, num_seeds)  # per-seed accuracy, model A\n",
    "acc_b = rng.normal(0.858, 0.010, num_seeds)  # model B: a real +0.008 gap\n",
    "observed = acc_b.mean() - acc_a.mean()       # test statistic: gap in means\n",
    "\n",
    "pooled = onp.concatenate([acc_a, acc_b])     # under H_0 the labels are arbitrary\n",
    "B = 10000\n",
    "gaps = onp.empty(B)\n",
    "for b in range(B):\n",
    "    perm = rng.permutation(pooled)           # shuffle the model labels\n",
    "    gaps[b] = perm[num_seeds:].mean() - perm[:num_seeds].mean()\n",
    "p_value = (1 + (onp.abs(gaps) >= abs(observed)).sum()) / (B + 1)  # two-sided\n",
    "print(f'observed gap        = {observed:.4f}')\n",
    "print(f'permutation p-value = {p_value:.4f}')\n",
    "\n",
    "d2l.set_figsize()\n",
    "d2l.plt.hist(gaps, bins=50, color='#1f77b4', alpha=0.8)\n",
    "d2l.plt.axvline(observed, color='#ff7f0e', lw=2)      # the observed gap\n",
    "d2l.plt.axvline(-observed, color='#ff7f0e', lw=2, ls='--')\n",
    "d2l.plt.xlabel('gap under label shuffling')\n",
    "d2l.plt.ylabel('count');"
   ],
   "execution_count": 6,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "observed gap        = 0.0073\npermutation p-value = 0.0197\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "image/svg+xml": "<?xml version=\"1.0\" encoding=\"utf-8\" standalone=\"no\"?>\n<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n  \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n<svg xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"249.465625pt\" height=\"183.35625pt\" viewBox=\"0 0 249.465625 183.35625\" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">\n <metadata>\n  <rdf:RDF xmlns:dc=\"http://purl.org/dc/elements/1.1/\" xmlns:cc=\"http://creativecommons.org/ns#\" xmlns:rdf=\"http://www.w3.org/1999/02/22-rdf-syntax-ns#\">\n   <cc:Work>\n    <dc:type rdf:resource=\"http://purl.org/dc/dcmitype/StillImage\"/>\n    <dc:date>2026-07-19T13:42:45.831653</dc:date>\n    <dc:format>image/svg+xml</dc:format>\n    <dc:creator>\n     <cc:Agent>\n      <dc:title>Matplotlib v3.10.8, https://matplotlib.org/</dc:title>\n     </cc:Agent>\n    </dc:creator>\n   </cc:Work>\n  </rdf:RDF>\n </metadata>\n <defs>\n  <style type=\"text/css\">*{stroke-linejoin: round; stroke-linecap: butt}</style>\n </defs>\n <g id=\"figure_1\">\n  <g id=\"patch_1\">\n   <path d=\"M 0 183.35625 \nL 249.465625 183.35625 \nL 249.465625 0 \nL 0 0 \nz\n\" style=\"fill: #ffffff\"/>\n  </g>\n  <g id=\"axes_1\">\n   <g id=\"patch_2\">\n    <path d=\"M 46.965625 145.8 \nL 242.265625 145.8 \nL 242.265625 7.2 \nL 46.965625 7.2 \nz\n\" style=\"fill: #ffffff\"/>\n   </g>\n   <g id=\"patch_3\">\n    <path d=\"M 55.842898 145.8 \nL 59.393807 145.8 \nL 59.393807 145.563016 \nL 55.842898 145.563016 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_4\">\n    <path d=\"M 59.393807 145.8 \nL 62.944716 145.8 \nL 62.944716 145.563016 \nL 59.393807 145.563016 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_5\">\n    <path d=\"M 62.944716 145.8 \nL 66.495625 145.8 \nL 66.495625 144.852065 \nL 62.944716 144.852065 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_6\">\n    <path d=\"M 66.495625 145.8 \nL 70.046534 145.8 \nL 70.046534 144.141113 \nL 66.495625 144.141113 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_7\">\n    <path d=\"M 70.046534 145.8 \nL 73.597443 145.8 \nL 73.597443 143.193178 \nL 70.046534 143.193178 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_8\">\n    <path d=\"M 73.597443 145.8 \nL 77.148352 145.8 \nL 77.148352 140.823339 \nL 73.597443 140.823339 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_9\">\n    <path d=\"M 77.148352 145.8 \nL 80.699261 145.8 \nL 80.699261 140.112388 \nL 77.148352 140.112388 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_10\">\n    <path d=\"M 80.699261 145.8 \nL 84.25017 145.8 \nL 84.25017 138.927469 \nL 80.699261 138.927469 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_11\">\n    <path d=\"M 84.25017 145.8 \nL 87.80108 145.8 \nL 87.80108 138.216517 \nL 84.25017 138.216517 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_12\">\n    <path d=\"M 87.80108 145.8 \nL 91.351989 145.8 \nL 91.351989 131.580969 \nL 87.80108 131.580969 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_13\">\n    <path d=\"M 91.351989 145.8 \nL 94.902898 145.8 \nL 94.902898 121.627648 \nL 91.351989 121.627648 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_14\">\n    <path d=\"M 94.902898 145.8 \nL 98.453807 145.8 \nL 98.453807 119.020826 \nL 94.902898 119.020826 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_15\">\n    <path d=\"M 98.453807 145.8 \nL 102.004716 145.8 \nL 102.004716 112.622262 \nL 98.453807 112.622262 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_16\">\n    <path d=\"M 102.004716 145.8 \nL 105.555625 145.8 \nL 105.555625 101.247038 \nL 102.004716 101.247038 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_17\">\n    <path d=\"M 105.555625 145.8 \nL 109.106534 145.8 \nL 109.106534 95.796409 \nL 105.555625 95.796409 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_18\">\n    <path d=\"M 109.106534 145.8 \nL 112.657443 145.8 \nL 112.657443 76.363734 \nL 109.106534 76.363734 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_19\">\n    <path d=\"M 112.657443 145.8 \nL 116.208352 145.8 \nL 116.208352 79.444524 \nL 112.657443 79.444524 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_20\">\n    <path d=\"M 116.208352 145.8 \nL 119.759261 145.8 \nL 119.759261 66.647397 \nL 116.208352 66.647397 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_21\">\n    <path d=\"M 119.759261 145.8 \nL 123.31017 145.8 \nL 123.31017 46.977738 \nL 119.759261 46.977738 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_22\">\n    <path d=\"M 123.31017 145.8 \nL 126.86108 145.8 \nL 126.86108 40.105206 \nL 123.31017 40.105206 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_23\">\n    <path d=\"M 126.86108 145.8 \nL 130.411989 145.8 \nL 130.411989 33.943627 \nL 126.86108 33.943627 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_24\">\n    <path d=\"M 130.411989 145.8 \nL 133.962898 145.8 \nL 133.962898 28.966966 \nL 130.411989 28.966966 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_25\">\n    <path d=\"M 133.962898 145.8 \nL 137.513807 145.8 \nL 137.513807 13.8 \nL 133.962898 13.8 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_26\">\n    <path d=\"M 137.513807 145.8 \nL 141.064716 145.8 \nL 141.064716 23.990305 \nL 137.513807 23.990305 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_27\">\n    <path d=\"M 141.064716 145.8 \nL 144.615625 145.8 \nL 144.615625 18.776661 \nL 141.064716 18.776661 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_28\">\n    <path d=\"M 144.615625 145.8 \nL 148.166534 145.8 \nL 148.166534 23.279354 \nL 144.615625 23.279354 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_29\">\n    <path d=\"M 148.166534 145.8 \nL 151.717443 145.8 \nL 151.717443 27.782047 \nL 148.166534 27.782047 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_30\">\n    <path d=\"M 151.717443 145.8 \nL 155.268352 145.8 \nL 155.268352 33.706643 \nL 151.717443 33.706643 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_31\">\n    <path d=\"M 155.268352 145.8 \nL 158.819261 145.8 \nL 158.819261 35.128546 \nL 155.268352 35.128546 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_32\">\n    <path d=\"M 158.819261 145.8 \nL 162.37017 145.8 \nL 162.37017 44.133932 \nL 158.819261 44.133932 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_33\">\n    <path d=\"M 162.37017 145.8 \nL 165.92108 145.8 \nL 165.92108 52.191382 \nL 162.37017 52.191382 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_34\">\n    <path d=\"M 165.92108 145.8 \nL 169.471989 145.8 \nL 169.471989 60.248833 \nL 165.92108 60.248833 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_35\">\n    <path d=\"M 169.471989 145.8 \nL 173.022898 145.8 \nL 173.022898 72.808977 \nL 169.471989 72.808977 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_36\">\n    <path d=\"M 173.022898 145.8 \nL 176.573807 145.8 \nL 176.573807 87.738959 \nL 173.022898 87.738959 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_37\">\n    <path d=\"M 176.573807 145.8 \nL 180.124716 145.8 \nL 180.124716 94.374506 \nL 176.573807 94.374506 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_38\">\n    <path d=\"M 180.124716 145.8 \nL 183.675625 145.8 \nL 183.675625 101.484022 \nL 180.124716 101.484022 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_39\">\n    <path d=\"M 183.675625 145.8 \nL 187.226534 145.8 \nL 187.226534 105.749731 \nL 183.675625 105.749731 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_40\">\n    <path d=\"M 187.226534 145.8 \nL 190.777443 145.8 \nL 190.777443 118.07289 \nL 187.226534 118.07289 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_41\">\n    <path d=\"M 190.777443 145.8 \nL 194.328352 145.8 \nL 194.328352 120.205745 \nL 190.777443 120.205745 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_42\">\n    <path d=\"M 194.328352 145.8 \nL 197.879261 145.8 \nL 197.879261 129.211131 \nL 194.328352 129.211131 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_43\">\n    <path d=\"M 197.879261 145.8 \nL 201.43017 145.8 \nL 201.43017 136.794614 \nL 197.879261 136.794614 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_44\">\n    <path d=\"M 201.43017 145.8 \nL 204.98108 145.8 \nL 204.98108 137.742549 \nL 201.43017 137.742549 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_45\">\n    <path d=\"M 204.98108 145.8 \nL 208.531989 145.8 \nL 208.531989 140.112388 \nL 204.98108 140.112388 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_46\">\n    <path d=\"M 208.531989 145.8 \nL 212.082898 145.8 \nL 212.082898 140.112388 \nL 208.531989 140.112388 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_47\">\n    <path d=\"M 212.082898 145.8 \nL 215.633807 145.8 \nL 215.633807 142.956194 \nL 212.082898 142.956194 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_48\">\n    <path d=\"M 215.633807 145.8 \nL 219.184716 145.8 \nL 219.184716 144.378097 \nL 215.633807 144.378097 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_49\">\n    <path d=\"M 219.184716 145.8 \nL 222.735625 145.8 \nL 222.735625 144.378097 \nL 219.184716 144.378097 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_50\">\n    <path d=\"M 222.735625 145.8 \nL 226.286534 145.8 \nL 226.286534 144.852065 \nL 222.735625 144.852065 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_51\">\n    <path d=\"M 226.286534 145.8 \nL 229.837443 145.8 \nL 229.837443 145.089048 \nL 226.286534 145.089048 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"patch_52\">\n    <path d=\"M 229.837443 145.8 \nL 233.388352 145.8 \nL 233.388352 145.326032 \nL 229.837443 145.326032 \nz\n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: #1f77b4; opacity: 0.8\"/>\n   </g>\n   <g id=\"matplotlib.axis_1\">\n    <g id=\"xtick_1\">\n     <g id=\"line2d_1\">\n      <defs>\n       <path id=\"mb143b2d7c8\" d=\"M 0 0 \nL 0 3.5 \n\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#mb143b2d7c8\" x=\"62.147897\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_1\">\n      <!-- −0.010 -->\n      <g transform=\"translate(43.64399 160.398438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-2212\" d=\"M 678 2272 \nL 4684 2272 \nL 4684 1741 \nL 678 1741 \nL 678 2272 \nz\n\" transform=\"scale(0.015625)\"/>\n        <path id=\"DejaVuSans-30\" d=\"M 2034 4250 \nQ 1547 4250 1301 3770 \nQ 1056 3291 1056 2328 \nQ 1056 1369 1301 889 \nQ 1547 409 2034 409 \nQ 2525 409 2770 889 \nQ 3016 1369 3016 2328 \nQ 3016 3291 2770 3770 \nQ 2525 4250 2034 4250 \nz\nM 2034 4750 \nQ 2819 4750 3233 4129 \nQ 3647 3509 3647 2328 \nQ 3647 1150 3233 529 \nQ 2819 -91 2034 -91 \nQ 1250 -91 836 529 \nQ 422 1150 422 2328 \nQ 422 3509 836 4129 \nQ 1250 4750 2034 4750 \nz\n\" transform=\"scale(0.015625)\"/>\n        <path id=\"DejaVuSans-2e\" d=\"M 684 794 \nL 1344 794 \nL 1344 0 \nL 684 0 \nL 684 794 \nz\n\" transform=\"scale(0.015625)\"/>\n        <path id=\"DejaVuSans-31\" d=\"M 794 531 \nL 1825 531 \nL 1825 4091 \nL 703 3866 \nL 703 4441 \nL 1819 4666 \nL 2450 4666 \nL 2450 531 \nL 3481 531 \nL 3481 0 \nL 794 0 \nL 794 531 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-2212\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(83.789062 0)\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(147.412109 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(179.199219 0)\"/>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(242.822266 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(306.445312 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_2\">\n     <g id=\"line2d_2\">\n      <g>\n       <use xlink:href=\"#mb143b2d7c8\" x=\"102.787539\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_2\">\n      <!-- −0.005 -->\n      <g transform=\"translate(84.283633 160.398438) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-35\" d=\"M 691 4666 \nL 3169 4666 \nL 3169 4134 \nL 1269 4134 \nL 1269 2991 \nQ 1406 3038 1543 3061 \nQ 1681 3084 1819 3084 \nQ 2600 3084 3056 2656 \nQ 3513 2228 3513 1497 \nQ 3513 744 3044 326 \nQ 2575 -91 1722 -91 \nQ 1428 -91 1123 -41 \nQ 819 9 494 109 \nL 494 744 \nQ 775 591 1075 516 \nQ 1375 441 1709 441 \nQ 2250 441 2565 725 \nQ 2881 1009 2881 1497 \nQ 2881 1984 2565 2268 \nQ 2250 2553 1709 2553 \nQ 1456 2553 1204 2497 \nQ 953 2441 691 2322 \nL 691 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-2212\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(83.789062 0)\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(147.412109 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(179.199219 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(242.822266 0)\"/>\n       <use xlink:href=\"#DejaVuSans-35\" transform=\"translate(306.445312 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_3\">\n     <g id=\"line2d_3\">\n      <g>\n       <use xlink:href=\"#mb143b2d7c8\" x=\"143.427182\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_3\">\n      <!-- 0.000 -->\n      <g transform=\"translate(129.113119 160.398438) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(95.410156 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(159.033203 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(222.65625 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_4\">\n     <g id=\"line2d_4\">\n      <g>\n       <use xlink:href=\"#mb143b2d7c8\" x=\"184.066825\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_4\">\n      <!-- 0.005 -->\n      <g transform=\"translate(169.752762 160.398438) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(95.410156 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(159.033203 0)\"/>\n       <use xlink:href=\"#DejaVuSans-35\" transform=\"translate(222.65625 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"xtick_5\">\n     <g id=\"line2d_5\">\n      <g>\n       <use xlink:href=\"#mb143b2d7c8\" x=\"224.706467\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_5\">\n      <!-- 0.010 -->\n      <g transform=\"translate(210.392405 160.398438) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n       <use xlink:href=\"#DejaVuSans-2e\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(95.410156 0)\"/>\n       <use xlink:href=\"#DejaVuSans-31\" transform=\"translate(159.033203 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(222.65625 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_6\">\n     <!-- gap under label shuffling -->\n     <g transform=\"translate(82.115625 174.076563) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-67\" d=\"M 2906 1791 \nQ 2906 2416 2648 2759 \nQ 2391 3103 1925 3103 \nQ 1463 3103 1205 2759 \nQ 947 2416 947 1791 \nQ 947 1169 1205 825 \nQ 1463 481 1925 481 \nQ 2391 481 2648 825 \nQ 2906 1169 2906 1791 \nz\nM 3481 434 \nQ 3481 -459 3084 -895 \nQ 2688 -1331 1869 -1331 \nQ 1566 -1331 1297 -1286 \nQ 1028 -1241 775 -1147 \nL 775 -588 \nQ 1028 -725 1275 -790 \nQ 1522 -856 1778 -856 \nQ 2344 -856 2625 -561 \nQ 2906 -266 2906 331 \nL 2906 616 \nQ 2728 306 2450 153 \nQ 2172 0 1784 0 \nQ 1141 0 747 490 \nQ 353 981 353 1791 \nQ 353 2603 747 3093 \nQ 1141 3584 1784 3584 \nQ 2172 3584 2450 3431 \nQ 2728 3278 2906 2969 \nL 2906 3500 \nL 3481 3500 \nL 3481 434 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-61\" d=\"M 2194 1759 \nQ 1497 1759 1228 1600 \nQ 959 1441 959 1056 \nQ 959 750 1161 570 \nQ 1363 391 1709 391 \nQ 2188 391 2477 730 \nQ 2766 1069 2766 1631 \nL 2766 1759 \nL 2194 1759 \nz\nM 3341 1997 \nL 3341 0 \nL 2766 0 \nL 2766 531 \nQ 2569 213 2275 61 \nQ 1981 -91 1556 -91 \nQ 1019 -91 701 211 \nQ 384 513 384 1019 \nQ 384 1609 779 1909 \nQ 1175 2209 1959 2209 \nL 2766 2209 \nL 2766 2266 \nQ 2766 2663 2505 2880 \nQ 2244 3097 1772 3097 \nQ 1472 3097 1187 3025 \nQ 903 2953 641 2809 \nL 641 3341 \nQ 956 3463 1253 3523 \nQ 1550 3584 1831 3584 \nQ 2591 3584 2966 3190 \nQ 3341 2797 3341 1997 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-70\" d=\"M 1159 525 \nL 1159 -1331 \nL 581 -1331 \nL 581 3500 \nL 1159 3500 \nL 1159 2969 \nQ 1341 3281 1617 3432 \nQ 1894 3584 2278 3584 \nQ 2916 3584 3314 3078 \nQ 3713 2572 3713 1747 \nQ 3713 922 3314 415 \nQ 2916 -91 2278 -91 \nQ 1894 -91 1617 61 \nQ 1341 213 1159 525 \nz\nM 3116 1747 \nQ 3116 2381 2855 2742 \nQ 2594 3103 2138 3103 \nQ 1681 3103 1420 2742 \nQ 1159 2381 1159 1747 \nQ 1159 1113 1420 752 \nQ 1681 391 2138 391 \nQ 2594 391 2855 752 \nQ 3116 1113 3116 1747 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-20\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-75\" d=\"M 544 1381 \nL 544 3500 \nL 1119 3500 \nL 1119 1403 \nQ 1119 906 1312 657 \nQ 1506 409 1894 409 \nQ 2359 409 2629 706 \nQ 2900 1003 2900 1516 \nL 2900 3500 \nL 3475 3500 \nL 3475 0 \nL 2900 0 \nL 2900 538 \nQ 2691 219 2414 64 \nQ 2138 -91 1772 -91 \nQ 1169 -91 856 284 \nQ 544 659 544 1381 \nz\nM 1991 3584 \nL 1991 3584 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6e\" d=\"M 3513 2113 \nL 3513 0 \nL 2938 0 \nL 2938 2094 \nQ 2938 2591 2744 2837 \nQ 2550 3084 2163 3084 \nQ 1697 3084 1428 2787 \nQ 1159 2491 1159 1978 \nL 1159 0 \nL 581 0 \nL 581 3500 \nL 1159 3500 \nL 1159 2956 \nQ 1366 3272 1645 3428 \nQ 1925 3584 2291 3584 \nQ 2894 3584 3203 3211 \nQ 3513 2838 3513 2113 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-64\" d=\"M 2906 2969 \nL 2906 4863 \nL 3481 4863 \nL 3481 0 \nL 2906 0 \nL 2906 525 \nQ 2725 213 2448 61 \nQ 2172 -91 1784 -91 \nQ 1150 -91 751 415 \nQ 353 922 353 1747 \nQ 353 2572 751 3078 \nQ 1150 3584 1784 3584 \nQ 2172 3584 2448 3432 \nQ 2725 3281 2906 2969 \nz\nM 947 1747 \nQ 947 1113 1208 752 \nQ 1469 391 1925 391 \nQ 2381 391 2643 752 \nQ 2906 1113 2906 1747 \nQ 2906 2381 2643 2742 \nQ 2381 3103 1925 3103 \nQ 1469 3103 1208 2742 \nQ 947 2381 947 1747 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-65\" d=\"M 3597 1894 \nL 3597 1613 \nL 953 1613 \nQ 991 1019 1311 708 \nQ 1631 397 2203 397 \nQ 2534 397 2845 478 \nQ 3156 559 3463 722 \nL 3463 178 \nQ 3153 47 2828 -22 \nQ 2503 -91 2169 -91 \nQ 1331 -91 842 396 \nQ 353 884 353 1716 \nQ 353 2575 817 3079 \nQ 1281 3584 2069 3584 \nQ 2775 3584 3186 3129 \nQ 3597 2675 3597 1894 \nz\nM 3022 2063 \nQ 3016 2534 2758 2815 \nQ 2500 3097 2075 3097 \nQ 1594 3097 1305 2825 \nQ 1016 2553 972 2059 \nL 3022 2063 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-72\" d=\"M 2631 2963 \nQ 2534 3019 2420 3045 \nQ 2306 3072 2169 3072 \nQ 1681 3072 1420 2755 \nQ 1159 2438 1159 1844 \nL 1159 0 \nL 581 0 \nL 581 3500 \nL 1159 3500 \nL 1159 2956 \nQ 1341 3275 1631 3429 \nQ 1922 3584 2338 3584 \nQ 2397 3584 2469 3576 \nQ 2541 3569 2628 3553 \nL 2631 2963 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6c\" d=\"M 603 4863 \nL 1178 4863 \nL 1178 0 \nL 603 0 \nL 603 4863 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-62\" d=\"M 3116 1747 \nQ 3116 2381 2855 2742 \nQ 2594 3103 2138 3103 \nQ 1681 3103 1420 2742 \nQ 1159 2381 1159 1747 \nQ 1159 1113 1420 752 \nQ 1681 391 2138 391 \nQ 2594 391 2855 752 \nQ 3116 1113 3116 1747 \nz\nM 1159 2969 \nQ 1341 3281 1617 3432 \nQ 1894 3584 2278 3584 \nQ 2916 3584 3314 3078 \nQ 3713 2572 3713 1747 \nQ 3713 922 3314 415 \nQ 2916 -91 2278 -91 \nQ 1894 -91 1617 61 \nQ 1341 213 1159 525 \nL 1159 0 \nL 581 0 \nL 581 4863 \nL 1159 4863 \nL 1159 2969 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-73\" d=\"M 2834 3397 \nL 2834 2853 \nQ 2591 2978 2328 3040 \nQ 2066 3103 1784 3103 \nQ 1356 3103 1142 2972 \nQ 928 2841 928 2578 \nQ 928 2378 1081 2264 \nQ 1234 2150 1697 2047 \nL 1894 2003 \nQ 2506 1872 2764 1633 \nQ 3022 1394 3022 966 \nQ 3022 478 2636 193 \nQ 2250 -91 1575 -91 \nQ 1294 -91 989 -36 \nQ 684 19 347 128 \nL 347 722 \nQ 666 556 975 473 \nQ 1284 391 1588 391 \nQ 1994 391 2212 530 \nQ 2431 669 2431 922 \nQ 2431 1156 2273 1281 \nQ 2116 1406 1581 1522 \nL 1381 1569 \nQ 847 1681 609 1914 \nQ 372 2147 372 2553 \nQ 372 3047 722 3315 \nQ 1072 3584 1716 3584 \nQ 2034 3584 2315 3537 \nQ 2597 3491 2834 3397 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-68\" d=\"M 3513 2113 \nL 3513 0 \nL 2938 0 \nL 2938 2094 \nQ 2938 2591 2744 2837 \nQ 2550 3084 2163 3084 \nQ 1697 3084 1428 2787 \nQ 1159 2491 1159 1978 \nL 1159 0 \nL 581 0 \nL 581 4863 \nL 1159 4863 \nL 1159 2956 \nQ 1366 3272 1645 3428 \nQ 1925 3584 2291 3584 \nQ 2894 3584 3203 3211 \nQ 3513 2838 3513 2113 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-66\" d=\"M 2375 4863 \nL 2375 4384 \nL 1825 4384 \nQ 1516 4384 1395 4259 \nQ 1275 4134 1275 3809 \nL 1275 3500 \nL 2222 3500 \nL 2222 3053 \nL 1275 3053 \nL 1275 0 \nL 697 0 \nL 697 3053 \nL 147 3053 \nL 147 3500 \nL 697 3500 \nL 697 3744 \nQ 697 4328 969 4595 \nQ 1241 4863 1831 4863 \nL 2375 4863 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-69\" d=\"M 603 3500 \nL 1178 3500 \nL 1178 0 \nL 603 0 \nL 603 3500 \nz\nM 603 4863 \nL 1178 4863 \nL 1178 4134 \nL 603 4134 \nL 603 4863 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-67\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(63.476562 0)\"/>\n      <use xlink:href=\"#DejaVuSans-70\" transform=\"translate(124.755859 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(188.232422 0)\"/>\n      <use xlink:href=\"#DejaVuSans-75\" transform=\"translate(220.019531 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(283.398438 0)\"/>\n      <use xlink:href=\"#DejaVuSans-64\" transform=\"translate(346.777344 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(410.253906 0)\"/>\n      <use xlink:href=\"#DejaVuSans-72\" transform=\"translate(471.777344 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(512.890625 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6c\" transform=\"translate(544.677734 0)\"/>\n      <use xlink:href=\"#DejaVuSans-61\" transform=\"translate(572.460938 0)\"/>\n      <use xlink:href=\"#DejaVuSans-62\" transform=\"translate(633.740234 0)\"/>\n      <use xlink:href=\"#DejaVuSans-65\" transform=\"translate(697.216797 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6c\" transform=\"translate(758.740234 0)\"/>\n      <use xlink:href=\"#DejaVuSans-20\" transform=\"translate(786.523438 0)\"/>\n      <use xlink:href=\"#DejaVuSans-73\" transform=\"translate(818.310547 0)\"/>\n      <use xlink:href=\"#DejaVuSans-68\" transform=\"translate(870.410156 0)\"/>\n      <use xlink:href=\"#DejaVuSans-75\" transform=\"translate(933.789062 0)\"/>\n      <use xlink:href=\"#DejaVuSans-66\" transform=\"translate(997.167969 0)\"/>\n      <use xlink:href=\"#DejaVuSans-66\" transform=\"translate(1032.373047 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6c\" transform=\"translate(1067.578125 0)\"/>\n      <use xlink:href=\"#DejaVuSans-69\" transform=\"translate(1095.361328 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(1123.144531 0)\"/>\n      <use xlink:href=\"#DejaVuSans-67\" transform=\"translate(1186.523438 0)\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"matplotlib.axis_2\">\n    <g id=\"ytick_1\">\n     <g id=\"line2d_6\">\n      <defs>\n       <path id=\"mb3a202c2f8\" d=\"M 0 0 \nL -3.5 0 \n\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </defs>\n      <g>\n       <use xlink:href=\"#mb3a202c2f8\" x=\"46.965625\" y=\"145.8\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_7\">\n      <!-- 0 -->\n      <g transform=\"translate(33.603125 149.599219) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-30\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_2\">\n     <g id=\"line2d_7\">\n      <g>\n       <use xlink:href=\"#mb3a202c2f8\" x=\"46.965625\" y=\"122.101616\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_8\">\n      <!-- 100 -->\n      <g transform=\"translate(20.878125 125.900835) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-31\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(127.246094 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_3\">\n     <g id=\"line2d_8\">\n      <g>\n       <use xlink:href=\"#mb3a202c2f8\" x=\"46.965625\" y=\"98.403232\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_9\">\n      <!-- 200 -->\n      <g transform=\"translate(20.878125 102.20245) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-32\" d=\"M 1228 531 \nL 3431 531 \nL 3431 0 \nL 469 0 \nL 469 531 \nQ 828 903 1448 1529 \nQ 2069 2156 2228 2338 \nQ 2531 2678 2651 2914 \nQ 2772 3150 2772 3378 \nQ 2772 3750 2511 3984 \nQ 2250 4219 1831 4219 \nQ 1534 4219 1204 4116 \nQ 875 4013 500 3803 \nL 500 4441 \nQ 881 4594 1212 4672 \nQ 1544 4750 1819 4750 \nQ 2544 4750 2975 4387 \nQ 3406 4025 3406 3419 \nQ 3406 3131 3298 2873 \nQ 3191 2616 2906 2266 \nQ 2828 2175 2409 1742 \nQ 1991 1309 1228 531 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-32\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(127.246094 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_4\">\n     <g id=\"line2d_9\">\n      <g>\n       <use xlink:href=\"#mb3a202c2f8\" x=\"46.965625\" y=\"74.704847\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_10\">\n      <!-- 300 -->\n      <g transform=\"translate(20.878125 78.504066) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-33\" d=\"M 2597 2516 \nQ 3050 2419 3304 2112 \nQ 3559 1806 3559 1356 \nQ 3559 666 3084 287 \nQ 2609 -91 1734 -91 \nQ 1441 -91 1130 -33 \nQ 819 25 488 141 \nL 488 750 \nQ 750 597 1062 519 \nQ 1375 441 1716 441 \nQ 2309 441 2620 675 \nQ 2931 909 2931 1356 \nQ 2931 1769 2642 2001 \nQ 2353 2234 1838 2234 \nL 1294 2234 \nL 1294 2753 \nL 1863 2753 \nQ 2328 2753 2575 2939 \nQ 2822 3125 2822 3475 \nQ 2822 3834 2567 4026 \nQ 2313 4219 1838 4219 \nQ 1578 4219 1281 4162 \nQ 984 4106 628 3988 \nL 628 4550 \nQ 988 4650 1302 4700 \nQ 1616 4750 1894 4750 \nQ 2613 4750 3031 4423 \nQ 3450 4097 3450 3541 \nQ 3450 3153 3228 2886 \nQ 3006 2619 2597 2516 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-33\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(127.246094 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_5\">\n     <g id=\"line2d_10\">\n      <g>\n       <use xlink:href=\"#mb3a202c2f8\" x=\"46.965625\" y=\"51.006463\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_11\">\n      <!-- 400 -->\n      <g transform=\"translate(20.878125 54.805682) scale(0.1 -0.1)\">\n       <defs>\n        <path id=\"DejaVuSans-34\" d=\"M 2419 4116 \nL 825 1625 \nL 2419 1625 \nL 2419 4116 \nz\nM 2253 4666 \nL 3047 4666 \nL 3047 1625 \nL 3713 1625 \nL 3713 1100 \nL 3047 1100 \nL 3047 0 \nL 2419 0 \nL 2419 1100 \nL 313 1100 \nL 313 1709 \nL 2253 4666 \nz\n\" transform=\"scale(0.015625)\"/>\n       </defs>\n       <use xlink:href=\"#DejaVuSans-34\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(127.246094 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"ytick_6\">\n     <g id=\"line2d_11\">\n      <g>\n       <use xlink:href=\"#mb3a202c2f8\" x=\"46.965625\" y=\"27.308079\" style=\"stroke: #000000; stroke-width: 0.8\"/>\n      </g>\n     </g>\n     <g id=\"text_12\">\n      <!-- 500 -->\n      <g transform=\"translate(20.878125 31.107298) scale(0.1 -0.1)\">\n       <use xlink:href=\"#DejaVuSans-35\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(63.623047 0)\"/>\n       <use xlink:href=\"#DejaVuSans-30\" transform=\"translate(127.246094 0)\"/>\n      </g>\n     </g>\n    </g>\n    <g id=\"text_13\">\n     <!-- count -->\n     <g transform=\"translate(14.798438 90.60625) rotate(-90) scale(0.1 -0.1)\">\n      <defs>\n       <path id=\"DejaVuSans-63\" d=\"M 3122 3366 \nL 3122 2828 \nQ 2878 2963 2633 3030 \nQ 2388 3097 2138 3097 \nQ 1578 3097 1268 2742 \nQ 959 2388 959 1747 \nQ 959 1106 1268 751 \nQ 1578 397 2138 397 \nQ 2388 397 2633 464 \nQ 2878 531 3122 666 \nL 3122 134 \nQ 2881 22 2623 -34 \nQ 2366 -91 2075 -91 \nQ 1284 -91 818 406 \nQ 353 903 353 1747 \nQ 353 2603 823 3093 \nQ 1294 3584 2113 3584 \nQ 2378 3584 2631 3529 \nQ 2884 3475 3122 3366 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-6f\" d=\"M 1959 3097 \nQ 1497 3097 1228 2736 \nQ 959 2375 959 1747 \nQ 959 1119 1226 758 \nQ 1494 397 1959 397 \nQ 2419 397 2687 759 \nQ 2956 1122 2956 1747 \nQ 2956 2369 2687 2733 \nQ 2419 3097 1959 3097 \nz\nM 1959 3584 \nQ 2709 3584 3137 3096 \nQ 3566 2609 3566 1747 \nQ 3566 888 3137 398 \nQ 2709 -91 1959 -91 \nQ 1206 -91 779 398 \nQ 353 888 353 1747 \nQ 353 2609 779 3096 \nQ 1206 3584 1959 3584 \nz\n\" transform=\"scale(0.015625)\"/>\n       <path id=\"DejaVuSans-74\" d=\"M 1172 4494 \nL 1172 3500 \nL 2356 3500 \nL 2356 3053 \nL 1172 3053 \nL 1172 1153 \nQ 1172 725 1289 603 \nQ 1406 481 1766 481 \nL 2356 481 \nL 2356 0 \nL 1766 0 \nQ 1100 0 847 248 \nQ 594 497 594 1153 \nL 594 3053 \nL 172 3053 \nL 172 3500 \nL 594 3500 \nL 594 4494 \nL 1172 4494 \nz\n\" transform=\"scale(0.015625)\"/>\n      </defs>\n      <use xlink:href=\"#DejaVuSans-63\"/>\n      <use xlink:href=\"#DejaVuSans-6f\" transform=\"translate(54.980469 0)\"/>\n      <use xlink:href=\"#DejaVuSans-75\" transform=\"translate(116.162109 0)\"/>\n      <use xlink:href=\"#DejaVuSans-6e\" transform=\"translate(179.541016 0)\"/>\n      <use xlink:href=\"#DejaVuSans-74\" transform=\"translate(242.919922 0)\"/>\n     </g>\n    </g>\n   </g>\n   <g id=\"line2d_12\">\n    <path d=\"M 202.812092 145.8 \nL 202.812092 7.2 \n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: none; stroke: #ff7f0e; stroke-width: 2; stroke-linecap: square\"/>\n   </g>\n   <g id=\"line2d_13\">\n    <path d=\"M 84.042271 145.8 \nL 84.042271 7.2 \n\" clip-path=\"url(#p835a1d5bfa)\" style=\"fill: none; stroke-dasharray: 7.4,3.2; stroke-dashoffset: 0; stroke: #ff7f0e; stroke-width: 2\"/>\n   </g>\n   <g id=\"patch_53\">\n    <path d=\"M 46.965625 145.8 \nL 46.965625 7.2 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_54\">\n    <path d=\"M 242.265625 145.8 \nL 242.265625 7.2 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_55\">\n    <path d=\"M 46.965625 145.8 \nL 242.265625 145.8 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n   <g id=\"patch_56\">\n    <path d=\"M 46.965625 7.2 \nL 242.265625 7.2 \n\" style=\"fill: none; stroke: #000000; stroke-width: 0.8; stroke-linejoin: miter; stroke-linecap: square\"/>\n   </g>\n  </g>\n </g>\n <defs>\n  <clipPath id=\"p835a1d5bfa\">\n   <rect x=\"46.965625\" y=\"7.2\" width=\"195.3\" height=\"138.6\"/>\n  </clipPath>\n </defs>\n</svg>\n"
     },
     "metadata": {}
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-7 fw=all -->\n",
    "\n",
    "The observed gap is $0.0073$ (model B looks better by about three quarters of an accuracy point), and only about $2\\%$ of label shuffles produce a gap that large, so $p \\approx 0.02 \\le \\alpha = 0.05$ and we reject $H_0$: the improvement is unlikely to be a fluke. The histogram of shuffled gaps makes the verdict visual: they pile up in a bell around zero (no Gaussian assumption put it there; the bell is the *finite-population* central limit theorem for means drawn without replacement from the pooled forty numbers, a relative of the i.i.d. CLT of that section), while the observed gap (solid line) sits far out in the right tail, with only a sliver of the null mass beyond the two lines that the two-sided test counts. Note how close the call is, though. A *real* $0.8\\%$ improvement, measured over twenty seeds, only just clears the bar, the power discussion above in action; with five seeds (a common budget) the same gap would usually go undetected. The permutation test assumes no Gaussian shape and works for any statistic we care to compute on the two groups; the same resample-and-recompute idea returns in the bootstrap below.\n",
    "\n",
    "One refinement matters in practice. When the comparison is *paired*, with models A and B trained and evaluated on the *same* seeds, the exchangeable objects are not the $40$ pooled numbers but the per-seed *differences*, and the right test flips the sign of each difference at random: the paired, or sign-flip, permutation test. It respects the pairing and typically gains power.\n",
    "\n",
    "## Confidence Intervals\n",
    "\n",
    "A point estimate $\\hat\\theta$ carries no notion of uncertainty: it is a single number that hides how much it would wobble on fresh data. A **confidence interval** repairs this by reporting an *interval* engineered to contain the true $\\theta$ with high probability. The idea is due to Jerzy Neyman [@Neyman.1937].\n",
    "\n",
    "### Definition and Interpretation\n",
    "\n",
    "A confidence interval for $\\theta$ is an interval $C_n$ computed from the data such that\n",
    "\n",
    "$$\n",
    "P_\\theta(C_n \\ni \\theta) \\ge 1 - \\alpha \\quad \\textrm{for all } \\theta,\n",
    "$$\n",
    "\n",
    "where $1-\\alpha$ is the *confidence level* or *coverage*. We write $C_n \\ni \\theta$ rather than $\\theta \\in C_n$ to stress where the randomness lives: $\\theta$ is a *fixed* unknown, and it is the *interval* $C_n$ that is random, redrawn with every dataset.\n",
    "\n",
    "This makes the correct interpretation subtle. A $95\\%$ confidence interval does *not* mean \"the true $\\theta$ lies in this particular interval with probability $95\\%$\": that particular interval is already drawn, and $\\theta$ either is or is not inside it. The right reading is *about the procedure*: if we generated many intervals this way, $95\\%$ of them would contain $\\theta$. The guarantee is on the long-run hit rate of the procedure, not on any single interval. the figure draws the reading literally: one hundred intervals from one hundred independent datasets, almost all trapping the fixed truth, a handful missing; nothing about any single interval announces which kind it is. This frequentist guarantee should not be confused with the Bayesian *credible interval*, which treats $\\theta$ itself as random with a prior (the machinery behind MAP estimation in that section) and therefore *can* assert \"$\\theta$ lies in this particular interval with probability $95\\%$,\" at the price of that probability depending on the chosen prior. For the Beta posterior of that section, for instance, a $95\\%$ credible interval for a coin's bias is read directly off the posterior's quantiles ($2.5\\%$ and $97.5\\%$ of $\\mathrm{Beta}(\\alpha+x,\\ \\beta+n-x)$), with no sampling-distribution argument needed.\n",
    "\n",
    "![One hundred $95\\%$ confidence intervals, each computed from its own independent dataset, plotted against the fixed true mean (dashed line). Most intervals trap the truth; the orange ones miss. The $95\\%$ guarantee is a statement about this long-run hit rate of the procedure, not about any single interval; a given interval carries no marker of whether it traps the truth.](../img/mdl-prob-coverage.svg)\n",
    "\n",
    "### A Gaussian Example\n",
    "\n",
    "The classic case is the mean of a Gaussian $\\mathcal{N}(\\mu,\\sigma^2)$ with both parameters unknown. From $n$ samples we form the usual estimators $\\hat\\mu_n=\\frac1n\\sum_i x_i$ and the unbiased $\\hat\\sigma_n^2=\\frac1{n-1}\\sum_i (x_i-\\hat\\mu_n)^2$ from the equation. The *studentized* statistic, centered and scaled by the *estimated* standard error $\\hat\\sigma_n/\\sqrt n$ in place of the unknown true one,\n",
    "\n",
    "$$\n",
    "T = \\frac{\\hat\\mu_n - \\mu}{\\hat\\sigma_n/\\sqrt n} ,\n",
    "$$\n",
    "\n",
    "follows *Student's $t$-distribution* on $n-1$ degrees of freedom [@Student.1908], a fact we take on faith here; the distribution itself is constructed in that section. As $n\\to\\infty$ the $t$-distribution approaches a standard Gaussian, and the attribution of that limit matters. For Gaussian data the numerator $\\hat\\mu_n-\\mu$ is *exactly* Gaussian at every $n$; the limit comes from the denominator, since $\\hat\\sigma_n$ is consistent for $\\sigma$ (the law of large numbers the equation again, applied to squared deviations) and replacing a factor that converges in probability to a constant leaves a limiting distribution unchanged (*Slutsky's theorem*, also taken on faith [@Wasserman.2013]). The *central limit theorem* of that section enters for *non-Gaussian* data: it makes the sampling distribution of the mean asymptotically Gaussian anyway, so the interval below remains asymptotically valid for the mean of any finite-variance distribution, which is why it gets used so freely in practice. For large $n$, then, $T$ lands in $[-1.96, 1.96]$ with probability $\\approx95\\%$ (the Gaussian's central $95\\%$): exactly $95\\%$ in the Gaussian limit, and slightly *less* at finite $n$, where the exact $t$-distribution has heavier tails and one should use the wider $t$-quantile in its place. Rearranging $-1.96 \\le T \\le 1.96$ for $\\mu$ yields the interval\n",
    "\n",
    "$$\n",
    "\\left[\\hat\\mu_n - 1.96\\,\\frac{\\hat\\sigma_n}{\\sqrt n},\\; \\hat\\mu_n + 1.96\\,\\frac{\\hat\\sigma_n}{\\sqrt n}\\right].\n",
    "$$\n",
    "\n",
    "This is one of the most-used formulas in statistics. The half-width $1.96\\,\\hat\\sigma_n/\\sqrt n$ shrinks like $1/\\sqrt n$: to halve the interval we need *four times* the data. Let us construct one for a standard-normal sample, taking the asymptotic $t_\\star=1.96$."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "statistics-confidence-interval"
   },
   "source": [
    "import numpy as onp\n",
    "N = 1000\n",
    "samples = onp.random.normal(loc=0, scale=1, size=(N,))\n",
    "t_star = 1.96  # asymptotic value; small N would look this up in a t-table\n",
    "mu_hat = onp.mean(samples)\n",
    "se = samples.std(ddof=1) / onp.sqrt(N)  # ddof=1: unbiased sigma_hat\n",
    "(mu_hat - t_star * se, mu_hat + t_star * se)"
   ],
   "execution_count": 7,
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(np.float64(-0.059154483584664394), np.float64(0.06315527195971658))"
     },
     "metadata": {},
     "execution_count": 7
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-8 fw=all -->\n",
    "\n",
    "The interval is narrow and brackets the true mean $0$, as it should roughly $95\\%$ of the time. The same $1/\\sqrt n$ scaling shows up everywhere uncertainty is reported (error bars on a learning curve, the spread of accuracies across random seeds), and the equation is the formula behind them.\n",
    "\n",
    "\"Roughly $95\\%$ of the time\" is the Neyman guarantee the equation, and we can audit it. The cell below repeats the whole procedure one thousand times: draw a fresh dataset, build the interval the equation, and check whether it traps the true mean. The hit count should land near $950$."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "mdl-statistics-a-gaussian-example-1"
   },
   "source": [
    "rng = onp.random.default_rng(0)\n",
    "mu, n, trials = 0.0, 100, 1000\n",
    "data = rng.normal(mu, 1.0, (trials, n))           # one row = one fresh dataset\n",
    "mu_hat = data.mean(axis=1)\n",
    "se = data.std(axis=1, ddof=1) / onp.sqrt(n)\n",
    "hits = (onp.abs(mu_hat - mu) <= 1.96 * se)        # interval traps the truth?\n",
    "print(f'{hits.sum()} of {trials} intervals contain the true mean '\n",
    "      f'(expected about {0.95 * trials:.0f})')"
   ],
   "execution_count": 8,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "937 of 1000 intervals contain the true mean (expected about 950)\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-9 fw=all -->\n",
    "\n",
    "This run counts $937$ hits. Its true target is in fact $\\approx947$ rather than $950$: at $n=100$ the exact $t$-quantile is a touch wider than $1.96$, the finite-$n$ effect noted above, so the procedure's true coverage here is $94.7\\%$. The printed $937$ sits $1.4$ standard deviations below that target (the count is binomial with standard deviation $\\approx7$), well within sampling noise. This is the coverage-strip picture of the figure rendered as a count, and it is the only sense in which any confidence interval is ever \"$95\\%$ sure.\"\n",
    "\n",
    "**Propagating error bars: the delta method.** Often the number we report is not $\\hat\\theta$ itself but a smooth function of it: the log-odds of an accuracy, the ratio of two losses, a percentage change. A first-order Taylor expansion, $g(\\hat\\theta)\\approx g(\\theta)+g'(\\theta)\\,(\\hat\\theta-\\theta)$, says the fluctuation of $g(\\hat\\theta)$ is the fluctuation of $\\hat\\theta$ scaled by the local slope, so standard errors propagate as\n",
    "\n",
    "$$\n",
    "\\operatorname{se}\\bigl(g(\\hat\\theta)\\bigr) \\;\\approx\\; \\bigl|g'(\\hat\\theta)\\bigr|\\,\\operatorname{se}(\\hat\\theta),\n",
    "$$\n",
    "\n",
    "the *delta method* [@Wasserman.2013]. For instance, an accuracy of $\\hat p=0.90$ measured on $n=1{,}000$ test examples has $\\operatorname{se}(\\hat p)=\\sqrt{\\hat p(1-\\hat p)/n}\\approx 0.0095$; reporting the log-odds $g(\\hat p)=\\log\\frac{\\hat p}{1-\\hat p}$ instead, the slope $g'(\\hat p)=1/(\\hat p(1-\\hat p))\\approx 11.1$ stretches the error bar to $\\approx 0.105$ on the log-odds scale. The same rule with a gradient in place of $g'$ handles vector parameters, and when no derivative is convenient the bootstrap below sidesteps the calculus entirely.\n",
    "\n",
    "### The Bootstrap\n",
    "\n",
    "The Gaussian interval the equation rests on a special feature of the mean: its sampling distribution is known, so its standard error has a closed form, $\\hat\\sigma_n/\\sqrt n$. Most quantities a practitioner actually cares about have no such formula. What is the standard error of a *median*, a *correlation*, a model's *test accuracy*, its *AUC* (the area under the ROC curve: the probability that a random positive example is scored above a random negative one) or its *BLEU* score (the sequence-overlap metric of that section)? These are complicated functions of the data with no textbook sampling distribution, and writing down their standard error analytically ranges from painful to impossible.\n",
    "\n",
    "The **bootstrap**, introduced by Bradley Efron [@Efron.1979], escapes this with a single substitution. We never had access to the true distribution $F$ that generated our $n$ data points; if we did, we could simulate the sampling distribution of a chosen statistic\n",
    "by drawing fresh datasets from $F$ and recomputing it. The bootstrap's move, the **plug-in principle**, is to substitute the *empirical* distribution $\\hat F_n$, which puts mass $1/n$ on each observed point, for the unknown $F$. Drawing $n$ points from $\\hat F_n$ is exactly *resampling our own data $n$ times with replacement*. Concretely:\n",
    "\n",
    "1. From the original sample of size $n$, draw $n$ points **with replacement** to form a *bootstrap resample*; some points appear several times, others not at all.\n",
    "2. Compute the statistic $\\hat\\theta^*$ on the resample.\n",
    "3. Repeat $B$ times to obtain $\\hat\\theta^*_1,\\ldots,\\hat\\theta^*_B$.\n",
    "\n",
    "The spread of these $B$ replicates approximates the sampling distribution of $\\hat\\theta$, so their standard deviation estimates its standard error, and the $\\alpha/2$ and $1-\\alpha/2$ empirical percentiles of $\\{\\hat\\theta^*_b\\}$ form a *percentile* confidence interval. How many replicates? A few hundred ($B\\approx200$) suffice for a standard error, but the percentile interval rests on estimated tail quantiles, so use $B$ of at least $1{,}000$--$2{,}000$; resampling is cheap, and below we simply take $B=10{,}000$. the figure shows the construction: one original sample fans out into many resamples, whose statistics pile up into a histogram standing in for the true (unknowable) sampling distribution, with the central band cut off at those percentiles. Two caveats matter. This resampling distribution is *centered at\n",
    "$\\hat\\theta$, not at $\\theta$*; the bootstrap estimates shape and width from\n",
    "the one sample we have. It also is not valid for every statistic or sampling\n",
    "scheme. Ordinary iid bootstrap theory works well for many smooth statistics,\n",
    "but extrema, parameters on a boundary, very small samples, and nonsmooth or\n",
    "non-identifiable problems can fail; dependent observations require block,\n",
    "cluster, or other structure-preserving resampling. Exercise 1 shows the classic\n",
    "failure for a sample maximum.\n",
    "\n",
    "![The bootstrap. From a single observed sample (top) we draw many resamples of the same size *with replacement* (middle); recomputing the statistic $\\hat\\theta$ on each gives the replicates $\\hat\\theta^\\ast_b$, whose histogram (bottom) approximates the sampling distribution. Its spread estimates the standard error, and the central $1-\\alpha$ percentile band is a confidence interval. The resampling distribution is centered at $\\hat\\theta$, the dashed estimate, rather than at the unknown true $\\theta$.](../img/mdl-prob-bootstrap.svg)\n",
    "\n",
    "Let us bootstrap a statistic with no closed-form standard error, the **median**, from a skewed sample, and contrast it with the Gaussian machinery. We index the data with a matrix of random positions to draw $B$ resamples at once."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "statistics-bootstrap"
   },
   "source": [
    "rng = onp.random.default_rng(0)\n",
    "data = rng.exponential(scale=1.0, size=200)  # skewed: median != mean, no SE formula\n",
    "n = len(data)\n",
    "theta_hat = onp.median(data)                 # statistic of interest\n",
    "\n",
    "B = 10000                                     # number of bootstrap resamples\n",
    "idx = rng.integers(0, n, size=(B, n))         # n positions per resample, WITH replacement\n",
    "boot = onp.median(data[idx], axis=1)          # one median per resample\n",
    "\n",
    "se_boot = boot.std(ddof=1)                    # bootstrap standard error of the median\n",
    "ci_pct = onp.percentile(boot, [2.5, 97.5])    # percentile 95% CI -- no formula needed\n",
    "print(f'sample median        = {theta_hat:.3f}')\n",
    "print(f'bootstrap SE         = {se_boot:.3f}')\n",
    "print(f'percentile 95% CI    = ({ci_pct[0]:.3f}, {ci_pct[1]:.3f})')"
   ],
   "execution_count": 9,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "sample median        = 0.808\nbootstrap SE         = 0.072\npercentile 95% CI    = (0.682, 0.955)\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-10 fw=all -->\n",
    "\n",
    "The bootstrap hands us a standard error and an interval for the median directly,\n",
    "without deriving its sampling distribution in closed form; its validity still\n",
    "rests on the regularity and iid-sampling conditions just stated. For contrast, the Gaussian formula the equation only knows how to handle the *mean*, a different target, which on this skewed data sits well above the median."
   ]
  },
  {
   "cell_type": "code",
   "metadata": {
    "id": "statistics-bootstrap-contrast"
   },
   "source": [
    "mu_hat = data.mean()\n",
    "se_mean = data.std(ddof=1) / n**0.5           # closed-form SE, but only for the mean\n",
    "ci_gauss = (mu_hat - 1.96 * se_mean, mu_hat + 1.96 * se_mean)\n",
    "print(f'Gaussian 95% CI (mean) = ({ci_gauss[0]:.3f}, {ci_gauss[1]:.3f})')"
   ],
   "execution_count": 10,
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "Gaussian 95% CI (mean) = (0.986, 1.272)\n"
     ]
    }
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-11 fw=all -->\n",
    "\n",
    "The two intervals answer different questions and do not overlap, a direct consequence of the skew. The percentile interval for the median is also slightly *asymmetric* about $\\hat\\theta$: it inherits the shape of the resampling distribution rather than forcing the symmetric $\\pm 1.96\\,\\widehat{\\operatorname{se}}$ of a Gaussian. This is exactly why the bootstrap is so widely used in machine learning: error bars on a held-out accuracy, an AUC, or a BLEU score have no closed-form standard error, but resampling the test set delivers one in a few lines [@Efron.Hastie.2016]; that section does exactly this, putting a bootstrap error bar on its MNIST classifier's test accuracy.\n",
    "\n",
    "## Summary\n",
    "\n",
    "* An *estimator* $\\hat\\theta_n$ is a function of the data; being random, it has a *sampling distribution* whose center and spread are summarized by *bias* $\\mathbb{E}[\\hat\\theta_n]-\\theta$ and *variance*. *Consistency* ($\\hat\\theta_n\\xrightarrow{P}\\theta$) follows when both shrink with $n$; the *weak law of large numbers* the equation is the sample-mean case; *efficiency* ranks unbiased estimators by their variance.\n",
    "* The *bias-variance decomposition* $\\operatorname{MSE}(\\hat\\theta_n)=\\operatorname{Bias}(\\hat\\theta_n)^2+\\operatorname{Var}(\\hat\\theta_n)$ splits the error exactly because, after centering at $\\mathbb{E}[\\hat\\theta_n]$, the cross term vanishes. This is the same U-curve as the under/overfitting trade-off (expected test error adds an irreducible noise floor $\\sigma^2$), and it explains why regularization trades bias for variance.\n",
    "* The unbiased sample variance divides by $n-1$, not $n$: estimating the mean from the same data costs one degree of freedom, and the $1/(n-1)$ factor corrects the resulting bias exactly.\n",
    "* *Hypothesis testing* weighs evidence against a null $H_0$ via a test statistic and its $p$-value $P_{H_0}(\\textrm{data this extreme})$; we control the type I error rate $\\alpha$ and want high power $1-\\beta$. A $p$-value is not $P(H_0\\mid\\textrm{data})$. A *permutation test* builds the null distribution by shuffling group labels, with no Gaussian assumptions; under many tests, control the family-wise error rate (Bonferroni) or the false discovery rate (Benjamini--Hochberg).\n",
    "* A *confidence interval* contains $\\theta$ with probability $\\ge 1-\\alpha$ over repeated datasets; the Gaussian interval $\\hat\\mu_n \\pm 1.96\\,\\hat\\sigma_n/\\sqrt n$, exact for Gaussian data in the large-$n$ limit and extended to other data by the *central limit theorem*, has half-width shrinking like $1/\\sqrt n$.\n",
    "* The ordinary iid *bootstrap* estimates the sampling distribution of many\n",
    "  regular statistics with no convenient closed-form standard error (a median,\n",
    "  an accuracy, an AUC) by resampling observations with replacement. It can fail\n",
    "  for extrema, boundaries, and nonregular statistics, and dependent data need a\n",
    "  resampling scheme that preserves their structure.\n",
    "\n",
    "## Exercises\n",
    "\n",
    "1. Let $X_1, \\ldots, X_n \\overset{\\textrm{iid}}{\\sim} \\textrm{Unif}(0,\\theta)$ and consider the estimators $\\hat\\theta = \\max\\{X_1,\\ldots,X_n\\}$ and $\\tilde\\theta = \\frac2n\\sum_i X_i$. Find the bias, variance, and MSE of each, and decide which is better. Is $\\hat\\theta$ biased? Is it consistent? Finally, try to bootstrap a confidence interval for $\\theta$ from $\\hat\\theta=\\max_i X_i$ and explain why it is poor: with what probability does a bootstrap resample contain the largest observation, so that $\\hat\\theta^* = \\hat\\theta$ exactly? (This is the classic example where the bootstrap fails: the true sampling distribution of $\\hat\\theta$ lives entirely below $\\theta$, while the resampling distribution puts a large point mass exactly at $\\hat\\theta$.)\n",
    "2. Prove the bias-variance decomposition the equation directly by expanding $\\mathbb{E}[(\\hat\\theta_n-\\theta)^2]$ into $\\mathbb{E}[\\hat\\theta_n^2]-2\\theta\\,\\mathbb{E}[\\hat\\theta_n]+\\theta^2$ and substituting $\\mathbb{E}[\\hat\\theta_n^2]=\\operatorname{Var}(\\hat\\theta_n)+\\mathbb{E}[\\hat\\theta_n]^2$. Confirm it agrees with the add-and-subtract proof in the text.\n",
    "3. The decomposition check computes the variance of the $10{,}000$ estimates with the plug-in estimator (`ddof=0`). Rerun it with the unbiased estimator (`ddof=1`) and compare the two sides of the equation again. Which variance estimator makes the identity exact to floating-point precision, and why? (Hint: the empirical distribution's expectations are plain averages; which denominator does that force?)\n",
    "4. Shrink the per-dataset size $n$ in the sampling-distribution simulation and confirm the spread of $\\hat\\theta$ widens like $\\sigma/\\sqrt n$. Repeat with the biased estimator $\\hat\\theta=\\max_i X_i$ for $\\textrm{Unif}(0,\\theta)$ and watch the center shift away from $\\theta$.\n",
    "5. A test reports $p = 0.5$. Is this evidence that $H_0$ is true? Explain in terms of $P(\\textrm{data}\\mid H_0)$ versus $P(H_0\\mid\\textrm{data})$, and describe a situation where a large $p$-value reflects only low power.\n",
    "6. Using the $1/(\\textrm{effect size})^2$ scaling, estimate how many times more samples are needed to detect an effect of size $0.1$ than one of size $0.5$ at the same $\\alpha$ and power.\n",
    "7. Run the confidence-interval code with $N=2$ and $\\alpha=0.5$ (so $t_\\star=1.0$) for $100$ independently generated datasets, and look at the resulting intervals. Some are extremely short and far from the true mean $0$. Does this contradict the $1-\\alpha$ coverage guarantee? Would you trust a short interval as a sign of a precise estimate?\n",
    "8. Modify the bootstrap code to target the *mean* instead of the median. Compare the bootstrap standard error and percentile interval against the closed-form Gaussian results $\\hat\\sigma_n/\\sqrt n$ and the equation on the same sample. Do they roughly agree? Why should they, given the central limit theorem?\n",
    "9. Suppose you compare $m=20$ models against a baseline, all in truth no better, each with an independent test at $\\alpha=0.05$. What is the probability that at least one clears $p\\le 0.05$ by chance? Recompute the per-test threshold the Bonferroni correction prescribes to hold the family-wide false-positive rate at $0.05$, and verify it brings the spurious-win probability back near $0.05$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- d2l:prose id=mdl-statistics-md-12 fw=pytorch -->\n",
    "\n",
    "[Discussions](https://d2l.discourse.group/t/1102)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "d2l-pytorch",
   "display_name": "d2l (PyTorch)",
   "language": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
