Variance — MrCLean

The previous chapter established that difference in means is unbiased under complete randomization. Unbiasedness is a one-line fact about a first moment, and it is cheap: many bad estimators are unbiased. The variance is where a design earns its keep, and in the design-based world it is not an asymptotic approximation. It is an exact identity about a finite sum over $(\binom{n}{n_{1}})$ assignments, and it is the identity Neyman wrote down in 1923.

This chapter has one theorem in it. Everything else — Bernoulli, complete randomization, Horvitz–Thompson, Neyman — falls out of the variance of a linear statistic $\sum_{i} a_{i} Z_{i}$ in terms of the inclusion probabilities $π_{i}$ and $π_{i j}$ . I will spend most of the space on that formula and then cash it in.

Three population variances

In a design-based setup there is no superpopulation, so the quantities

S_{1}^{2} = \frac{1}{n - 1} \sum_{i} (y_{i} (1) - \bar{y} (1))^{2}, S_{0}^{2} = \frac{1}{n - 1} \sum_{i} (y_{i} (0) - \bar{y} (0))^{2}, S_{τ}^{2} = \frac{1}{n - 1} \sum_{i} (τ_{i} - τ)^{2}

are not parameters of a distribution. They are fixed numbers attached to the fixed population, exactly like $τ$ itself, and the only reason they carry an $n - 1$ is that it is the convention under which Neyman’s formula comes out without stray factors. Nothing here is a degrees-of-freedom argument.

The finite-population variance and covarianceNew tabOpens with 1,053 lines of library preamble — the code above is at the bottom of the editor.

/-- The **finite-population variance** of a vector `v`, with the `n - 1` denominator:
`fpVar v = (∑ i, (v i - v̄)²) / (n - 1)`.

This is a fixed number attached to the population, not an estimate.  Following the convention
of `MrCLean/Basic.lean`, Lean's `x / 0 = 0` makes `fpVar` equal `0` when `n ≤ 1`, so the
definition needs no side condition. -/
noncomputable def fpVar (v : Fin n → ℝ) : ℝ := (∑ i, (v i - mean v) ^ 2) / ((n : ℝ) - 1)

/-- The **finite-population covariance** of two vectors,
`fpCov u v = (∑ i, (u i - ū)(v i - v̄)) / (n - 1)`.

For the causal problem this is the quantity that can never be estimated: `fpCov y1 y0` needs
both potential outcomes of the same unit, and an experiment reveals only one of them. -/
noncomputable def fpCov (u v : Fin n → ℝ) : ℝ :=
  (∑ i, (u i - mean u) * (v i - mean v)) / ((n : ℝ) - 1)

Two things to read off the Lean rather than the mathematics.

First, fpVar takes a plain vector v : Fin n → ℝ, not a population and not a random variable. There is no probability anywhere in this definition, which is the whole point: fpVar P.y1 is as deterministic as P.y1.

Second, the definition has no hypothesis 1 < n, even though the displayed formula divides by $n - 1$ . That is Lean’s x / 0 = 0 convention doing its usual work: at $n = 1$ the denominator is $0$ and fpVar v = 0, and at $n = 0$ the numerator is an empty sum. The definition is total, and the hypothesis appears later only in the lemmas that actually need it.

The three named variances of the theorem, and the covariance that is the source of all the trouble:

S₁², S₀², S_τ² and S₁₀New tabOpens with 1,224 lines of library preamble — the code above is at the bottom of the editor.

/-- `S₁²`, the finite-population variance of the treated potential outcomes. -/
noncomputable def S1sq : ℝ := fpVar P.y1

/-- `S₀²`, the finite-population variance of the control potential outcomes. -/
noncomputable def S0sq : ℝ := fpVar P.y0

/-- `S_τ²`, the finite-population variance of the individual treatment effects.

This is the quantity that no experiment can estimate, and the one that makes the standard
variance estimator conservative. -/
noncomputable def Stausq : ℝ := fpVar P.effect

/-- `S₁₀`, the finite-population covariance of the two potential outcomes. -/
noncomputable def S10 : ℝ := fpCov P.y1 P.y0

$S_{10} = fpCov (y (1), y (0))$ is the quantity no experiment can estimate. It needs both potential outcomes of the same unit, and the fundamental problem of causal inference says you never see both. It is not that $S_{10}$ is hard to estimate; it is that nothing in the observed data moves when $S_{10}$ changes. And it is tied to $S_{τ}^{2}$ by a schoolbook identity:

S_τ² = S₁² + S₀² − 2 S₁₀New tabOpens with 1,239 lines of library preamble — the code above is at the bottom of the editor.

/-- **`S_τ² = S₁² + S₀² - 2 S₁₀`.**

The first two terms are estimable from an experiment; the third is not, because it needs both
potential outcomes of the same unit. -/
theorem Stausq_eq : P.Stausq = P.S1sq + P.S0sq - 2 * P.S10 := by
  have h : P.effect = fun i => P.y1 i - P.y0 i := rfl
  rw [Stausq, h, fpVar_sub]
  rfl

$S_{1}^{2}$ and $S_{0}^{2}$ are estimable, $S_{10}$ is not, and therefore $S_{τ}^{2}$ is not. Remember that sentence; it is the entire reason the last section of this chapter exists.

ExerciseScaling a vector (warm-up)

mean_const_mul pulls the c out of the mean; after that every summand is c ^ 2 * (v i - mean v) ^ 2, so the goal is Finset.mul_sum plus a ring on each term. Finset.sum_congr rfl fun i _ => … is the idiom for “prove this sum equals that sum by proving it termwise”.

New tabOpens with 1,084 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  Finite-population variance is quadratic: `fpVar (c·v) = c² fpVar v`.

Hint: `mean_const_mul` moves the constant out of the mean, and then every summand is
`c² (v i - v̄)²`. -/
theorem fpVar_const_mul (c : ℝ) (v : Fin n → ℝ) :
    fpVar (fun i => c * v i) = c ^ 2 * fpVar v := by
  sorry

Show solution

New tabOpens with 1,084 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  Finite-population variance is quadratic: `fpVar (c·v) = c² fpVar v`.

Hint: `mean_const_mul` moves the constant out of the mean, and then every summand is
`c² (v i - v̄)²`. -/
theorem fpVar_const_mul (c : ℝ) (v : Fin n → ℝ) :
    fpVar (fun i => c * v i) = c ^ 2 * fpVar v := by
  simp only [fpVar, mean_const_mul]
  rw [← mul_div_assoc]
  congr 1
  rw [Finset.mul_sum]
  exact Finset.sum_congr rfl fun i _ => by ring

ExerciseVariance of a difference (core)

Do not expand anything. Write u i - v i as 1 * u i + (-1) * v i with funext i; ring, then quote fpVar_linear_comb and let ring finish. This is the lemma that produces Stausq_eq above.

New tabOpens with 1,116 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (core).  The variance of a difference:
`fpVar (u - v) = fpVar u + fpVar v - 2 fpCov u v`.

Specialised to `u = y1`, `v = y0` this is the identity `S_τ² = S₁² + S₀² - 2 S₁₀` that makes
the un-estimable term in Neyman's formula visible.

Hint: `u i - v i = 1 * u i + (-1) * v i`, then quote `fpVar_linear_comb`. -/
theorem fpVar_sub (u v : Fin n → ℝ) :
    fpVar (fun i => u i - v i) = fpVar u + fpVar v - 2 * fpCov u v := by
  sorry

Show solution

New tabOpens with 1,116 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (core).  The variance of a difference:
`fpVar (u - v) = fpVar u + fpVar v - 2 fpCov u v`.

Specialised to `u = y1`, `v = y0` this is the identity `S_τ² = S₁² + S₀² - 2 S₁₀` that makes
the un-estimable term in Neyman's formula visible.

Hint: `u i - v i = 1 * u i + (-1) * v i`, then quote `fpVar_linear_comb`. -/
theorem fpVar_sub (u v : Fin n → ℝ) :
    fpVar (fun i => u i - v i) = fpVar u + fpVar v - 2 * fpCov u v := by
  have h : (fun i => u i - v i) = fun i => (1 : ℝ) * u i + (-1 : ℝ) * v i := by
    funext i; ring
  rw [h, fpVar_linear_comb]
  ring

ExerciseShifting a vector (warm-up)

mean_add and mean_const give mean (fun i => v i + c) = mean v + c; note that mean_const is where 0 < n is needed, so the hypothesis is not decoration. After rewriting the mean, each summand is (v i + c - (mean v + c)) ^ 2, which ring recognizes.

New tabOpens with 2,030 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  Shifting every entry of a vector by the same constant
leaves the finite-population variance alone.

