Statistical inference

In a trial of a malaria vaccine in West Africa, some children receive the vaccine and some a control injection, and the numbers of malaria cases in the two groups are counted. The counts are data; the question is about the world: how well does the vaccine protect children in general? A soil scientist in Punjab measures wheat yields on forty plots treated with a new fertiliser. A survey team in Nairobi asks four hundred households whether they use mobile money. In each case we see a finite, random sample and want to say something about an unknown quantity — a protection rate, a mean yield, a population proportion — with an honest statement of how uncertain we are.

Probability theory goes forwards: given the model (a coin with $P(\text{heads}) = p$), what data should we expect? Statistical inference goes backwards: given the data, what can we say about $p$? It does so by imagining the study repeated: an estimate computed from the data is itself a random variable, and its distribution across hypothetical repetitions — its sampling distribution — tells us how far to trust it.

This lesson builds on probability theory (expectation, variance, the law of large numbers and the central limit theorem) and on the school-level ideas of samples, bias and margins of error in statistical inference: first steps. Matrix notation from linear algebra appears in the section on regression.

Suppose $7$ of $20$ seeds germinate. If the germination probability were $p = 0.9$, seeing only $7$ would be astonishing; if it were $p = 0.05$, also astonishing. The likelihood $L(p) = p^7(1 - p)^{13}$ measures, for each candidate $p$, how probable the observed data would be. It is largest at $p = 7/20 = 0.35$: the maximum likelihood estimate.

Two features of the picture carry the whole subject. The peak is our estimate. The width tells us how well the data pin the parameter down; with $n$ observations it shrinks like $1/\sqrt n$, the same square-root law as the standard deviation of a sample mean. Confidence intervals, tests and Bayesian posteriors are different ways of turning that width into a statement of uncertainty.

So an estimator can be bad in two ways — systematically off (bias) or erratic (variance) — and reducing one may increase the other.

The standard error of a mean is therefore $\sigma/\sqrt n$: to halve it, quadruple the sample.

How should we estimate $\sigma^2$? The natural guess $\frac1n\sum(X_i - \bar X)^2$ is too small on average, because the data are closer to their own mean $\bar X$ than to the true mean $\mu$. The key is an algebraic identity.

Geometrically, the residual vector $(X_i - \bar X)_i$ is confined to the $(n - 1)$-dimensional subspace of vectors whose entries sum to zero: one "degree of freedom" was used up in estimating $\mu$ (see the viewpoints).

How good can an estimator be? Define the score $U(\theta) = \ell'(\theta) = \sum_i \partial_\theta \log f(X_i; \theta)$ and the Fisher information $I(\theta) = E_\theta\big[(\partial_\theta \log f(X_1; \theta))^2\big]$.

For the Bernoulli model, $\partial_p \log f(x; p) = \frac xp - \frac{1 - x}{1 - p}$ and $I(p) = \frac{1}{p(1 - p)}$, so the bound is $\frac{p(1-p)}{n} = \operatorname{Var}\hat p$: the sample proportion is efficient. The uniform model violates the regularity conditions (its support depends on $\theta$), and indeed its unbiased estimator has variance of order $1/n^2$, beating the "bound". Under regularity conditions the MLE is consistent and asymptotically normal: $\sqrt n(\hat\theta_n - \theta) \to N\big(0, 1/I(\theta)\big)$ in distribution (Fisher, 1925; made rigorous by Cramér, 1946) — it attains the Cramér–Rao bound in the limit, which is the main justification for using it.

The pattern — find a pivot, a function of the data and $\theta$ whose distribution does not depend on $\theta$, and invert it — gives most classical intervals:

• Unknown $\sigma$, normal data. Replacing $\sigma$ by $S$ changes the pivot to $T = \frac{\bar X - \mu}{S/\sqrt n}$, which has Student's $t$ distribution with $n - 1$ degrees of freedom (Gosset, writing as "Student", 1908; proved by Fisher). The interval $\bar X \pm t_{n-1}\,S/\sqrt n$ has exact coverage; $t_{n-1} > 1.96$ (for example $t_4 = 2.776$) accounts for the extra uncertainty of estimating $\sigma$.
• Large samples. By the central limit theorem, $\bar X \pm 1.96\,S/\sqrt n$ has coverage approximately $95\%$ for large $n$, for any distribution with finite variance.
• A proportion. The Wald interval $\hat p \pm 1.96\sqrt{\hat p(1 - \hat p)/n}$ comes from the CLT. Its actual coverage can be far below $95\%$ when $n\hat p$ is small (with $n = 10$ and $p = 0.05$ it is $40\%$); the Wilson interval (1927), obtained by inverting the test rather than plugging in $\hat p$, behaves much better. The first experiment below computes these coverages exactly.

What "95% confidence" means. Before the data are collected, the procedure has a $95\%$ chance of producing an interval that covers $p$. After the data are in, $[0.603, 0.697]$ either contains $p$ or it does not; $p$ is not random in this framework, so "$P(0.603 \le p \le 0.697) = 0.95$" is not a frequentist statement. If the survey were repeated many times, about $95\%$ of the intervals would cover the true $p$. (The Bayesian credible interval below does support a probability statement about $p$, at the price of a prior.)

The power of a test depends on the true effect, the noise and the sample size. For a one-sided $z$-test of $H_0 : \mu = 0$ at level $0.05$ with known $\sigma$, we reject when $\sqrt n\,\bar X/\sigma > 1.645$; if the true mean is $\mu = \delta\sigma$, the power is $1 - \Phi(1.645 - \delta\sqrt n)$. For an effect of half a standard deviation ($\delta = 0.5$), the power is $0.47$ with $n = 10$, $0.80$ with $n = 25$ and $0.95$ with $n = 44$.

Three facts complete the picture.

1. Tests and intervals are two views of one thing. The two-sided $z$-test rejects $H_0 : \mu = \mu_0$ at level $\alpha$ exactly when $\mu_0$ lies outside the $(1 - \alpha)$ $z$-interval. A confidence interval is the set of parameter values that the data do not reject — and it shows the size of the effect, which a p-value alone hides.
2. Not rejecting is not accepting. A study with low power will usually fail to reject even when the effect is real. "$p = 0.3$" means the data are compatible with no effect; they may be equally compatible with a large one. Report the interval.
3. Significance is not importance. With enough data, a trivially small effect becomes "significant"; with too little, an important one does not.

Run $20$ independent tests of true null hypotheses at level $0.05$. The chance that at least one is "significant" is $1 - 0.95^{20} \approx 0.64$. Search through enough variables, subgroups or outcome measures, and a spurious discovery is almost guaranteed — and if only the significant results are reported, the literature fills with false positives. The family-wise error rate (FWER) is the probability of at least one false rejection.

For independent tests the exact FWER is $1 - (1 - \alpha/m)^m$, and Bernoulli's inequality $(1 - a)^m \ge 1 - ma$ shows how close the union bound is. Bonferroni is safe but conservative: with many tests it loses power. The Benjamini–Hochberg procedure (1995) controls instead the false discovery rate, the expected proportion of false discoveries among all discoveries: sort the p-values $p_{(1)} \le \cdots \le p_{(m)}$, find the largest $k$ with $p_{(k)} \le k\alpha/m$, and reject the $k$ smallest. For independent tests it guarantees false discovery rate at most $\alpha m_0/m \le \alpha$, where $m_0$ is the number of true nulls (stated here without proof). The second experiment compares both procedures by simulation.

Multiplicity also hides in choices that are never written down — which outcome, which subgroup, which covariates, when to stop collecting data. Pre-registering the analysis, reporting every test performed, and publishing code, data and random seeds so that others can rerun the analysis (see computational reproducibility) are the practical defences.

In the linear model $y = X\beta + \varepsilon$, the $n \times p$ matrix $X$ holds the explanatory variables (a first column of ones gives an intercept), $\beta \in \R^p$ is unknown, and $\varepsilon$ is random noise. The least-squares estimate minimises the residual sum of squares $\|y - X\beta\|^2 = \sum_i \big(y_i - (X\beta)_i\big)^2$.