Hint: `mean_add` and `mean_const` give `mean (v + c) = mean v + c`; after that
every summand `(v i + c - (v̄ + c))²` is `(v i - v̄)²`, which `ring` sees. -/
theorem fpVar_add_const {n : ℕ} (hn : 0 < n) (v : Fin n → ℝ) (c : ℝ) :
    fpVar (fun i => v i + c) = fpVar v := by
  sorry

Show solution

New tabOpens with 2,030 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  Shifting every entry of a vector by the same constant
leaves the finite-population variance alone.

Hint: `mean_add` and `mean_const` give `mean (v + c) = mean v + c`; after that
every summand `(v i + c - (v̄ + c))²` is `(v i - v̄)²`, which `ring` sees. -/
theorem fpVar_add_const {n : ℕ} (hn : 0 < n) (v : Fin n → ℝ) (c : ℝ) :
    fpVar (fun i => v i + c) = fpVar v := by
  have hm : mean (fun i => v i + c) = mean v + c := by
    rw [mean_add v (fun _ => c), mean_const hn c]
  simp only [fpVar, hm]
  congr 1
  exact Finset.sum_congr rfl fun i _ => by ring

The one formula: variance of a linear statistic

Every estimator in this library has the shape

\hat{θ} (z) = \sum_{i} a_{i} Z_{i} (z) + c

for coefficients $a$ and a constant $c$ that depend on the population and the design but not on the realized assignment. Difference in means has this shape (on the support of complete randomization), Horvitz–Thompson has it on every assignment, and the blocked and clustered estimators of chapter 10 have it too. So one variance formula does all the work.

The variance and covariance being computed are the design’s own, defined in chapter 7 as finite sums over assignments:

The design-based varianceNew tabOpens with 426 lines of library preamble — the code above is at the bottom of the editor.

/-- The **variance** of a statistic under the design, `Var f = E[(f - E f)²]`. -/
noncomputable def var : ℝ := D.expect (fun z => (f z - D.expect f) ^ 2)

D and f do not appear in the signature because they come from a variable line further up the file — Lean’s way of declaring parameters once for a run of declarations and letting each one pick up the ones it mentions. Design.cov is the same thing with two statistics. Both are expectations, so both are sums over all $2^{n}$ assignments weighted by D.prob; there is no measure theory to check.

The derivation is two lines of second-moment algebra. $Var$ of a sum is the double sum of covariances; constants pull out; and since $Z_{i} \in {0, 1}$ we have $Z_{i}^{2} = Z_{i}$ , so

Cov (Z_{i}, Z_{j}) = E [Z_{i} Z_{j}] - E [Z_{i}] E [Z_{j}] = π_{i j} - π_{i} π_{j} .

TheoremVariance of a linear statistic

For any design $D$ and any fixed coefficients $a_{1}, \dots, a_{n}$ ,

Var (\sum_{i} a_{i} Z_{i}) = \sum_{i} \sum_{j} a_{i} a_{j} (π_{i j} - π_{i} π_{j}) .

The Lean says exactly this, and the types are what make it precise. D.var : (Assignment n → ℝ) → ℝ is the design-based variance $E [(f - E f)^{2}]$ , a finite sum over assignments (defined in MrCLean/Design.lean). The statistic is fun z => ∑ i, a i * Z z i: z is the random object, a i is not — and you can see that in the type, because a : Fin n → ℝ has no z argument to take. A formalization in which the coefficients could depend on the assignment would be a different, and false, theorem.

New tabOpens with 1,367 lines of library preamble — the code above is at the bottom of the editor.

theorem var_linear (D : Design n) (a : Fin n → ℝ) :
    D.var (fun z => ∑ i, a i * Z z i)
      = ∑ i, ∑ j, a i * a j * (D.jointPropensity i j - D.propensity i * D.propensity j) := by
  rw [var_sum]
  refine Finset.sum_congr rfl fun i _ => Finset.sum_congr rfl fun j _ => ?_
  rw [cov_const_mul_left, cov_const_mul_right, cov_Z]
  ring

The proof is the four tactic lines shown beside the statement, and it is worth walking through because it is the template for every variance computation in the library.

rw [var_sum] replaces $Var (\sum_{i} F_{i})$ by $\sum_{i} \sum_{j} Cov (F_{i}, F_{j})$ — this is where linearity of expectation is spent, and it is proved once for arbitrary index types. Finset.sum_congr rfl fun i _ => … is the workhorse for sums: to show two sums over the same index set (that is the rfl) are equal, prove the summands equal one index at a time; the _ is the unused proof that i is in the set. The remaining goal is about a single pair $(i, j)$ , so cov_const_mul_left and cov_const_mul_right pull a i and a j out of the covariance, cov_Z names what is left, and ring rearranges. No design has been mentioned yet, which is the point.

The textbook display splits the diagonal off:

Diagonal and off-diagonalNew tabOpens with 1,383 lines of library preamble — the code above is at the bottom of the editor.

/-- The same formula with the diagonal split off, which is the textbook display:

`Var(∑ aᵢ Zᵢ) = ∑ᵢ aᵢ² πᵢ(1 - πᵢ) + ∑ᵢ ∑_{j ≠ i} aᵢ aⱼ (πᵢⱼ - πᵢ πⱼ)`.

The first sum is what one would get if the units were independent; the second is the price of
whatever dependence the design creates. -/
theorem var_linear_split (D : Design n) (a : Fin n → ℝ) :
    D.var (fun z => ∑ i, a i * Z z i)
      = (∑ i, (a i) ^ 2 * (D.propensity i * (1 - D.propensity i)))
        + ∑ i, ∑ j ∈ univ.erase i,
            a i * a j * (D.jointPropensity i j - D.propensity i * D.propensity j) := by
  rw [var_linear, ← Finset.sum_add_distrib]
  refine Finset.sum_congr rfl fun i _ => ?_
  rw [← Finset.add_sum_erase _ _ (Finset.mem_univ i)]
  congr 1
  rw [jointPropensity_self]
  ring

The first sum, $\sum_{i} a_{i}^{2} π_{i} (1 - π_{i})$ , is what you would get if the units were assigned independently. The second is the price of whatever dependence the design imposes. Note univ.erase i — Lean’s way of writing “all $j \neq i$ ” as a Finset — rather than a filter on a proposition; the two are equal, but the erase form has the cardinality lemmas attached.

Exerciseπᵢᵢ = πᵢ (warm-up)

This is Z_mul_self : Z z i * Z z i = Z z i and nothing else. Unfold jointPropensity to an expectation and use Design.expect_congr, which says two statistics that agree pointwise have the same expectation.

New tabOpens with 1,352 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  The joint propensity of a unit with itself is its propensity,
because `Z² = Z`. -/
theorem jointPropensity_self (D : Design n) (i : Fin n) :
    D.jointPropensity i i = D.propensity i := by
  sorry

Show solution

New tabOpens with 1,352 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  The joint propensity of a unit with itself is its propensity,
because `Z² = Z`. -/
theorem jointPropensity_self (D : Design n) (i : Fin n) :
    D.jointPropensity i i = D.propensity i := by
  rw [jointPropensity]
  exact D.expect_congr fun z => Z_mul_self z i

ExerciseAdding a constant (warm-up)

Compute the new expectation with expect_add and expect_const first, as a have. Then simp only [var, hE] puts the goal in the form E[(f z + c - (E f + c))²] = E[(f z - E f)²], and expect_congr reduces it to an identity ring can close. This lemma is what lets Neyman’s proof throw away the - (∑ y0)/n₀ tail of difference in means.

New tabOpens with 1,308 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  Adding a constant to a statistic does not change its variance.

This is what lets us ignore the `- (∑ y0)/n₀` tail of difference-in-means. -/
theorem var_add_const (D : Design n) (f : Assignment n → ℝ) (c : ℝ) :
    D.var (fun z => f z + c) = D.var f := by
  sorry

Show solution

New tabOpens with 1,308 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  Adding a constant to a statistic does not change its variance.

This is what lets us ignore the `- (∑ y0)/n₀` tail of difference-in-means. -/
theorem var_add_const (D : Design n) (f : Assignment n → ℝ) (c : ℝ) :
    D.var (fun z => f z + c) = D.var f := by
  have hE : D.expect (fun z => f z + c) = D.expect f + c := by
    rw [D.expect_add f (fun _ => c), D.expect_const c]
  simp only [var, hE]
  exact D.expect_congr fun z => by ring

ExerciseScaling a statistic (warm-up)

Three rewrites and a ring. Go through var_eq_cov_self so that the two cov_const_mul_* lemmas apply, then come back.

New tabOpens with 1,318 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  Scaling a statistic scales its variance by the square:
`Var(c · f) = c² Var f`. -/
theorem var_const_mul (D : Design n) (c : ℝ) (f : Assignment n → ℝ) :
    D.var (fun z => c * f z) = c ^ 2 * D.var f := by
  sorry

Show solution

New tabOpens with 1,318 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  Scaling a statistic scales its variance by the square:
`Var(c · f) = c² Var f`. -/
theorem var_const_mul (D : Design n) (c : ℝ) (f : Assignment n → ℝ) :
    D.var (fun z => c * f z) = c ^ 2 * D.var f := by
  rw [var_eq_cov_self, cov_const_mul_left, cov_const_mul_right, ← var_eq_cov_self]
  ring

Exchangeable designs, and the two we have

Both designs in the library are exchangeable: every unit has the same propensity $p$ , and every ordered pair of distinct units the same joint propensity $q$ . Then the double sum collapses, because the kernel takes only two values:

Exchangeable designsNew tabOpens with 1,440 lines of library preamble — the code above is at the bottom of the editor.

/-- **The variance of a linear statistic under an exchangeable design.**

If every propensity equals `p` and every distinct pair has joint propensity `q`, then

`Var(∑ aᵢ Zᵢ) = (p - p²) ∑ aᵢ² + (q - p²) ((∑ aᵢ)² - ∑ aᵢ²)`.

For Bernoulli, `q = p²` and the second term vanishes — that is exactly independence.  For
complete randomisation `q < p²`, and the negative second term is the finite-population
correction. -/
theorem var_linear_of_exchangeable (D : Design n) (a : Fin n → ℝ) (p q : ℝ)
    (hp : ∀ i, D.propensity i = p) (hq : ∀ i j, i ≠ j → D.jointPropensity i j = q) :
    D.var (fun z => ∑ i, a i * Z z i)
      = (p - p ^ 2) * (∑ i, (a i) ^ 2) + (q - p ^ 2) * ((∑ i, a i) ^ 2 - ∑ i, (a i) ^ 2) := by
  rw [var_linear]
  exact sum_sum_diag_offdiag a (p - p ^ 2) (q - p ^ 2)
    (fun i j => D.jointPropensity i j - D.propensity i * D.propensity j)
    (fun i => by rw [jointPropensity_self, hp i]; ring)
    (fun i j hij => by rw [hq i j hij, hp i, hp j]; ring)

The combinatorial content is $\sum_{i} \sum_{j \neq i} a_{i} a_{j} = (\sum_{i} a_{i})^{2} - \sum_{i} a_{i}^{2}$ , isolated as sum_sum_diag_offdiag. Specializing is now mechanical.

Bernoulli: independenceNew tabOpens with 1,463 lines of library preamble — the code above is at the bottom of the editor.

/-- **Under the Bernoulli design the units are independent**, so the variance of a linear
statistic is the plain sum of the individual variances:

`Var(∑ aᵢ Zᵢ) = p(1 - p) ∑ aᵢ²`. -/
theorem bernoulli_var_linear {p : ℝ} (hp0 : 0 ≤ p) (hp1 : p ≤ 1) (a : Fin n → ℝ) :
    (bernoulliDesign n p hp0 hp1).var (fun z => ∑ i, a i * Z z i)
      = p * (1 - p) * ∑ i, (a i) ^ 2 := by
  rw [Design.var_linear_of_exchangeable _ a p (p ^ 2)
    (fun i => bernoulli_propensity i) (fun i j hij => bernoulli_jointPropensity hij)]
  ring

For Bernoulli $q = p^{2}$ , so the off-diagonal term is exactly zero and the variance is the plain sum $p (1 - p) \sum_{i} a_{i}^{2}$ . That vanishing is what “independent assignment” means here; it is not assumed, it is computed from bernoulli_jointPropensity.

Complete randomizationNew tabOpens with 1,474 lines of library preamble — the code above is at the bottom of the editor.

/-- **The variance of a linear statistic under complete randomisation.**

`Var(∑ aᵢ Zᵢ) = (n₁ n₀ / n) · fpVar a`, where `n₀ = n - n₁`.

The whole finite-population variance `fpVar a` appears because fixing the number of treated
units makes the indicators negatively dependent by exactly the amount that turns
`∑ aᵢ²` into the *centred* sum of squares. -/
theorem completeRandomization_var_linear {n₁ : ℕ} (hle : n₁ ≤ n) (hpos : 0 < n₁) (hlt : n₁ < n)
    (a : Fin n → ℝ) :
    (completeRandomization n n₁ hle).var (fun z => ∑ i, a i * Z z i)
      = ((n₁ : ℝ) * ((n : ℝ) - (n₁ : ℝ)) / (n : ℝ)) * fpVar a := by
  have hn2 : 2 ≤ n := by omega
  have hn : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr (by omega)
  have hn1 : (n : ℝ) - 1 ≠ 0 := by
    have : (2 : ℝ) ≤ (n : ℝ) := by exact_mod_cast hn2
    linarith
  rw [Design.var_linear_of_exchangeable _ a ((n₁ : ℝ) / (n : ℝ))
      ((n₁ : ℝ) * ((n₁ : ℝ) - 1) / ((n : ℝ) * ((n : ℝ) - 1)))
      (fun i => completeRandomization_propensity i)
      (fun i j hij => completeRandomization_jointPropensity hij),
    fpVar_eq_sum_sq hn]
  field_simp
  ring

For complete randomization $p = n_{1} / n$ and $q = n_{1} (n_{1} - 1) / (n (n - 1))$ , and

q - p^{2} = \frac{n_{1} (n_{1} - 1)}{n (n - 1)} - \frac{n_{1}^{2}}{n^{2}} < 0 whenever 0 < n_{1} < n,

so the off-diagonal term is negative: fixing the number of treated units makes the indicators negatively dependent. And the negative term is exactly the right size to convert $\sum_{i} a_{i}^{2}$ into the centered sum of squares, which is why the whole thing collapses to

Var (\sum_{i} a_{i} Z_{i}) = \frac{n_{1} n_{0}}{n} fpVar (a) .

The field_simp; ring at the end of that proof is doing the cancellation; the mathematical content is entirely in the two propensity lemmas from chapter 7.

Neyman’s theorem

TheoremNeyman (1923)

Under complete randomization with $n_{1}$ treated and $n_{0} = n - n_{1}$ control units, $0 < n_{1} < n$ ,

Var ({\hat{τ}}_{DiM}) = \frac{S_{1}^{2}}{n_{1}} + \frac{S_{0}^{2}}{n_{0}} - \frac{S_{τ}^{2}}{n} .

An exact identity for finite $n$ , not a limit. The third term is subtracted: the randomization makes the treated and control averages negatively dependent, so difference in means is less variable than two independent samples of sizes $n_{1}$ and $n_{0}$ would be, and less by exactly the variance of the individual treatment effects.

New tabOpens with 1,549 lines of library preamble — the code above is at the bottom of the editor.