The fitted values $\hat y = X\hat\beta$ are the orthogonal projection of $y$ onto the column space of $X$ — the viewpoints section takes this picture further. For a straight line $y = \beta_0 + \beta_1 x$, the two normal equations $\sum_i(y_i - \beta_0 - \beta_1x_i) = 0$ and $\sum_i x_i(y_i - \beta_0 - \beta_1x_i) = 0$ give

so the line passes through $(\bar x, \bar y)$, and the residuals sum to zero and are uncorrelated with $x$.

Statistical properties. If $X$ is fixed with independent columns, $E[\varepsilon] = 0$ and $\operatorname{Cov}(\varepsilon) = \sigma^2 I$, then $E[\hat\beta] = (X^\top X)^{-1}X^\top X\beta = \beta$ (unbiased) and $\operatorname{Cov}(\hat\beta) = \sigma^2(X^\top X)^{-1}$. The Gauss–Markov theorem (stated without proof) says $\hat\beta$ has the smallest variance among all linear unbiased estimators. The residual sum of squares divided by $n - p$ is unbiased for $\sigma^2$ — Bessel's correction is the special case $p = 1$, $X$ a column of ones. When an intercept is included, $\sum(y_i - \bar y)^2 = \sum(\hat y_i - \bar y)^2 + \sum(y_i - \hat y_i)^2$ (Pythagoras), and $R^2 = 1 - \frac{\sum(y_i - \hat y_i)^2}{\sum(y_i - \bar y)^2}$ is the fraction of variation "explained"; for a straight line, $R^2 = r^2$.

Check the residuals. Every one of these properties rests on assumptions — linearity, constant variance, independence, and for intervals and tests often normality of the noise. They are checked by looking at the residuals $y_i - \hat y_i$: plotted against the fitted values they should show no pattern (a curve signals non-linearity, a funnel signals non-constant variance); a normal quantile plot checks normality; and points with large leverage $h_{ii}$ (the diagonal of the projection matrix $X(X^\top X)^{-1}X^\top$) deserve scrutiny, because they pull the line towards themselves. Anscombe's quartet (1973) makes the point unforgettable: four data sets with the same means, variances, correlation ($0.816$) and fitted line $y = 3.00 + 0.500x$ — one roughly linear, one a smooth curve, one linear with a single outlier, and one where a single point determines the slope. And no residual check turns correlation into causation: that needs design, such as randomisation.

The Bayesian approach treats the unknown parameter as random, with a prior distribution $\pi(\theta)$ expressing what is believed before the data. Bayes' theorem gives the posterior

and all conclusions are read off it: a posterior mean as an estimate, an interval of posterior probability $0.95$ (a credible interval) as a statement of uncertainty.

The posterior mean is a compromise between the prior mean and the data, and the data's weight $1 - w$ tends to $1$ as $n$ grows. With the uniform prior Beta$(1, 1)$ it is $\frac{k + 1}{n + 2}$, Laplace's rule of succession. For the survey above ($260$ of $400$), the uniform prior gives the posterior Beta$(261, 141)$, with mean $0.649$ and central $95\%$ credible interval $[0.602, 0.695]$ — numerically almost the same as the confidence intervals, with a different meaning: given the model and the prior, $P(0.602 \le p \le 0.695 \mid \text{data}) = 0.95$. This agreement is no accident: under regularity conditions, the Bernstein–von Mises theorem says that for large $n$ the posterior is approximately $N\big(\hat\theta, 1/(nI(\hat\theta))\big)$, whatever the (smooth, positive) prior. With little data, or many parameters, the prior matters, and should be reported and examined like any other assumption.

Mathlib has no "statistics" library as such, but its probability theory — random variables as functions X : Ω → ℝ, expectation μ[X], variance Var[X; μ], independence X ⟂ᵢ[μ] Y — is exactly what the theory of estimation needs, and its linear algebra handles least squares. Data sets are functions summed over a Finset.

Estimators. The sample mean as a definition; the shift identity; unbiasedness and the variance $\sigma^2/n$ of the sample mean; and Bessel's correction, proved as in the lesson. The hypotheses are pairwise independence (Mathlib's form of "uncorrelated" in IndepFun.variance_sum) with common mean and variance:

/-- The **sample mean** `x̄ = (∑ᵢ xᵢ)/|s|` of the data `xᵢ`, `i ∈ s`. (For the empty data set
Lean's convention `a / 0 = 0` gives `0`.) -/
noncomputable def sampleMean {ι : Type*} (s : Finset ι) (x : ι → ℝ) : ℝ :=
  (∑ i ∈ s, x i) / s.card

/-- The total of the data is `n · x̄`. -/
theorem sum_eq_card_mul_mean {ι : Type*} (s : Finset ι) (x : ι → ℝ) :
    ∑ i ∈ s, x i = s.card * sampleMean s x := by
  unfold sampleMean
  rcases Nat.eq_zero_or_pos s.card with h | h
  · rw [Finset.card_eq_zero.mp h]
    simp
  · have : (s.card : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr h.ne'
    field_simp

/-- **The shift identity.** For any number `c`,
`∑ (xᵢ − x̄)² = ∑ (xᵢ − c)² − n (x̄ − c)²`. In particular (taking `c` arbitrary) the sample
mean minimises `c ↦ ∑ (xᵢ − c)²`, and (taking `c = μ`) it is the key step in Bessel's
correction. Proof: expand `((xᵢ − c) − (x̄ − c))²` and use `∑ (xᵢ − c) = n (x̄ − c)`. -/
theorem sum_sq_dev_eq {ι : Type*} (s : Finset ι) (x : ι → ℝ) (c : ℝ) :
    ∑ i ∈ s, (x i - sampleMean s x) ^ 2 =
      ∑ i ∈ s, (x i - c) ^ 2 - s.card * (sampleMean s x - c) ^ 2 := by
  set m := sampleMean s x
  have hsum : ∑ i ∈ s, (x i - c) = s.card * (m - c) := by
    rw [sum_sub_distrib, sum_eq_card_mul_mean s x, sum_const, nsmul_eq_mul]
    ring
  have hpt : ∀ i ∈ s, (x i - m) ^ 2 = (x i - c) ^ 2 - 2 * (m - c) * (x i - c) + (m - c) ^ 2 := by
    intro i _
    ring
  rw [sum_congr rfl hpt, sum_add_distrib, sum_sub_distrib, ← mul_sum, hsum, sum_const,
    nsmul_eq_mul]
  ring

/-- **The sample mean is unbiased:** if each `Xᵢ` has mean `m`, then `E[X̄ₙ] = m` (`n ≥ 1`).
Only linearity of expectation is used — no independence. -/
theorem sample_mean_unbiased (hint : ∀ i, Integrable (X i) μ) (hmean : ∀ i, μ[X i] = m)
    {n : ℕ} (hn : 0 < n) : μ[fun ω => sampleMean (range n) (fun i => X i ω)] = m := by
  have hnR : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr hn.ne'
  simp only [sampleMean, card_range]
  rw [integral_div, integral_finsetSum _ fun i _ => hint i]
  simp only [hmean, sum_const, card_range, nsmul_eq_mul]
  field_simp

/-- **The variance of the sample mean is `σ²/n`** for pairwise independent observations with
common variance `σ² = v`: variances of independent variables add (`IndepFun.variance_sum`),
and `Var(X/n) = Var(X)/n²`. -/
theorem variance_sample_mean (hL2 : ∀ i, MemLp (X i) 2 μ)
    (hindep : Pairwise fun i j => X i ⟂ᵢ[μ] X j) (hvar : ∀ i, Var[X i; μ] = v) (n : ℕ) :
    Var[fun ω => sampleMean (range n) (fun i => X i ω); μ] = v / n := by
  have h : (fun ω => sampleMean (range n) (fun i => X i ω)) =
      fun ω => (n : ℝ)⁻¹ * (∑ i ∈ range n, X i) ω := by
    funext ω
    simp [sampleMean, Finset.sum_apply, div_eq_inv_mul]
  rw [h, variance_const_mul,
    IndepFun.variance_sum (fun i _ => hL2 i) (fun i _ j _ hij => hindep hij)]
  simp only [hvar, sum_const, card_range, nsmul_eq_mul]
  rcases Nat.eq_zero_or_pos n with rfl | hn
  · simp
  · have hnR : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr hn.ne'
    field_simp

/-- The sample mean of finitely many observations with finite variance has finite variance. -/
theorem sample_mean_memLp [IsProbabilityMeasure μ] (hL2 : ∀ i, MemLp (X i) 2 μ) (n : ℕ) :
    MemLp (fun ω => sampleMean (range n) (fun i => X i ω)) 2 μ := by
  have hsum : MemLp (∑ i ∈ range n, X i) 2 μ := memLp_finsetSum' _ fun i _ => hL2 i
  have := hsum.const_mul ((n : ℝ)⁻¹)
  convert this using 1
  funext ω
  simp [sampleMean, Finset.sum_apply, div_eq_inv_mul]

/-- **Bessel's correction.** For pairwise independent observations with common mean `m` and
common variance `v`, and `n ≥ 1`, `E[∑ (Xᵢ − X̄)²] = (n − 1) v`. Hence the sample variance
`S² = ∑ (Xᵢ − X̄)²/(n − 1)` is an unbiased estimator of `v` (for `n ≥ 2`), while dividing by
`n` underestimates it on average. Proof, as in the lesson: by the shift identity with `c = m`,
`∑ (Xᵢ − X̄)² = ∑ (Xᵢ − m)² − n (X̄ − m)²`; the expectations are `n v` and `n · v/n`. -/
theorem bessel [IsProbabilityMeasure μ] (hL2 : ∀ i, MemLp (X i) 2 μ) (hindep : Pairwise fun i j => X i ⟂ᵢ[μ] X j)
    (hmean : ∀ i, μ[X i] = m) (hvar : ∀ i, Var[X i; μ] = v) {n : ℕ} (hn : 0 < n) :
    ∫ ω, ∑ i ∈ range n, (X i ω - sampleMean (range n) (fun j => X j ω)) ^ 2 ∂μ =
      (n - 1) * v := by
  set Xbar : Ω → ℝ := fun ω => sampleMean (range n) (fun j => X j ω) with hXbar
  have hnR : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr hn.ne'
  -- Pointwise, the shift identity with `c = m`.
  have hpt : (fun ω => ∑ i ∈ range n, (X i ω - Xbar ω) ^ 2) =
      fun ω => ∑ i ∈ range n, (X i ω - m) ^ 2 - n * (Xbar ω - m) ^ 2 := by
    funext ω
    rw [hXbar]
    simp only
    rw [sum_sq_dev_eq (range n) (fun j => X j ω) m, card_range]
  -- Integrability of the two pieces.
  have hint_i : ∀ i, Integrable (fun ω => (X i ω - m) ^ 2) μ := fun i =>
    ((hL2 i).sub (memLp_const m)).integrable_sq
  have hint_bar : Integrable (fun ω => (Xbar ω - m) ^ 2) μ :=
    ((sample_mean_memLp hL2 n).sub (memLp_const m)).integrable_sq
  -- Each `E[(Xᵢ − m)²] = v`, and `E[(X̄ − m)²] = Var X̄ = v/n`.
  have hEi : ∀ i, ∫ ω, (X i ω - m) ^ 2 ∂μ = v := by
    intro i
    rw [← hvar i, variance_eq_integral (hL2 i).aestronglyMeasurable.aemeasurable, hmean i]
  have hEbar : ∫ ω, (Xbar ω - m) ^ 2 ∂μ = v / n := by
    have hm : μ[Xbar] = m :=
      sample_mean_unbiased (fun i => (hL2 i).integrable one_le_two) hmean hn
    rw [← variance_sample_mean hL2 hindep hvar n, ← hXbar,
      variance_eq_integral (sample_mean_memLp hL2 n).aestronglyMeasurable.aemeasurable, hm]
  show ∫ ω, (fun ω => ∑ i ∈ range n, (X i ω - Xbar ω) ^ 2) ω ∂μ = (n - 1) * v
  rw [hpt, integral_sub (integrable_finsetSum _ fun i _ => hint_i i) (hint_bar.const_mul _),
    integral_finsetSum _ fun i _ => hint_i i, integral_const_mul, hEbar]
  simp only [hEi, sum_const, card_range, nsmul_eq_mul]
  field_simp

Consistency of the sample mean (Mathlib's Chebyshev inequality applied to $\bar X_n$), the computing formula for the variance, and the bias–variance decomposition:

/-- **Consistency of the sample mean**, with a rate: `P(|X̄ₙ − m| ≥ ε) ≤ v / (n ε²)`. Proof:
Mathlib's Chebyshev inequality (`meas_ge_le_variance_div_sq`) applied to `X̄ₙ`, whose mean is
`m` and whose variance is `v/n`. As `n → ∞` the bound tends to `0`. -/
theorem sample_mean_consistent [IsProbabilityMeasure μ] (hL2 : ∀ i, MemLp (X i) 2 μ)
    (hindep : Pairwise fun i j => X i ⟂ᵢ[μ] X j) (hmean : ∀ i, μ[X i] = m)
    (hvar : ∀ i, Var[X i; μ] = v) {n : ℕ} (hn : 0 < n) {ε : ℝ} (hε : 0 < ε) :
    μ {ω | ε ≤ |sampleMean (range n) (fun i => X i ω) - m|} ≤
      ENNReal.ofReal (v / (n * ε ^ 2)) := by
  have h := meas_ge_le_variance_div_sq (sample_mean_memLp hL2 n) hε
  rw [sample_mean_unbiased (fun i => (hL2 i).integrable one_le_two) hmean hn,
    variance_sample_mean hL2 hindep hvar n, div_div] at h
  exact h

/-- The computing formula `Var X = E[X²] − (E X)²`, for `X` with finite variance. This is
Mathlib's `ProbabilityTheory.variance_eq_sub`. -/
theorem variance_computing_formula [IsProbabilityMeasure μ] {X : Ω → ℝ} (hX : MemLp X 2 μ) :
    Var[X; μ] = μ[X ^ 2] - μ[X] ^ 2 :=
  variance_eq_sub hX

/-- **Mean squared error = variance + bias².** For an estimator `T` with finite variance and
a parameter value `θ`, `E[(T − θ)²] = Var T + (E T − θ)²`. Proof: write
`T − θ = (T − E T) + (E T − θ)`; the cross term has expectation `0`. -/
theorem mse_decomposition [IsProbabilityMeasure μ] {T : Ω → ℝ} (hT : MemLp T 2 μ) (θ : ℝ) :
    ∫ ω, (T ω - θ) ^ 2 ∂μ = Var[T; μ] + (μ[T] - θ) ^ 2 := by
  set a := μ[T]
  have hTi : Integrable T μ := hT.integrable one_le_two
  have hpt : (fun ω => (T ω - θ) ^ 2) =
      fun ω => (T ω - a) ^ 2 + (2 * (a - θ) * (T ω - a) + (a - θ) ^ 2) := by
    funext ω
    ring
  have h1 : Integrable (fun ω => (T ω - a) ^ 2) μ := (hT.sub (memLp_const a)).integrable_sq
  have h2 : Integrable (fun ω => 2 * (a - θ) * (T ω - a)) μ :=
    (hTi.sub (integrable_const a)).const_mul _
  have h3 : Integrable (fun ω => 2 * (a - θ) * (T ω - a) + (a - θ) ^ 2) μ :=
    h2.add (integrable_const _)
  rw [hpt, integral_add h1 h3, integral_add h2 (integrable_const _),
    integral_const_mul, integral_sub hTi (integrable_const a),
    variance_eq_integral hT.aestronglyMeasurable.aemeasurable]
  simp [a]

Likelihood. The calculus-free proof that $k/n$ maximises the Bernoulli log-likelihood, using Mathlib's Real.log_le_sub_one_of_pos:

/-- **The Bernoulli maximum likelihood estimator is `k/n`.** With `k` successes in `n`
independent trials (`0 < k < n`), the log-likelihood `ℓ(p) = k log p + (n − k) log(1 − p)` is
largest at `p̂ = k/n`. Proof without calculus, as in the lesson: with `q = k/n`,
`ℓ(p) − ℓ(q) = k log(p/q) + (n − k) log((1 − p)/(1 − q))
  ≤ k (p/q − 1) + (n − k)((1 − p)/(1 − q) − 1) = 0`, using `log x ≤ x − 1`
(Mathlib's `Real.log_le_sub_one_of_pos`). -/
theorem bernoulli_mle {n k : ℕ} (hk : 0 < k) (hkn : k < n) {p : ℝ} (hp0 : 0 < p) (hp1 : p < 1) :
    k * Real.log p + (n - k) * Real.log (1 - p) ≤
      k * Real.log (k / n) + (n - k) * Real.log (1 - k / n) := by
  have hkR : (0 : ℝ) < k := by exact_mod_cast hk
  have hknR : (k : ℝ) < n := by exact_mod_cast hkn
  have hnR : (0 : ℝ) < n := hkR.trans hknR
  set q : ℝ := k / n with hq
  have hq0 : 0 < q := div_pos hkR hnR
  have hq1 : q < 1 := (div_lt_one hnR).mpr hknR
  have h1 : Real.log p - Real.log q ≤ p / q - 1 := by
    rw [← Real.log_div hp0.ne' hq0.ne']
    exact Real.log_le_sub_one_of_pos (div_pos hp0 hq0)
  have h2 : Real.log (1 - p) - Real.log (1 - q) ≤ (1 - p) / (1 - q) - 1 := by
    rw [← Real.log_div (by linarith) (by linarith)]
    exact Real.log_le_sub_one_of_pos (div_pos (by linarith) (by linarith))
  have hzero : (k : ℝ) * (p / q - 1) + (n - k) * ((1 - p) / (1 - q) - 1) = 0 := by
    have hnk : (n : ℝ) - k ≠ 0 := by linarith
    rw [hq]
    field_simp
    ring
  have e1 := mul_le_mul_of_nonneg_left h1 hkR.le
  have e2 := mul_le_mul_of_nonneg_left h2 (show (0 : ℝ) ≤ n - k by linarith)
  linarith

Least squares. Both directions of the normal-equations theorem, for any real matrix X : Matrix n p ℝ, with ⬝ᵥ the dot product and *ᵥ matrix–vector multiplication:

/-- **Least squares via the normal equations.** If `b` satisfies `(XᵀX) b = Xᵀy`, then `b`
minimises the residual sum of squares `‖y − Xβ‖²` over all `β`. Proof (Pythagoras): with
`r = y − Xb` and `w = b − β`, `y − Xβ = r + Xw`, and the cross term `r · Xw = (Xᵀr) · w`
vanishes because `Xᵀr = Xᵀy − XᵀXb = 0`; so `‖y − Xβ‖² = ‖r‖² + ‖Xw‖² ≥ ‖r‖²`. -/
theorem normal_equations_minimize {n p : Type*} [Fintype n] [Fintype p] (X : Matrix n p ℝ)
    (y : n → ℝ) (b : p → ℝ) (hb : (Xᵀ * X) *ᵥ b = Xᵀ *ᵥ y) (β : p → ℝ) :
    (y - X *ᵥ b) ⬝ᵥ (y - X *ᵥ b) ≤ (y - X *ᵥ β) ⬝ᵥ (y - X *ᵥ β) := by
  set r := y - X *ᵥ b
  set w := b - β
  have hr : Xᵀ *ᵥ r = 0 := by
    simp only [r, mulVec_sub, mulVec_mulVec, hb, sub_self]
  have hcross : r ⬝ᵥ (X *ᵥ w) = 0 := by
    rw [dotProduct_mulVec, ← mulVec_transpose, hr, zero_dotProduct]
  have hsplit : y - X *ᵥ β = r + X *ᵥ w := by
    simp only [r, w, mulVec_sub]
    abel
  have hsq : 0 ≤ (X *ᵥ w) ⬝ᵥ (X *ᵥ w) :=
    Finset.sum_nonneg fun i _ => mul_self_nonneg _
  rw [hsplit, add_dotProduct, dotProduct_add, dotProduct_add, hcross, dotProduct_comm (X *ᵥ w) r,
    hcross]
  linarith

/-- **Conversely, every least-squares solution satisfies the normal equations.** If `b`
minimises `‖y − Xβ‖²`, then `(XᵀX) b = Xᵀy`. Proof: for each direction `w` and real `t`,
`‖y − X(b + tw)‖² = ‖r‖² − 2t (Xᵀr)·w + t² ‖Xw‖²` must be at least `‖r‖²`; a quadratic
`−2ta + t²c` with `c ≥ 0` that is never negative has `a = 0`. Taking `w = Xᵀr` gives
`‖Xᵀr‖² = 0`. -/
theorem minimizer_normal_equations {n p : Type*} [Fintype n] [Fintype p] (X : Matrix n p ℝ)
    (y : n → ℝ) (b : p → ℝ)
    (hmin : ∀ β, (y - X *ᵥ b) ⬝ᵥ (y - X *ᵥ b) ≤ (y - X *ᵥ β) ⬝ᵥ (y - X *ᵥ β)) :
    (Xᵀ * X) *ᵥ b = Xᵀ *ᵥ y := by
  set r := y - X *ᵥ b
  set g := Xᵀ *ᵥ r
  -- Expand the residual sum of squares at `b + t g`.
  have hexp : ∀ t : ℝ, (y - X *ᵥ (b + t • g)) ⬝ᵥ (y - X *ᵥ (b + t • g)) =
      r ⬝ᵥ r - 2 * t * (g ⬝ᵥ g) + t ^ 2 * ((X *ᵥ g) ⬝ᵥ (X *ᵥ g)) := by
    intro t
    have hsplit : y - X *ᵥ (b + t • g) = r - t • (X *ᵥ g) := by
      simp only [r, mulVec_add, mulVec_smul]
      abel
    have hrg : r ⬝ᵥ (X *ᵥ g) = g ⬝ᵥ g := by
      rw [dotProduct_mulVec, ← mulVec_transpose]
    rw [hsplit]
    simp only [sub_dotProduct, dotProduct_sub, smul_dotProduct, dotProduct_smul, smul_eq_mul,
      hrg, dotProduct_comm (X *ᵥ g) r]
    ring
  set a := g ⬝ᵥ g
  set c := (X *ᵥ g) ⬝ᵥ (X *ᵥ g)
  have hc : 0 ≤ c := Finset.sum_nonneg fun i _ => mul_self_nonneg _
  have ha : 0 ≤ a := Finset.sum_nonneg fun i _ => mul_self_nonneg _
  -- Choose `t = a/(c + 1)`: the quadratic becomes `−a²(c + 2)/(c + 1)² ≤ 0`, so `a = 0`.
  have key := hmin (b + (a / (c + 1)) • g)
  rw [hexp] at key
  have hc1 : 0 < c + 1 := by linarith
  have hq : -2 * (a / (c + 1)) * a + (a / (c + 1)) ^ 2 * c = -(a ^ 2 * (c + 2)) / (c + 1) ^ 2 := by
    field_simp
    ring
  have h0 : 0 ≤ -(a ^ 2 * (c + 2)) / (c + 1) ^ 2 := by
    rw [← hq]
    linarith
  have h0' : a ^ 2 * (c + 2) ≤ 0 := by
    by_contra hneg
    have hpos : 0 < a ^ 2 * (c + 2) := not_le.mp hneg
    have : -(a ^ 2 * (c + 2)) / (c + 1) ^ 2 < 0 := div_neg_of_neg_of_pos (by linarith) (by positivity)
    linarith
  have ha0 : a = 0 := by
    have ha2 : a ^ 2 = 0 := by nlinarith [sq_nonneg a]
    exact pow_eq_zero_iff (two_ne_zero) |>.mp ha2
  -- `g · g = 0` forces `g = 0`, i.e. `Xᵀ(y − Xb) = 0`.
  have hg : g = 0 := by
    funext i
    have hsum : ∑ j, g j * g j = 0 := ha0
    have := (Finset.sum_eq_zero_iff_of_nonneg fun j _ => mul_self_nonneg (g j)).mp hsum i
      (Finset.mem_univ i)
    simpa using this
  have : Xᵀ *ᵥ y - (Xᵀ * X) *ᵥ b = 0 := by
    rw [← mulVec_mulVec, ← mulVec_sub]
    exact hg
  exact (sub_eq_zero.mp this).symm