theorem neyman_variance {n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n) (hpos : 0 < n₁)
    (hlt : n₁ < n) :
    (completeRandomization n n₁ hle).var (DiM P)
      = P.S1sq / (n₁ : ℝ) + P.S0sq / ((n : ℝ) - (n₁ : ℝ)) - P.Stausq / (n : ℝ) := by
  have hn : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr (by omega)
  have hn1 : (n₁ : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr (by omega)
  have hn0 : (n : ℝ) - (n₁ : ℝ) ≠ 0 := by
    have : (n₁ : ℝ) < (n : ℝ) := by exact_mod_cast hlt
    linarith
  -- Step 1: difference-in-means is a linear statistic plus a constant, on the support.
  have hlin : (completeRandomization n n₁ hle).var (DiM P)
      = (completeRandomization n n₁ hle).var
          (fun z => ∑ i, dimCoef P (n₁ : ℝ) ((n : ℝ) - (n₁ : ℝ)) i * Z z i) := by
    rw [(completeRandomization n n₁ hle).var_congr_of_support
      (g := fun z => (∑ i, dimCoef P (n₁ : ℝ) ((n : ℝ) - (n₁ : ℝ)) i * Z z i)
        + -((∑ i, P.y0 i) / ((n : ℝ) - (n₁ : ℝ))))
      (fun z hz => DiM_eq_linear_of_card P hlt (card_eq_of_prob_ne_zero hz))]
    exact (completeRandomization n n₁ hle).var_add_const _ _
  -- Step 2: the general linear-statistic variance, then pure algebra.
  rw [hlin, completeRandomization_var_linear hle hpos hlt, fpVar_dimCoef,
    P.Stausq_eq]
  field_simp
  ring

Read the hypotheses. hle : n₁ ≤ n is not a mathematical assumption at all — it is an argument to the design: completeRandomization n n₁ hle cannot be formed without it, because there are no subsets of size $n_{1} > n$ to be uniform over. The mathematical assumptions are hpos : 0 < n₁ and hlt : n₁ < n, and they are there because both groups must be non-empty for the two sample averages to mean anything. Drop hpos and the statement becomes false rather than undefined: with $n_{1} = 0$ the design puts all its mass on the empty assignment, difference in means is the constant $- \bar{y} (0)$ and has variance $0$ , while the right-hand side reads $0 + S_{0}^{2} / n - S_{τ}^{2} / n$ because S1sq / 0 = 0.

Read the arithmetic, too. The control group size appears as (n : ℝ) - (n₁ : ℝ), real subtraction, not ((n - n₁ : ℕ) : ℝ). Natural subtraction truncates, and mixing the two is the single most common way to get a formalized finite-population statement subtly wrong. neyman_variance_nat is the same theorem with a natural n₀ and a hypothesis n₁ + n₀ = n; it does the cast once, correctly, so that downstream users do not have to.

The proof, in four steps

Step 1 — make difference in means linear. On assignments with $# z = n_{1}$ ,

{\hat{τ}}_{DiM} (z) = \sum_{i} a_{i} Z_{i} (z) - \frac{1}{n_{0}} \sum_{i} y_{i} (0), a_{i} = \frac{y_{i} (1)}{n_{1}} + \frac{y_{i} (0)}{n_{0}} .

The coefficient vectorNew tabOpens with 1,505 lines of library preamble — the code above is at the bottom of the editor.

/-- The coefficient vector that makes difference-in-means a linear statistic:
`aᵢ = y1ᵢ/n₁ + y0ᵢ/n₀`.

It is worth staring at: the *control* outcomes enter with a **plus** sign, because a unit that
is treated is simultaneously removed from the control average. -/
noncomputable def dimCoef (P : Population n) (n₁ n₀ : ℝ) (i : Fin n) : ℝ :=
  (1 / n₁) * P.y1 i + (1 / n₀) * P.y0 i

/-- **On the support of complete randomisation, difference-in-means is linear in the
indicators**: `DiM P z = ∑ i, aᵢ Zᵢ + c` with `a = dimCoef` and a constant `c` that does not
depend on the assignment. -/
theorem DiM_eq_linear_of_card {n₁ : ℕ} (P : Population n) (hlt : n₁ < n)
    {z : Assignment n} (hz : #z = n₁) :
    DiM P z
      = (∑ i, dimCoef P (n₁ : ℝ) ((n : ℝ) - (n₁ : ℝ)) i * Z z i)
        + -((∑ i, P.y0 i) / ((n : ℝ) - (n₁ : ℝ))) := by
  have hcompl : ((#zᶜ : ℕ) : ℝ) = (n : ℝ) - (n₁ : ℝ) := by
    rw [card_compl, hz, Nat.cast_sub hlt.le]
  calc DiM P z
      = (∑ i, Z z i * P.y1 i) / (n₁ : ℝ)
          - (∑ i, (1 - Z z i) * P.y0 i) / ((n : ℝ) - (n₁ : ℝ)) := by
        unfold DiM
        rw [sum_mem_eq_sum_Z_mul, sum_compl_eq_sum_one_sub_Z_mul, hz, hcompl]
    _ = ∑ i, (Z z i * P.y1 i / (n₁ : ℝ)
          - (1 - Z z i) * P.y0 i / ((n : ℝ) - (n₁ : ℝ))) := by
        rw [Finset.sum_div, Finset.sum_div, Finset.sum_sub_distrib]
    _ = ∑ i, (dimCoef P (n₁ : ℝ) ((n : ℝ) - (n₁ : ℝ)) i * Z z i
          + -(P.y0 i / ((n : ℝ) - (n₁ : ℝ)))) :=
        Finset.sum_congr rfl fun i _ => by unfold dimCoef; ring
    _ = (∑ i, dimCoef P (n₁ : ℝ) ((n : ℝ) - (n₁ : ℝ)) i * Z z i)
          + ∑ i, -(P.y0 i / ((n : ℝ) - (n₁ : ℝ))) := Finset.sum_add_distrib
    _ = (∑ i, dimCoef P (n₁ : ℝ) ((n : ℝ) - (n₁ : ℝ)) i * Z z i)
          + -((∑ i, P.y0 i) / ((n : ℝ) - (n₁ : ℝ))) := by
        rw [Finset.sum_neg_distrib, Finset.sum_div]

The sign is worth staring at: the control outcome enters $a_{i}$ with a plus. Treating unit $i$ both adds $y_{i} (1) / n_{1}$ to the treated average and removes $y_{i} (0) / n_{0}$ from the control average, and the estimator is a difference, so the two effects add.

Step 2 — get onto the support. The identity above holds only when $# z = n_{1}$ , but Design.var sums over all $2^{n}$ assignments, including the ones with probability zero on which DiM divides by zero and returns garbage. This is what Design.var_congr_of_support is for: two statistics that agree wherever D.prob z ≠ 0 have the same variance. card_eq_of_prob_ne_zero supplies #z = n₁ from D.prob z ≠ 0. Then var_add_const deletes the constant tail.

Step 3 — quote the general formula. completeRandomization_var_linear turns the variance into $\frac{n_{1} n_{0}}{n} fpVar (a)$ .

Step 4 — expand. fpVar_dimCoef expands $fpVar (a)$ into $S_{1}^{2} / n_{1}^{2} + S_{0}^{2} / n_{0}^{2} + 2 S_{10} / (n_{1} n_{0})$ , Stausq_eq replaces $S_{10}$ by $(S_{1}^{2} + S_{0}^{2} - S_{τ}^{2}) / 2$ , and field_simp; ring does the algebra with the three non-vanishing denominators established as haves at the top of the proof.

A worked example you can check by hand

Four units, two treated. Two of the four units gain $2$ from treatment and two are unaffected, so $τ = 1$ and the effect is heterogeneous.

A four-unit populationNew tabOpens with 1,936 lines of library preamble — the code above is at the bottom of the editor.

/-- A concrete population of `n = 4` units, written with Mathlib's vector
notation `![a, b, c, d] : Fin 4 → ℝ`.

Unit 0 is unaffected, units 1 and 2 gain 2, unit 3 is unaffected: the individual
effects are `(0, 2, 2, 0)`, so `τ = 1` and the effect is *not* constant. -/
noncomputable def exP : Population 4 where
  y1 := ![0, 2, 2, 2]
  y0 := ![0, 0, 0, 2]

/-- The average treatment effect of `exP` is `1`. -/
theorem exP_tau : exP.tau = 1 := by
  norm_num [Population.tau, Population.effect, exP, mean, Fin.sum_univ_four, Matrix.cons_val_two, Matrix.cons_val_three,
      Matrix.head_cons, Matrix.tail_cons]

ExerciseThe three variances of the example (warm-up)

$\bar{y} (1) = 3 / 2$ and $\bar{y} (0) = 1 / 2$ , so both centered sums of squares are $3$ and $S_{1}^{2} = S_{0}^{2} = 3 / 3 = 1$ . The effect vector is $(0, 2, 2, 0)$ with mean $1$ , centered sum of squares $4$ , so $S_{τ}^{2} = 4 / 3$ . In Lean: everything in sight is a definition, so put Population.S1sq, fpVar, mean and exP in the simp set, expand the sum with Fin.sum_univ_four, and let norm_num finish. ![a, b, c, d] 2 needs Matrix.cons_val_two to reduce — that one is not a simp lemma.

New tabOpens with 1,950 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  The three finite-population variances of `exP`.

`S₁² = S₀² = 1` and `S_τ² = 4/3`: the effect vector `(0, 2, 2, 0)` is more spread
out around its mean than either potential-outcome vector is around its own.

Hint: unfold to a sum over `Fin 4` with `Fin.sum_univ_four` and let `norm_num`
do the arithmetic.  `Population.S1sq` and friends are definitions, so they have
to appear in the simp set. -/
theorem exP_variances : exP.S1sq = 1 ∧ exP.S0sq = 1 ∧ exP.Stausq = 4 / 3 := by
  sorry

Show solution

New tabOpens with 1,950 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (warm-up).  The three finite-population variances of `exP`.

`S₁² = S₀² = 1` and `S_τ² = 4/3`: the effect vector `(0, 2, 2, 0)` is more spread
out around its mean than either potential-outcome vector is around its own.

Hint: unfold to a sum over `Fin 4` with `Fin.sum_univ_four` and let `norm_num`
do the arithmetic.  `Population.S1sq` and friends are definitions, so they have
to appear in the simp set. -/
theorem exP_variances : exP.S1sq = 1 ∧ exP.S0sq = 1 ∧ exP.Stausq = 4 / 3 := by
  refine ⟨?_, ?_, ?_⟩ <;>
    norm_num [Population.S1sq, Population.S0sq, Population.Stausq, Population.effect,
      fpVar, mean, exP, Fin.sum_univ_four, Matrix.cons_val_two, Matrix.cons_val_three,
      Matrix.head_cons, Matrix.tail_cons]

ExerciseNeyman's formula on the example (core)

$\frac{1}{2} + \frac{1}{2} - \frac{4 / 3}{4} = \frac{2}{3}$ . Quote neyman_variance with the four arguments it wants (all four discharged by norm_num), rewrite with the three variances from the previous exercise, and finish with norm_num. Note that the (by norm_num : 2 ≤ 4) in the goal and the one you pass are different terms but the same proof — proofs of a Prop are interchangeable, so rw does not care.

New tabOpens with 1,995 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (core).  **Neyman's formula on the worked example.**

With `n = 4`, `n₁ = n₀ = 2`,

`Var(DiM) = S₁²/2 + S₀²/2 - S_τ²/4 = 1/2 + 1/2 - 1/3 = 2/3`.

Hint: `neyman_variance` reduces the goal to arithmetic in the three population
variances; `exP_variances` supplies those. -/
theorem exP_var_DiM :
    (completeRandomization 4 2 (by norm_num)).var (DiM exP) = 2 / 3 := by
  sorry

Show solution

New tabOpens with 1,995 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (core).  **Neyman's formula on the worked example.**

With `n = 4`, `n₁ = n₀ = 2`,

`Var(DiM) = S₁²/2 + S₀²/2 - S_τ²/4 = 1/2 + 1/2 - 1/3 = 2/3`.

Hint: `neyman_variance` reduces the goal to arithmetic in the three population
variances; `exP_variances` supplies those. -/
theorem exP_var_DiM :
    (completeRandomization 4 2 (by norm_num)).var (DiM exP) = 2 / 3 := by
  obtain ⟨h1, h0, hτ⟩ := exP_variances
  rw [neyman_variance (by norm_num) exP (by norm_num) (by norm_num), h1, h0, hτ]
  norm_num

There are $(\binom{4}{2}) = 6$ equally likely assignments, and the estimator takes a different value on each. Here is one of them:

One draw from the designNew tabOpens with 1,964 lines of library preamble — the code above is at the bottom of the editor.

/-- One draw from the design.  If the coin puts units 2 and 3 in treatment, the
difference in means is `2` — twice the true `τ = 1`.  Nothing is wrong: a single
realisation of an unbiased estimator has no obligation to be close. -/
theorem exP_DiM_of_treated_two_three : DiM exP {2, 3} = 2 := by
  have hcompl : ({2, 3} : Assignment 4)ᶜ = {0, 1} := by decide
  unfold DiM
  rw [hcompl, Finset.sum_pair (by decide : (2 : Fin 4) ≠ 3),
    Finset.sum_pair (by decide : (0 : Fin 4) ≠ 1),
    Finset.card_pair (by decide : (2 : Fin 4) ≠ 3),
    Finset.card_pair (by decide : (0 : Fin 4) ≠ 1)]
  norm_num [exP, Matrix.cons_val_two, Matrix.cons_val_three, Matrix.head_cons,
    Matrix.tail_cons]

/-- On that same draw the Neyman variance estimator reports **zero**: both
treated units have `y1 = 2` and both control units have `y0 = 0`, so both sample
variances vanish.

So the estimator is *not* conservative draw by draw — `0 < 2/3 = Var(DiM)`.
Conservativeness is a statement about its expectation and nothing else. -/
theorem exP_neymanVarEst_of_treated_two_three : neymanVarEst exP {2, 3} = 0 := by
  have hcompl : ({2, 3} : Assignment 4)ᶜ = {0, 1} := by decide
  unfold neymanVarEst sampleVar
  rw [hcompl, Finset.sum_pair (by decide : (2 : Fin 4) ≠ 3),
    Finset.sum_pair (by decide : (0 : Fin 4) ≠ 1),
    Finset.sum_pair (by decide : (2 : Fin 4) ≠ 3),
    Finset.sum_pair (by decide : (0 : Fin 4) ≠ 1),
    Finset.card_pair (by decide : (2 : Fin 4) ≠ 3),
    Finset.card_pair (by decide : (0 : Fin 4) ≠ 1)]
  norm_num [exP, Matrix.cons_val_two, Matrix.cons_val_three, Matrix.head_cons,
    Matrix.tail_cons]

If the randomization puts units 2 and 3 in treatment, both treated units have $y (1) = 2$ and both control units have $y (0) = 0$ , so the point estimate is $2$ — twice the true $τ$ — and both sample variances, and hence the reported variance, are $0$ . Hold on to that pair of numbers; the last section needs them.

The variance estimator, and what “conservative” means

Nobody reports $S_{1}^{2} / n_{1} + S_{0}^{2} / n_{0} - S_{τ}^{2} / n$ , because nobody knows $S_{τ}^{2}$ . What gets reported is the plug-in estimator built from the two sample variances.

The Neyman variance estimatorNew tabOpens with 1,602 lines of library preamble — the code above is at the bottom of the editor.

/-- The **sample variance** of a vector `v` over a group `s` of units, with the usual `#s - 1`
denominator.  Applied to the treated set it is `s₁²`; applied to the control set, `s₀²`. -/
noncomputable def sampleVar (v : Fin n → ℝ) (s : Finset (Fin n)) : ℝ :=
  (∑ i ∈ s, (v i - (∑ j ∈ s, v j) / (#s : ℝ)) ^ 2) / ((#s : ℝ) - 1)

/-- The **Neyman variance estimator**: the treated sample variance over `n₁` plus the control
sample variance over `n₀`.

It is a genuine statistic — it uses only the outcomes the assignment reveals, since `y1` is
summed over the treated set and `y0` over the control set. -/
noncomputable def neymanVarEst (P : Population n) (z : Assignment n) : ℝ :=
  sampleVar P.y1 z / (#z : ℝ) + sampleVar P.y0 zᶜ / ((#zᶜ : ℕ) : ℝ)

This is a genuine statistic, and you can check that by reading the definition rather than trusting the name: P.y1 is summed only over z, and P.y0 only over zᶜ. Every outcome it touches is one the assignment reveals.

Why the sample variances are unbiased

The engine is a single design-agnostic lemma, and it covers both groups at once by taking the group indicator as a parameter — W = Z for the treated group, W = fun z i => 1 - Z z i for the control group.

Sample variance of a simple random sample without replacementNew tabOpens with 1,622 lines of library preamble — the code above is at the bottom of the editor.

/-- **The unbiasedness engine for sample variances.**

Let `W z i ∈ {0, 1}` be the indicator of "unit `i` is in the group" under assignment `z`
(so `W = Z` for the treated group and `W = 1 - Z` for the control group).  If the group has
fixed size `m`, every unit is included with probability `m/n`, and every pair with probability
`m(m-1)/(n(n-1))` — the moments of simple random sampling without replacement — then

`E[ ∑ᵢ Wᵢ vᵢ² - (∑ᵢ Wᵢ vᵢ)²/m ] = (m - 1) · fpVar v`.

Dividing by `m - 1` says the sample variance of a simple random sample without replacement is
unbiased for the finite-population variance.  The `n - 1` in `fpVar` is exactly what makes
this come out clean. -/
theorem expect_ss_of_exchangeable (D : Design n) (v : Fin n → ℝ)
    (W : Assignment n → Fin n → ℝ) (m : ℝ)
    (hn : (n : ℝ) ≠ 0) (hn1 : (n : ℝ) - 1 ≠ 0) (hm : m ≠ 0)
    (hidem : ∀ z i, W z i * W z i = W z i)
    (h1 : ∀ i, D.expect (fun z => W z i) = m / (n : ℝ))
    (h2 : ∀ i j, i ≠ j →
      D.expect (fun z => W z i * W z j) = m * (m - 1) / ((n : ℝ) * ((n : ℝ) - 1))) :
    D.expect (fun z => (∑ i, W z i * (v i) ^ 2) - (∑ i, W z i * v i) ^ 2 / m)
      = (m - 1) * fpVar v := by
  -- first moment: the expected sum of squares over the group
  have hA : D.expect (fun z => ∑ i, W z i * (v i) ^ 2) = (m / (n : ℝ)) * ∑ i, (v i) ^ 2 := by
    rw [D.expect_sum univ (fun i z => W z i * (v i) ^ 2), Finset.mul_sum]
    exact Finset.sum_congr rfl fun i _ => by
      rw [D.expect_mul_const (fun z => W z i) ((v i) ^ 2), h1 i]
  -- second moment: the expected squared group total
  have hB : D.expect (fun z => (∑ i, W z i * v i) ^ 2)
      = (m / (n : ℝ)) * (∑ i, (v i) ^ 2)
        + (m * (m - 1) / ((n : ℝ) * ((n : ℝ) - 1))) * ((∑ i, v i) ^ 2 - ∑ i, (v i) ^ 2) := by
    have hsq : ∀ z : Assignment n,
        (∑ i, W z i * v i) ^ 2 = ∑ i, ∑ j, v i * v j * (W z i * W z j) := by
      intro z
      rw [sq, Finset.sum_mul_sum]
      exact Finset.sum_congr rfl fun i _ =>
        Finset.sum_congr rfl fun j _ => by ring
    rw [D.expect_congr hsq,
      D.expect_sum univ (fun i z => ∑ j, v i * v j * (W z i * W z j)),
      Finset.sum_congr rfl fun i _ =>
        D.expect_sum univ (fun j z => v i * v j * (W z i * W z j)),
      Finset.sum_congr rfl fun i _ => Finset.sum_congr rfl fun j _ =>
        D.expect_const_mul (fun z => W z i * W z j) (v i * v j)]
    exact sum_sum_diag_offdiag v (m / (n : ℝ)) (m * (m - 1) / ((n : ℝ) * ((n : ℝ) - 1)))
      (fun i j => D.expect (fun z => W z i * W z j))
      (fun i => by rw [D.expect_congr (fun z => hidem z i), h1 i])
      (fun i j hij => h2 i j hij)
  rw [D.expect_sub _ _, hA, D.expect_div_const _ m, hB, fpVar_eq_sum_sq hn]
  field_simp
  ring

In words: the sample variance of a simple random sample without replacement is unbiased for the finite-population variance, which is the classical reason the $n - 1$ is where it is. The hypotheses are exactly the first two moments of sampling without replacement — h1 says each unit is in the group with probability $m / n$ , h2 says each pair is with probability $m (m - 1) / (n (n - 1))$ — plus hidem, which says $W$ takes values in ${0, 1}$ by asserting $W^{2} = W$ . Stating idempotence instead of " $W z i = 0 \lor W z i = 1$ " is what keeps the proof inside ring.

Feeding W = Z gives $E [s_{1}^{2}] = S_{1}^{2}$ , feeding W = 1 - Z gives $E [s_{0}^{2}] = S_{0}^{2}$ — the control group of a completely randomized experiment is itself a simple random sample of size $n_{0}$ — and adding them:

The expectation of the estimatorNew tabOpens with 1,772 lines of library preamble — the code above is at the bottom of the editor.

/-- **The expectation of the Neyman variance estimator is `S₁²/n₁ + S₀²/n₀`.**

Notice what is missing: the `- S_τ²/n` of Neyman's formula.  The estimator has no way to
recover it, because `S_τ²` involves `fpCov y1 y0` and no unit ever shows both outcomes. -/
theorem completeRandomization_expect_neymanVarEst {n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n)
    (hn₁ : 1 < n₁) (hn₀ : 1 < n - n₁) :
    (completeRandomization n n₁ hle).expect (neymanVarEst P)
      = P.S1sq / (n₁ : ℝ) + P.S0sq / ((n : ℝ) - (n₁ : ℝ)) := by
  set D := completeRandomization n n₁ hle with hD
  have hlt : n₁ < n := by omega
  have hsupp : ∀ z, D.prob z ≠ 0 →
      neymanVarEst P z
        = sampleVar P.y1 z / (n₁ : ℝ) + sampleVar P.y0 zᶜ / ((n : ℝ) - (n₁ : ℝ)) := by
    intro z hz
    have hcard : #z = n₁ := card_eq_of_prob_ne_zero hz
    have hc : ((#zᶜ : ℕ) : ℝ) = (n : ℝ) - (n₁ : ℝ) := by
      rw [card_compl, hcard, Nat.cast_sub hle]
    rw [neymanVarEst, hcard, hc]
  rw [D.expect_congr_of_support hsupp,
    D.expect_add (fun z => sampleVar P.y1 z / (n₁ : ℝ))
      (fun z => sampleVar P.y0 zᶜ / ((n : ℝ) - (n₁ : ℝ))),
    D.expect_div_const _ (n₁ : ℝ), D.expect_div_const _ ((n : ℝ) - (n₁ : ℝ)), hD,
    completeRandomization_expect_sampleVar_treated hle P hn₁ hlt,
    completeRandomization_expect_sampleVar_control hle P (by omega) hn₀]

Notice what is missing from the right-hand side: the $- S_{τ}^{2} / n$ . There is no way to recover it. Subtract Neyman’s formula and the bias is not merely bounded but identified:

The bias, exactlyNew tabOpens with 1,797 lines of library preamble — the code above is at the bottom of the editor.

/-- **The bias of the Neyman variance estimator is exactly `S_τ²/n`.**

`E[v̂] = Var(DiM) + S_τ²/n`.  This is the sharpest possible statement of the classical
"conservativeness": we do not merely bound the bias, we identify it. -/
theorem neymanVarEst_bias {n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n)
    (hn₁ : 1 < n₁) (hn₀ : 1 < n - n₁) :
    (completeRandomization n n₁ hle).expect (neymanVarEst P)
      = (completeRandomization n n₁ hle).var (DiM P) + P.Stausq / (n : ℝ) := by
  rw [completeRandomization_expect_neymanVarEst hle P hn₁ hn₀,
    neyman_variance hle P (by omega) (by omega)]
  ring

Note the hypotheses have changed. Neyman’s theorem needed only $0 < n_{1} < n$ ; everything about the estimator needs 1 < n₁ and 1 < n - n₁, because both sample variances divide by $n_{i} - 1$ . This is the sort of difference that is invisible in a paper and glaring in a formalization.

Conservative, as an inequality between two numbers

ExerciseThe estimator is conservative (core)

One rewrite and one linarith: neymanVarEst_bias turns the goal into Var ≤ Var + Stausq / n, and fpVar_nonneg says the added term is a non-negative real divided by a non-negative real. Note that fpVar_nonneg needs no hypothesis, so no side condition leaks into this theorem beyond the two the estimator already required.

New tabOpens with 1,809 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (core).  **The Neyman variance estimator is conservative**: on average it is at
least as large as the true variance of difference-in-means.

Hint: quote `neymanVarEst_bias` and then `fpVar_nonneg`. -/
theorem var_DiM_le_expect_neymanVarEst {n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n)
    (hn₁ : 1 < n₁) (hn₀ : 1 < n - n₁) :
    (completeRandomization n n₁ hle).var (DiM P)
      ≤ (completeRandomization n n₁ hle).expect (neymanVarEst P) := by
  sorry

Show solution

New tabOpens with 1,809 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (core).  **The Neyman variance estimator is conservative**: on average it is at
least as large as the true variance of difference-in-means.

Hint: quote `neymanVarEst_bias` and then `fpVar_nonneg`. -/
theorem var_DiM_le_expect_neymanVarEst {n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n)
    (hn₁ : 1 < n₁) (hn₀ : 1 < n - n₁) :
    (completeRandomization n n₁ hle).var (DiM P)
      ≤ (completeRandomization n n₁ hle).expect (neymanVarEst P) := by
  rw [neymanVarEst_bias hle P hn₁ hn₀]
  have hS : (0 : ℝ) ≤ P.Stausq := fpVar_nonneg _
  have hnn : (0 : ℝ) ≤ (n : ℝ) := Nat.cast_nonneg n
  have : (0 : ℝ) ≤ P.Stausq / (n : ℝ) := div_nonneg hS hnn
  linarith

It is worth being precise about what has just been proved, because “conservative” is used loosely in the literature and the Lean is not loose at all. The statement is

Var ({\hat{τ}}_{DiM}) \leq E [\hat{V}],

an inequality between two real numbers, both of them fixed once the population and the design are fixed. It is not a statement about a limit, not a statement about coverage, and — this is the one people slide past — not a statement about any particular realization of $\hat{V}$ . The four-unit example makes all three numbers concrete:

The gap on the worked exampleNew tabOpens with 2,009 lines of library preamble — the code above is at the bottom of the editor.

/-- The Neyman variance estimator on the same example has expectation
`S₁²/2 + S₀²/2 = 1`, which overshoots the truth `2/3` by `S_τ²/n = (4/3)/4 = 1/3`.

The gap is not an artefact of small `n`: it is the un-identified `S_τ²` term, and
no estimator built from the observed data can remove it. -/
theorem exP_expect_neymanVarEst :
    (completeRandomization 4 2 (by norm_num)).expect (neymanVarEst exP) = 1 := by
  obtain ⟨h1, h0, _⟩ := exP_variances
  rw [completeRandomization_expect_neymanVarEst (by norm_num) exP (by norm_num) (by norm_num),
    h1, h0]
  norm_num

/-- The over-estimation on the worked example, spelled out: `1 - 2/3 = 1/3`. -/
theorem exP_neyman_gap :
    (completeRandomization 4 2 (by norm_num)).expect (neymanVarEst exP)
      - (completeRandomization 4 2 (by norm_num)).var (DiM exP) = 1 / 3 := by
  rw [exP_expect_neymanVarEst, exP_var_DiM]
  norm_num

$E [\hat{V}] = 1$ , $Var ({\hat{τ}}_{DiM}) = 2 / 3$ , and the gap is $S_{τ}^{2} / n = (4 / 3) / 4 = 1 / 3$ exactly — exP_neyman_gap is that subtraction, machine-checked. Meanwhile the assignment that treats units 2 and 3 reports $\hat{V} = 0$ , which is below the true variance. Conservativeness lives entirely on the left of $E$ ; draw by draw the estimator is free to be anything.

Finally, the inequality is sharp, with an equality condition that is exactly the classical one:

Exactly unbiased iff constant effectsNew tabOpens with 1,823 lines of library preamble — the code above is at the bottom of the editor.

/-- **The Neyman estimator is exactly unbiased iff the treatment effect is constant.**

The over-estimation `S_τ²/n` vanishes precisely when every unit has the same individual
treatment effect — the "constant effects" or "sharp null of no *heterogeneous* effect" case.
Under any effect heterogeneity at all, the reported variance is strictly too large, so
confidence intervals built from it are conservative but valid. -/
theorem expect_neymanVarEst_eq_var_iff {n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n)
    (hn₁ : 1 < n₁) (hn₀ : 1 < n - n₁) :
    (completeRandomization n n₁ hle).expect (neymanVarEst P)
        = (completeRandomization n n₁ hle).var (DiM P)
      ↔ ∀ i j, P.effect i = P.effect j := by
  have hn : 1 < n := by omega
  have hnR : (0 : ℝ) < (n : ℝ) := by
    have : (1 : ℝ) < (n : ℝ) := by exact_mod_cast hn
    linarith
  rw [neymanVarEst_bias hle P hn₁ hn₀]
  constructor
  · intro h
    have hzero : P.Stausq = 0 := by
      have hq : P.Stausq / (n : ℝ) = 0 := by linarith
      exact (div_eq_zero_iff.mp hq).resolve_right hnR.ne'
    exact (fpVar_eq_zero_iff hn P.effect).mp hzero
  · intro h
    have hzero : P.Stausq = 0 := (fpVar_eq_zero_iff hn P.effect).mpr h
    rw [hzero, zero_div, add_zero]

∀ i j, P.effect i = P.effect j is “the treatment effect is the same for every unit” — the constant-effects assumption. It is strictly weaker than Fisher’s sharp null of no effect at all, which is the special case $τ_{i} = 0$ , and it is not comparable to the weak null $τ = 0$ : constant effects with $τ_{i} = 5$ says nothing about the average being zero, and an average of zero says nothing about homogeneity. Getting that relationship backwards is easy in prose and impossible in Lean, where the two hypotheses are visibly different terms. The ↔ is the interesting part: under any effect heterogeneity whatsoever, the reported variance is strictly too large, so intervals built from it over-cover. Conservative, but never anti-conservative.

ExerciseConstant effects (stretch)

If $y_{i} (1) = y_{i} (0) + c$ for all $i$ then $S_{1}^{2} = S_{0}^{2}$ and $S_{τ}^{2} = 0$ , so Neyman’s formula collapses to $S_{0}^{2} \cdot n / (n_{1} n_{0})$ — the two-sample expression a first course would write. In Lean: get P.y1 = fun i => P.y0 i + c by funext and linarith from the hypothesis, then fpVar_add_const for $S_{1}^{2}$ and fpVar_eq_zero_iff for $S_{τ}^{2}$ , then field_simp; ring. Note that 1 < n has to be derived from 0 < n₁ and n₁ < n by omega before fpVar_eq_zero_iff will apply.

New tabOpens with 2,043 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (stretch).  **Under a constant treatment effect, Neyman's formula
collapses to the textbook two-sample expression.**

If `y1 i = y0 i + c` for every unit then `S₁² = S₀²` and `S_τ² = 0`, so

`Var(DiM) = S₀²/n₁ + S₀²/n₀ = S₀² · n / (n₁ n₀)`,

and by `expect_neymanVarEst_eq_var_iff` the reported variance is then exactly
right rather than merely conservative.

Hint: rewrite `P.y1` as `fun i => P.y0 i + c` with `funext`, use
`fpVar_add_const` for `S₁²` and `fpVar_eq_zero_iff` for `S_τ²`, then
`field_simp`. -/
theorem var_DiM_of_constant_effect {n n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n)
    (hpos : 0 < n₁) (hlt : n₁ < n) (c : ℝ) (hc : ∀ i, P.effect i = c) :
    (completeRandomization n n₁ hle).var (DiM P)
      = P.S0sq * ((n : ℝ) / ((n₁ : ℝ) * ((n : ℝ) - (n₁ : ℝ)))) := by
  sorry

Show solution

New tabOpens with 2,043 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (stretch).  **Under a constant treatment effect, Neyman's formula
collapses to the textbook two-sample expression.**

If `y1 i = y0 i + c` for every unit then `S₁² = S₀²` and `S_τ² = 0`, so

`Var(DiM) = S₀²/n₁ + S₀²/n₀ = S₀² · n / (n₁ n₀)`,

and by `expect_neymanVarEst_eq_var_iff` the reported variance is then exactly
right rather than merely conservative.

Hint: rewrite `P.y1` as `fun i => P.y0 i + c` with `funext`, use
`fpVar_add_const` for `S₁²` and `fpVar_eq_zero_iff` for `S_τ²`, then
`field_simp`. -/
theorem var_DiM_of_constant_effect {n n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n)
    (hpos : 0 < n₁) (hlt : n₁ < n) (c : ℝ) (hc : ∀ i, P.effect i = c) :
    (completeRandomization n n₁ hle).var (DiM P)
      = P.S0sq * ((n : ℝ) / ((n₁ : ℝ) * ((n : ℝ) - (n₁ : ℝ)))) := by
  have hn1 : 1 < n := by omega
  have hnR : (n₁ : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr (by omega)
  have hn0R : (n : ℝ) - (n₁ : ℝ) ≠ 0 := by
    have : (n₁ : ℝ) < (n : ℝ) := by exact_mod_cast hlt
    linarith
  have hy1 : P.y1 = fun i => P.y0 i + c := by
    funext i
    have := hc i
    unfold Population.effect at this
    linarith
  have hS1 : P.S1sq = P.S0sq := by
    rw [Population.S1sq, Population.S0sq, hy1, fpVar_add_const (by omega)]
  have hSt : P.Stausq = 0 :=
    (fpVar_eq_zero_iff hn1 P.effect).mpr fun i j => by rw [hc i, hc j]
  rw [neyman_variance hle P hpos hlt, hS1, hSt]
  field_simp
  ring

The general Horvitz–Thompson variance

Nothing above used a property of complete randomization until step 3, so the same route gives the Horvitz–Thompson variance for an arbitrary design. The difference is step 2. Horvitz–Thompson is linear in the indicators on every assignment, not just on a support, because it involves no group sizes: HT_eq_linear has no hypothesis on z at all. The proof of var_HT still goes through var_congr_of_support, but look at how it is called — fun z _ => HT_eq_linear D P z throws the support hypothesis away with an underscore. That is the shape to look for. A _ where a proof of D.prob z ≠ 0 was expected tells you the identity is unconditional and the support lemma is doing nothing but the rewriting.

Var(HT) for any designNew tabOpens with 1,855 lines of library preamble — the code above is at the bottom of the editor.

/-- The coefficient vector that makes Horvitz–Thompson a linear statistic in the indicators:
`aᵢ = (1/n)(y1ᵢ/πᵢ + y0ᵢ/(1 - πᵢ))`. -/
noncomputable def htCoef (D : Design n) (P : Population n) (i : Fin n) : ℝ :=
  (1 / (n : ℝ)) * (P.y1 i / D.propensity i + P.y0 i / (1 - D.propensity i))

/-- Horvitz–Thompson is a linear statistic in the indicators plus a constant — on *every*
assignment, not just on a support, because no group sizes are involved. -/
theorem HT_eq_linear (D : Design n) (P : Population n) (z : Assignment n) :
    HT D P z
      = (∑ i, htCoef D P i * Z z i)
        + -((1 / (n : ℝ)) * ∑ i, P.y0 i / (1 - D.propensity i)) := by
  have h1 : HT D P z = ∑ i, (1 / (n : ℝ)) * (Z z i * P.y1 i / D.propensity i
      - (1 - Z z i) * P.y0 i / (1 - D.propensity i)) := Finset.mul_sum _ _ _
  calc HT D P z
      = ∑ i, (1 / (n : ℝ)) * (Z z i * P.y1 i / D.propensity i
          - (1 - Z z i) * P.y0 i / (1 - D.propensity i)) := h1
    _ = ∑ i, (htCoef D P i * Z z i
          + -((1 / (n : ℝ)) * (P.y0 i / (1 - D.propensity i)))) :=
        Finset.sum_congr rfl fun i _ => by unfold htCoef; ring
    _ = (∑ i, htCoef D P i * Z z i)
          + ∑ i, -((1 / (n : ℝ)) * (P.y0 i / (1 - D.propensity i))) := Finset.sum_add_distrib
    _ = (∑ i, htCoef D P i * Z z i)
          + -((1 / (n : ℝ)) * ∑ i, P.y0 i / (1 - D.propensity i)) := by
        rw [Finset.sum_neg_distrib, ← Finset.mul_sum]

/-- **The Horvitz–Thompson variance for an arbitrary design**, in terms of the first- and
second-order inclusion probabilities:

`Var(HT) = ∑ᵢ ∑ⱼ aᵢ aⱼ (πᵢⱼ - πᵢ πⱼ)`,  `aᵢ = (1/n)(y1ᵢ/πᵢ + y0ᵢ/(1-πᵢ))`.

Specialising `πᵢⱼ` gives Bernoulli (`πᵢⱼ = p²`, so only the diagonal survives) and complete
randomisation (`πᵢⱼ = n₁(n₁-1)/(n(n-1))`, which is `neyman_variance` again). -/
theorem var_HT (D : Design n) (P : Population n) :
    D.var (HT D P)
      = ∑ i, ∑ j, htCoef D P i * htCoef D P j
          * (D.jointPropensity i j - D.propensity i * D.propensity j) := by
  have h := D.var_congr_of_support (f := HT D P)
    (g := fun z => (∑ i, htCoef D P i * Z z i)
      + -((1 / (n : ℝ)) * ∑ i, P.y0 i / (1 - D.propensity i)))
    (fun z _ => HT_eq_linear D P z)
  rw [h, D.var_add_const, D.var_linear]

Specializing $π_{i j}$ recovers everything: $π_{i j} = p^{2}$ kills the off-diagonal and gives Bernoulli, $π_{i j} = n_{1} (n_{1} - 1) / (n (n - 1))$ gives Neyman’s formula back.

ExerciseHorvitz–Thompson under Bernoulli (stretch)

Start from var_HT, then apply sum_sum_diag_offdiag with c = p - p² (from Design.jointPropensity_self plus bernoulli_propensity) and d = 0 (from bernoulli_jointPropensity). The d = 0 is independence, and it is why the answer is a single sum.

New tabOpens with 1,897 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (stretch).  The Horvitz–Thompson variance under the Bernoulli design collapses to
a single sum, because distinct units are independent.

Hint: start from `var_HT`, then use `sum_sum_diag_offdiag` with `c = p - p²` (the diagonal,
by `Design.jointPropensity_self`) and `d = 0` (the off-diagonal, by
`bernoulli_jointPropensity`). -/
theorem bernoulli_var_HT {p : ℝ} (hp0 : 0 ≤ p) (hp1 : p ≤ 1) (P : Population n) :
    (bernoulliDesign n p hp0 hp1).var (HT (bernoulliDesign n p hp0 hp1) P)
      = p * (1 - p) * ∑ i, (htCoef (bernoulliDesign n p hp0 hp1) P i) ^ 2 := by
  sorry

Show solution

New tabOpens with 1,897 lines of library preamble — the code above is at the bottom of the editor.

/-- Exercise (stretch).  The Horvitz–Thompson variance under the Bernoulli design collapses to
a single sum, because distinct units are independent.

Hint: start from `var_HT`, then use `sum_sum_diag_offdiag` with `c = p - p²` (the diagonal,
by `Design.jointPropensity_self`) and `d = 0` (the off-diagonal, by
`bernoulli_jointPropensity`). -/
theorem bernoulli_var_HT {p : ℝ} (hp0 : 0 ≤ p) (hp1 : p ≤ 1) (P : Population n) :
    (bernoulliDesign n p hp0 hp1).var (HT (bernoulliDesign n p hp0 hp1) P)
      = p * (1 - p) * ∑ i, (htCoef (bernoulliDesign n p hp0 hp1) P i) ^ 2 := by
  rw [var_HT, sum_sum_diag_offdiag (htCoef (bernoulliDesign n p hp0 hp1) P) (p - p ^ 2) 0
    (fun i j => (bernoulliDesign n p hp0 hp1).jointPropensity i j
      - (bernoulliDesign n p hp0 hp1).propensity i * (bernoulliDesign n p hp0 hp1).propensity j)
    (fun i => by rw [Design.jointPropensity_self, bernoulli_propensity]; ring)
    (fun i j hij => by
      rw [bernoulli_jointPropensity hij, bernoulli_propensity, bernoulli_propensity]; ring)]
  ring

Reading autoformalized Lean

Two formalizations of this chapter’s results that a language model will happily produce, and that a careful reader should reject.

Conservativeness, quantified over the wrong thing

-- WRONG: a plausible-looking "the estimator is conservative"
theorem neyman_conservative_bad {n n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n)
    (hn₁ : 1 < n₁) (hn₀ : 1 < n - n₁) (z : Assignment n) :
    (completeRandomization n n₁ hle).var (DiM P) ≤ neymanVarEst P z := by
  sorry

This reads well in English — “the Neyman estimator is at least as large as the true variance” — and it is false. The z in the binder makes it a claim about every realization; the true theorem is a claim about the expectation, and Design.expect is nowhere in the statement. The example population refutes it outright: exP_neymanVarEst_of_treated_two_three says $\hat{V} = 0$ on the assignment ${2, 3}$ , while exP_var_DiM says the variance is $2 / 3$ .

An “estimator” that is not a statistic

-- WRONG: unbiased, conservative, and useless
noncomputable def varEstBad {n : ℕ} (P : Population n) (n₁ : ℕ)
    (_z : Assignment n) : ℝ :=
  P.S1sq / (n₁ : ℝ) + P.S0sq / ((n : ℝ) - (n₁ : ℝ))

-- and then, trivially:
theorem varEstBad_conservative … :
    (completeRandomization n n₁ hle).var (DiM P)
      ≤ (completeRandomization n n₁ hle).expect (varEstBad P n₁) := by
  sorry

Every property you would check holds: it has the right type, its expectation is $S_{1}^{2} / n_{1} + S_{0}^{2} / n_{0}$ (it is constant, so its expectation is itself), and the conservativeness inequality is true of it. It is nonetheless not an estimator, because it is a function of the population, not of the data: computing it requires $y_{i} (1)$ for control units.

Lean will not stop you. Population is in scope inside any definition, and the type Assignment n → ℝ does not encode “depends only on what $z$ reveals”. Two tells: the assignment argument is unused — the underscore in _z is the conventional way to say so, and, less helpfully for a reader, it also silences the unused-variable linter that would otherwise have flagged it — and the body mentions P.y1 with no Z z i or ∑ i ∈ z restricting which units it is read for. The library’s honest point estimators carry a lemma making this checkable (HT_eq_observed and DiM_eq_observed, each showing the estimator factors through Yobs); neymanVarEst has no such lemma, so there the check has to be done by reading the definition — y1 summed over z, y0 over zᶜ.

Two habits worth keeping from this chapter. First, count the hypotheses and ask what each one is for: n₁ ≤ n was data for the design, 0 < n₁ < n was needed for Neyman, 1 < n₁ and 1 < n - n₁ were needed only once sample variances appeared. Second, when a statement about an estimator holds unconditionally, look for a division by zero that is quietly making a degenerate case true.