Variance — MrCLean

Under complete randomization, the variance of difference in means is an exact finite-population quantity rather than an asymptotic approximation. This chapter derives Neyman’s formula by first proving a general expression for the variance of a linear statistic $\sum_{i} a_{i} Z_{i}$ in terms of the inclusion probabilities $π_{i}$ and $π_{i j}$ . Bernoulli and complete-randomization formulas then follow as special cases.

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 812 lines of library preamble — the code above is at the bottom of the editor.

noncomputable def fpVar (v : Fin n → ℝ) : ℝ := (∑ i, (v i - mean v) ^ 2) / ((n : ℝ) - 1)

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 946 lines of library preamble — the code above is at the bottom of the editor.

noncomputable def S1sq : ℝ := fpVar P.y1

noncomputable def S0sq : ℝ := fpVar P.y0

noncomputable def Stausq : ℝ := fpVar P.effect

noncomputable def S10 : ℝ := fpCov P.y1 P.y0

$S_{10}$ , the finite-population covariance of $y (1)$ and $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.

Proposition. $S_{τ}^{2} = S_{1}^{2} + S_{0}^{2} - 2 S_{10}$ .

Proof.

The individual effects are $τ_{i} = y_{i} (1) - y_{i} (0)$ , so this is the usual “variance of a difference” identity applied to the two outcome vectors: for any $u, v$ ,

\frac{1}{n - 1} \sum_{i} ((u_{i} - v_{i}) - (\bar{u} - \bar{v}))^{2} = \frac{1}{n - 1} \sum_{i} ((u_{i} - \bar{u}) - (v_{i} - \bar{v}))^{2},

and expanding the square termwise gives

\frac{1}{n - 1} \sum_{i} (u_{i} - \bar{u})^{2} + \frac{1}{n - 1} \sum_{i} (v_{i} - \bar{v})^{2} - \frac{2}{n - 1} \sum_{i} (u_{i} - \bar{u}) (v_{i} - \bar{v}),

which is the claim with $u = y (1)$ and $v = y (0)$ .∎

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

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

```
have h : P.effect = fun i => P.y1 i - P.y0 i := rfl
```
effect is defined as the pointwise difference, so the two functions are equal by rfl. The have exists only so that the next rw has an equation to work with: rw needs a proof term, and rfl proofs of definitional identities are the cheapest way to produce one.
```
rw [Stausq, h]
```
The first two rewrites of one rw list. Stausq unfolds the left-hand side to fpVar P.effect, and h replaces P.effect by the explicit difference, so the goal is now about the variance of a difference of two outcome vectors.
Goal
```
⊢ (fpVar fun i => P.y1 i - P.y0 i) = P.S1sq + P.S0sq - 2 * P.S10
```
```
fpVar_sub
```
Third rewrite. fpVar_sub — the “variance of a difference” identity, itself proved as an exercise below by writing $u - v$ as $1 \cdot u + (- 1) \cdot v$ and quoting the linear-combination lemma — turns that into fpVar P.y1 + fpVar P.y0 - 2 * fpCov P.y1 P.y0.
Goal
```
⊢ fpVar P.y1 + fpVar P.y0 - 2 * fpCov P.y1 P.y0 = P.S1sq + P.S0sq - 2 * P.S10
```
```
rfl
```
The remaining goal compares that with P.S1sq + P.S0sq - 2 * P.S10, and those three names are definitions of exactly the three expressions on the left. The proof therefore ends where a textbook’s does: with the observation that both sides are the same thing under different names.

$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)

Scaling every entry of a vector by c multiplies its finite-population variance by c ^ 2: the variance is quadratic, not linear, in the data.

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

theorem fpVar_const_mul (c : ℝ) (v : Fin n → ℝ) :
    fpVar (fun i => c * v i) = c ^ 2 * fpVar v := by
  sorry

Hint

simp only [fpVar, mean_const_mul] pulls the c out of the mean; after ← mul_div_assoc and congr 1 the goal is a sum equal to c ^ 2 times a sum, so Finset.mul_sum followed by Finset.sum_congr rfl fun i _ => by ring proves it termwise — every summand is c ^ 2 * (v i - mean v) ^ 2.

Show solution

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

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)

The variance of a difference: fpVar (u - v) = fpVar u + fpVar v - 2 fpCov u v. This is the lemma that produces Stausq_eq above.

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

theorem fpVar_sub (u v : Fin n → ℝ) :
    fpVar (fun i => u i - v i) = fpVar u + fpVar v - 2 * fpCov u v := by
  sorry

Hint

Do not expand anything. Write u i - v i as 1 * u i + (-1) * v i as a have closed by funext i; ring, rewrite with it and with fpVar_linear_comb, and ring finishes.

Show solution

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

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)

Shifting every entry of a vector by the same constant leaves its finite-population variance unchanged. The hypothesis 0 < n is not decoration; part of the exercise is seeing where it gets spent.

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

theorem fpVar_add_const {n : ℕ} (hn : 0 < n) (v : Fin n → ℝ) (c : ℝ) :
    fpVar (fun i => v i + c) = fpVar v := by
  sorry

Hint

mean_add and mean_const give mean (fun i => v i + c) = mean v + c as a have — mean_const is where 0 < n is needed. After simp only [fpVar, hm] and congr 1, each summand is (v i + c - (mean v + c)) ^ 2, which Finset.sum_congr rfl fun i _ => by ring matches to (v i - mean v) ^ 2.

Show solution

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

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 299 lines of library preamble — the code above is at the bottom of the editor.

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.

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 chapter 7). 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,070 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

Proof.

The derivation is two lines of second-moment algebra. The variance of a sum is the double sum of covariances,

Var (\sum_{i} a_{i} Z_{i}) = \sum_{i} \sum_{j} Cov (a_{i} Z_{i}, a_{j} Z_{j}),

which is linearity of expectation applied to the expanded square. Covariance is bilinear, so the constants come out: $Cov (a_{i} Z_{i}, a_{j} Z_{j}) = a_{i} a_{j} Cov (Z_{i}, Z_{j})$ . Finally, for indicators,

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

by the definitions of the first- and second-order inclusion probabilities.∎

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

```
rw [var_sum]
```
Replace $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 an arbitrary index type, from cov_sum_sum — which in turn is the expansion of a product of two sums, Finset.sum_mul_sum, under an expectation. No design and no indicator has been mentioned yet.
Goal
```
⊢ (∑ i, ∑ j, D.cov (fun z => a i * Z z i) fun z => a j * Z z j) =
    ∑ i, ∑ j, a i * a j * (D.jointPropensity i j - D.propensity i * D.propensity j)
```
```
refine Finset.sum_congr rfl fun i _ => Finset.sum_congr rfl fun j _ => ?_
```
Finset.sum_congr 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. Applying it twice, once per layer, leaves a goal about a single pair $(i, j)$ .
Goal
```
i j : Fin n
⊢ (D.cov (fun z => a i * Z z i) fun z => a j * Z z j) =
    a i * a j * (D.jointPropensity i j - D.propensity i * D.propensity j)
```
```
rw [cov_const_mul_left, cov_const_mul_right]
```
Bilinearity, one argument at a time: cov_const_mul_left : Cov(c·f, g) = c·Cov(f, g) pulls out a i, and its right-hand sibling pulls out a j. What is left inside the covariance is the pair of bare indicators.
Goal
```
⊢ a i * (a j * D.cov (fun z => Z z i) fun z => Z z j) =
    a i * a j * (D.jointPropensity i j - D.propensity i * D.propensity j)
```
```
rw [cov_Z]
```
Name what is left. cov_Z : Cov(Z i, Z j) = πᵢⱼ - πᵢπⱼ is proved by cov_eq_expect_mul_sub_mul followed by rfl, because jointPropensity and propensity are defined as the two expectations that appear.
```
ring
```
The goal is now a i * (a j * X) = a i * a j * X with X an opaque real. Reassociating is exactly what ring is for, and using it here rather than mul_assoc by hand keeps the proof stable if the statement is later regrouped.

The textbook display splits the diagonal off.

Proposition. For any design and any fixed coefficients,

Var (\sum_{i} a_{i} Z_{i}) = \sum_{i} a_{i}^{2} π_{i} (1 - π_{i}) + \sum_{i} \sum_{j \neq i} a_{i} a_{j} (π_{i j} - π_{i} π_{j}) .

Proof.

Split each inner sum of the previous theorem at $j = i$ . The diagonal term is $a_{i}^{2} (π_{i i} - π_{i}^{2})$ , and $π_{i i} = E [Z_{i}^{2}] = E [Z_{i}] = π_{i}$ because an indicator is idempotent, so it equals $a_{i}^{2} π_{i} (1 - π_{i})$ — the variance of a single Bernoulli indicator. The rest of the inner sum is the off-diagonal term.∎

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

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

```
rw [var_linear, ← Finset.sum_add_distrib]
```
Start from the double sum just proved, then run Finset.sum_add_distrib : ∑ (f i + g i) = ∑ f i + ∑ g i backwards, merging the two sums on the right-hand side into one sum whose summand is diagonal i + offdiagonal i. Combining first means the rest of the proof compares two sums over the same index, one unit at a time.

refine Finset.sum_congr rfl fun i _ => ?_

One unit at a time, as before. The goal is now: the $i$ -th row of the double sum equals its diagonal entry plus the row with $j = i$ removed.

Goal

i : Fin n
⊢ ∑ j, a i * a j * (D.jointPropensity i j - D.propensity i * D.propensity j) =
    a i ^ 2 * (D.propensity i * (1 - D.propensity i)) +
      ∑ j ∈ univ.erase i, a i * a j * (D.jointPropensity i j - D.propensity i * D.propensity j)

```
rw [← Finset.add_sum_erase _ _ (Finset.mem_univ i)]
```
Finset.add_sum_erase : a ∈ s → f a + ∑ x ∈ s.erase a, f x = ∑ x ∈ s, f x, used right to left, splits the left-hand sum at $j = i$ . univ.erase i is 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. The membership side condition is Finset.mem_univ i, which is free.

congr 1

Both sides are now something + (the same erased sum), so congr 1 discards the matching summand and leaves only the diagonal entries to be compared.

Goal

case e_a
⊢ a i * a i * (D.jointPropensity i i - D.propensity i * D.propensity i) =
    a i ^ 2 * (D.propensity i * (1 - D.propensity i))

```
rw [jointPropensity_self]; ring
```
jointPropensity_self : πᵢᵢ = πᵢ is the idempotence step of the informal proof (it is the warm-up exercise below, and its content is $Z_{i}^{2} = Z_{i}$ ). After it the goal is a i * a i * (π i - π i * π i) = a i ^ 2 * (π i * (1 - π i)), which is a polynomial identity, so 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.

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

The joint propensity of a unit with itself is its propensity, because $Z_{i}^{2} = Z_{i}$ .

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

theorem jointPropensity_self (D : Design n) (i : Fin n) :
    D.jointPropensity i i = D.propensity i := by
  sorry

Hint

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

Show solution

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

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)

Adding a constant to a statistic does not change its variance. This lemma is what lets Neyman’s proof throw away the - (∑ y0)/n₀ tail of difference in means.

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

theorem var_add_const (D : Design n) (f : Assignment n → ℝ) (c : ℝ) :
    D.var (fun z => f z + c) = D.var f := by
  sorry

Hint

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.

Show solution

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

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)

Scaling a statistic scales its variance by the square: $Var (c \cdot f) = c^{2} Var (f)$ .

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

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

Hint

Four rewrites and a ring. Go through var_eq_cov_self so that cov_const_mul_left and cov_const_mul_right apply, then come back with ← var_eq_cov_self.

Show solution

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

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.

Proposition. If $π_{i} = p$ for every $i$ and $π_{i j} = q$ for every $i \neq j$ , then

Var (\sum_{i} a_{i} Z_{i}) = (p - p^{2}) \sum_{i} a_{i}^{2} + (q - p^{2}) ((\sum_{i} a_{i})^{2} - \sum_{i} a_{i}^{2}) .

Proof.

In the double sum $\sum_{i} \sum_{j} a_{i} a_{j} (π_{i j} - π_{i} π_{j})$ the kernel $π_{i j} - π_{i} π_{j}$ equals $p - p^{2}$ on the diagonal and $q - p^{2}$ off it, so

\sum_{i} \sum_{j} a_{i} a_{j} (π_{i j} - π_{i} π_{j}) = (p - p^{2}) \sum_{i} a_{i}^{2} + (q - p^{2}) \sum_{i} \sum_{j \neq i} a_{i} a_{j} .

For the second double sum, the total $\sum_{j} a_{j}$ does not depend on $i$ , so $\sum_{j \neq i} a_{i} a_{j} = a_{i} (\sum_{j} a_{j} - a_{i})$ , and summing over $i$ gives $(\sum_{i} a_{i})^{2} - \sum_{i} a_{i}^{2}$ . That purely combinatorial identity is isolated in the library as sum_sum_diag_offdiag, and it is the only thing the two designs below need.∎

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

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)

```
rw [var_linear]
```
Get to the double sum. Everything specific to an exchangeable design happens after this line; before it, the general theorem is doing all the work.
```
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) …
```
Apply the combinatorial identity, instantiating its kernel F with the design’s own $π_{i j} - π_{i} π_{j}$ , its diagonal value c with $p - p^{2}$ and its off-diagonal value d with $q - p^{2}$ . The lemma then is the goal; there is nothing left but its two side conditions.
```
fun i => by rw [jointPropensity_self, hp i]; ring
```
The diagonal side condition, ∀ i, F i i = c. Rewriting $π_{i i}$ to $π_{i}$ and then $π_{i}$ to $p$ leaves p - p * p = p - p ^ 2, which is ring.
```
fun i j hij => by rw [hq i j hij, hp i, hp j]; ring
```
The off-diagonal side condition, ∀ i j, i ≠ j → F i j = d. The hypothesis hij : i ≠ j is exactly what hq needs; after the three rewrites the goal is q - p * p = q - p ^ 2. Note that the statement of the theorem quantifies hq over distinct pairs only — asking for $π_{i j} = q$ at $i = j$ as well would be asking for $p = q$ , which holds for no interesting design.

Corollary (Bernoulli). Under independent assignment with common probability $p$ , $q = p^{2}$ , the off-diagonal term vanishes and $Var (\sum_{i} a_{i} Z_{i}) = p (1 - p) \sum_{i} a_{i}^{2}$ .

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

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

```
rw [Design.var_linear_of_exchangeable _ a p (p ^ 2) (fun i => bernoulli_propensity i) (fun i j hij => bernoulli_jointPropensity hij)]
```
One rewrite with the exchangeable formula, telling it that this design’s two constants are $p$ and $p^{2}$ . The two fun arguments are the proofs — bernoulli_propensity : π i = p and bernoulli_jointPropensity : i ≠ j → πᵢⱼ = p ^ 2, both from chapter 7. Everything design-specific is in those two lemmas.
```
ring
```
The goal is (p - p ^ 2) * ∑ a i ^ 2 + (p ^ 2 - p ^ 2) * (…) = p * (1 - p) * ∑ a i ^ 2. The second coefficient is literally $p^{2} - p^{2}$ , so ring deletes that whole term and matches the first. That vanishing is what “independent assignment” means here; it is not assumed, it is computed.

Corollary (complete randomization). With $n_{1}$ treated out of $n$ , $0 < n_{1} < n$ ,

Var (\sum_{i} a_{i} Z_{i}) = \frac{n_{1} n_{0}}{n} S_{a}^{2}, n_{0} = n - n_{1},

where $S_{a}^{2} = \frac{1}{n - 1} \sum_{i} (a_{i} - \bar{a})^{2}$ is the finite-population variance of the coefficient vector.

Proof.

Here $p = n_{1} / n$ and $q = \frac{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. The negative term is exactly the right size to convert $\sum_{i} a_{i}^{2}$ into the centered sum of squares. Concretely, with $T = \sum_{i} a_{i}$ and $Q = \sum_{i} a_{i}^{2}$ , the exchangeable formula reads $(p - p^{2}) Q + (q - p^{2}) (T^{2} - Q)$ , and substituting and clearing denominators turns it into $\frac{n_{1} n_{0}}{n} \cdot \frac{n Q - T^{2}}{n (n - 1)}$ , which is $\frac{n_{1} n_{0}}{n} S_{a}^{2}$ because $\sum_{i} (a_{i} - \bar{a})^{2} = Q - T^{2} / n$ .∎

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

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

```
have hn2 : 2 ≤ n := by omega
```
From 0 < n₁ and n₁ < n there are at least two units, which is what makes $n - 1 \neq 0$ below. omega decides goals of linear arithmetic over ℕ and ℤ from the hypotheses in context, and this is the kind of bookkeeping it exists for.
```
have hn : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr (by omega); have hn1 : (n : ℝ) - 1 ≠ 0 := …
```
The two non-vanishing denominators, in $R$ . The second has to be done by casting hn2 and calling linarith, because (n : ℝ) - 1 is real subtraction: the natural-number statement n - 1 ≠ 0 is a different claim, and at $n = 0$ it would be false while the real one would not.

rw [Design.var_linear_of_exchangeable _ a ((n₁ : ℝ) / (n : ℝ)) ((n₁ : ℝ) * ((n₁ : ℝ) - 1) / ((n : ℝ) * ((n : ℝ) - 1))) …]

The exchangeable formula again, this time with $p = n_{1} / n$ and $q = n_{1} (n_{1} - 1) / (n (n - 1))$ — the two propensity lemmas proved in chapter 7 by counting subsets. Everything from here on is algebra with no probability in it.

Goal

⊢ (↑n₁ / ↑n - (↑n₁ / ↑n) ^ 2) * ∑ i, a i ^ 2 +
      (↑n₁ * (↑n₁ - 1) / (↑n * (↑n - 1)) - (↑n₁ / ↑n) ^ 2) * ((∑ i, a i) ^ 2 - ∑ i, a i ^ 2) =
    ↑n₁ * (↑n - ↑n₁) / ↑n * fpVar a

fpVar_eq_sum_sq hn

Second entry of the same rw list: rewrite the target’s fpVar a into the computational form $\frac{n \sum a_{i}^{2} - (\sum a_{i})^{2}}{n (n - 1)}$ . That form has no inner mean, hence no nested / n, which is what makes the next line possible.

Goal

⊢ (↑n₁ / ↑n - (↑n₁ / ↑n) ^ 2) * ∑ i, a i ^ 2 +
      (↑n₁ * (↑n₁ - 1) / (↑n * (↑n - 1)) - (↑n₁ / ↑n) ^ 2) * ((∑ i, a i) ^ 2 - ∑ i, a i ^ 2) =
    ↑n₁ * (↑n - ↑n₁) / ↑n * ((↑n * ∑ i, a i ^ 2 - (∑ i, a i) ^ 2) / (↑n * (↑n - 1)))

```
field_simp
```
Clear all four denominators, using hn and hn1 from the context. What is left is a polynomial identity in $n$ , $n_{1}$ , $\sum a_{i}$ and $\sum a_{i}^{2}$ .
```
ring
```
Check that identity. This is the cancellation the informal proof waved at, and it is the whole reason $\frac{n_{1} n_{0}}{n}$ comes out in front of a centered sum of squares rather than a raw one.

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,203 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
  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 _ _
  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

Proof.

Write $n_{0} = n - n_{1}$ and let $z$ be any assignment with $| z | = n_{1}$ . The treated average is $\frac{1}{n_{1}} \sum_{i} Z_{i} y_{i} (1)$ and the control average is $\frac{1}{n_{0}} \sum_{i} (1 - Z_{i}) y_{i} (0) = \frac{1}{n_{0}} \sum_{i} y_{i} (0) - \frac{1}{n_{0}} \sum_{i} Z_{i} y_{i} (0)$ , so

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

Complete randomization gives probability zero to every other assignment, and adding a constant does not change a variance, so $Var ({\hat{τ}}_{DiM}) = Var (\sum_{i} a_{i} Z_{i}) = \frac{n_{1} n_{0}}{n} S_{a}^{2}$ by the corollary above.

It remains to expand $S_{a}^{2}$ . The finite-population variance is a quadratic form, so for $a = \frac{1}{n_{1}} y (1) + \frac{1}{n_{0}} y (0)$ ,

S_{a}^{2} = \frac{S_{1}^{2}}{n_{1}^{2}} + \frac{S_{0}^{2}}{n_{0}^{2}} + \frac{2 S_{10}}{n_{1} n_{0}},

and multiplying by $\frac{n_{1} n_{0}}{n}$ gives $\frac{n_{0} S_{1}^{2}}{n n_{1}} + \frac{n_{1} S_{0}^{2}}{n n_{0}} + \frac{2 S_{10}}{n}$ . Substituting $2 S_{10} = S_{1}^{2} + S_{0}^{2} - S_{τ}^{2}$ and collecting the $S_{1}^{2}$ terms, their coefficient is $\frac{n_{0}}{n n_{1}} + \frac{1}{n} = \frac{n_{0} + n_{1}}{n n_{1}} = \frac{1}{n_{1}}$ ; likewise for $S_{0}^{2}$ ; and $- S_{τ}^{2} / n$ is left over.∎

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,167 lines of library preamble — the code above is at the bottom of the editor.

noncomputable def dimCoef (P : Population n) (n₁ n₀ : ℝ) (i : Fin n) : ℝ :=
  (1 / n₁) * P.y1 i + (1 / n₀) * P.y0 i

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.

```
have hcompl : ((#zᶜ : ℕ) : ℝ) = (n : ℝ) - (n₁ : ℝ) := by rw [card_compl, hz, Nat.cast_sub hlt.le]
```
The size of the control group as a real number, and the one place natural subtraction has to be handled. card_compl : #zᶜ = n - #z lives in ℕ, where subtraction truncates; hz replaces #z by n₁; and Nat.cast_sub, whose hypothesis is n₁ ≤ n, is what licenses moving that subtraction across the cast.
```
unfold DiM; rw [sum_mem_eq_sum_Z_mul, sum_compl_eq_sum_one_sub_Z_mul, hz, hcompl]
```
First calc step: rewrite both group sums as indicator-weighted sums over the whole population. sum_mem_eq_sum_Z_mul : ∑ i ∈ z, v i = ∑ i, Z z i * v i and its control-side counterpart sum_compl_eq_sum_one_sub_Z_mul are the two directions this chapter keeps using; hz and hcompl then replace the two group sizes by the constants n₁ and (n : ℝ) - (n₁ : ℝ). After this step the assignment appears only through Z z i.
Goal
```
⊢ (∑ i ∈ z, P.y1 i) / ↑(#z) - (∑ i ∈ zᶜ, P.y0 i) / ↑(#zᶜ) =
    (∑ i, Z z i * P.y1 i) / ↑n₁ - (∑ i, (1 - Z z i) * P.y0 i) / (↑n - ↑n₁)
```
```
rw [Finset.sum_div, Finset.sum_div, Finset.sum_sub_distrib]
```
Second calc step: push the two divisions inside their sums (Finset.sum_div : (∑ f i) / c = ∑ (f i / c)), then match the two sides with Finset.sum_sub_distrib : ∑ (f i - g i) = ∑ f i - ∑ g i — rw rewrites wherever the pattern occurs, and the sum of differences it finds is on the right. Everything is now a single sum over units, which is the shape a linear statistic has to have.
```
Finset.sum_congr rfl fun i _ => by unfold dimCoef; ring
```
Third step, and the only one with content: regroup the $i$ -th summand as dimCoef … i * Z z i + -(y0 i / n₀). Unfolding dimCoef makes the goal a polynomial identity in Z z i, y1 i, y0 i and the two group sizes, so ring closes it. This is where the plus sign in the coefficient comes from: the - (1 - Z z i) * y0 i / n₀ of the previous line is -y0 i / n₀ + Z z i * y0 i / n₀.
```
Finset.sum_add_distrib
```
Fourth step: split the sum of two things into two sums — the linear part and a sum of copies of the constant.
```
rw [Finset.sum_neg_distrib, Finset.sum_div]
```
Fifth step: pull the negation out of the second sum and recombine ∑ i, y0 i / n₀ into (∑ i, y0 i) / n₀, which is the additive constant as the statement writes it. The calc chain ends syntactically equal to the goal.

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} S_{a}^{2}$ .

Step 4 — expand. fpVar_dimCoef expands $S_{a}^{2}$ 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.

```
have hn : (n : ℝ) ≠ 0 := …; have hn1 : (n₁ : ℝ) ≠ 0 := …; have hn0 : (n : ℝ) - (n₁ : ℝ) ≠ 0 := …
```
The three denominators of the statement, shown to be non-zero before anything else happens. The first two are Nat.cast_ne_zero.mpr (by omega); the third has to be done in ℝ, by casting hlt and calling linarith.
```
rw [(completeRandomization n n₁ hle).var_congr_of_support (g := fun z => (∑ i, dimCoef P … i * Z z i) + -(…)) (fun z hz => DiM_eq_linear_of_card P hlt (card_eq_of_prob_ne_zero hz))]
```
Step 2, first half. var_congr_of_support needs the replacement statistic g spelled out — Lean cannot guess it — and a proof that DiM P equals it wherever the probability is non-zero. That proof is the composition of the support lemma with step 1: card_eq_of_prob_ne_zero turns D.prob z ≠ 0 into #z = n₁, which is precisely the hypothesis DiM_eq_linear_of_card wants.
```
exact (completeRandomization n n₁ hle).var_add_const _ _
```
Step 2, second half: var_add_const : Var(f + c) = Var(f) throws the constant tail away. That finishes the have hlin, which now equates the variance of DiM P with the variance of the linear statistic.

rw [hlin]

Substitute it into the main goal. What is left is the variance of the linear statistic, with no DiM in sight.

Goal

⊢ ((completeRandomization n n₁ hle).var fun z =>
      ∑ i, dimCoef P (↑n₁) (↑n - ↑n₁) i * Z z i) =
    P.S1sq / ↑n₁ + P.S0sq / (↑n - ↑n₁) - P.Stausq / ↑n

completeRandomization_var_linear hle hpos hlt

Step 3. The complete-randomization formula replaces that variance by $\frac{n_{1} n_{0}}{n} S_{a}^{2}$ , with $a$ the coefficient vector dimCoef P n₁ (n - n₁). This is the last time the design is mentioned.

Goal

⊢ ↑n₁ * (↑n - ↑n₁) / ↑n * fpVar (dimCoef P (↑n₁) (↑n - ↑n₁)) =
    P.S1sq / ↑n₁ + P.S0sq / (↑n - ↑n₁) - P.Stausq / ↑n

fpVar_dimCoef

Step 4, first half. fpVar_dimCoef is the quadratic-form expansion: it applies fpVar_linear_comb to $a = \frac{1}{n_{1}} y (1) + \frac{1}{n_{0}} y (0)$ and then recognizes the three pieces as S1sq, S0sq and S10 by rfl.

Goal

⊢ ↑n₁ * (↑n - ↑n₁) / ↑n *
      ((1 / ↑n₁) ^ 2 * P.S1sq + (1 / (↑n - ↑n₁)) ^ 2 * P.S0sq +
        2 * (1 / ↑n₁) * (1 / (↑n - ↑n₁)) * P.S10) =
    P.S1sq / ↑n₁ + P.S0sq / (↑n - ↑n₁) - P.Stausq / ↑n

```
P.Stausq_eq
```
Step 4, second half: eliminate the unestimable covariance by rewriting $S_{τ}^{2} = S_{1}^{2} + S_{0}^{2} - 2 S_{10}$ — the identity proved at the top of the chapter — in the goal’s right-hand side. Both sides are now built from the same three quantities.

field_simp

Clear the denominators $n$ , $n_{1}$ and $n - n_{1}$ , which is possible only because hn, hn1 and hn0 are in context. What is left is a polynomial identity in those three numbers and the three variances.

Goal

⊢ (↑n - ↑n₁) ^ 2 * P.S1sq + ↑n₁ ^ 2 * P.S0sq + ↑n₁ * (↑n - ↑n₁) * 2 * P.S10 =
    ↑n * ((↑n - ↑n₁) * P.S1sq + ↑n₁ * P.S0sq) -
      ↑n₁ * (↑n - ↑n₁) * (P.S1sq + P.S0sq - 2 * P.S10)

```
ring
```
Check it. Reading the two sides of the goal above against the informal proof: the left-hand side is $\frac{n_{1} n_{0}}{n} S_{a}^{2}$ with everything multiplied out, the right-hand side is the claimed formula with $S_{τ}^{2}$ substituted, and the identity is the collection of coefficients performed there by hand.

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,507 lines of library preamble — the code above is at the bottom of the editor.

noncomputable def exP : Population 4 where
  y1 := ![0, 2, 2, 2]
  y0 := ![0, 0, 0, 2]

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]

```
y1 := ![0, 2, 2, 2]
```
![…] is Mathlib’s notation for a vector literal, i.e. a function Fin 4 → ℝ given by listing its values. It is built from Matrix.vecCons, which is why the lemmas needed to evaluate it below are called Matrix.cons_val_*.
```
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]
```
One call, with everything it needs named. Population.tau, effect, exP and mean are unfolded to expose the arithmetic; Fin.sum_univ_four replaces ∑ i : Fin 4, f i by f 0 + f 1 + f 2 + f 3, since norm_num cannot evaluate a sum over a general Finset; and the four Matrix lemmas evaluate ![a, b, c, d] i at each index. The last two are needed because cons_val_two rewrites ![a, b, c, d] 2 to vecHead (vecTail ![b, c, d]) rather than straight to c, and tail_cons/head_cons are what peel those off; none of the four is a simp lemma by default, which is why they are listed. What is left, $(0 + 2 + 2 + 0) / 4 = 1$ , is arithmetic.

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$ . Make Lean agree.

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

theorem exP_variances : exP.S1sq = 1 ∧ exP.S0sq = 1 ∧ exP.Stausq = 4 / 3 := by
  sorry

Hint

Split the conjunction with refine ⟨?_, ?_, ?_⟩ <;> and hand all three goals to one norm_num. Everything in sight is a definition, so Population.S1sq, Population.S0sq, Population.Stausq, Population.effect, fpVar, mean and exP all go in its simp set, with Fin.sum_univ_four to expand the sums; ![a, b, c, d] 2 needs Matrix.cons_val_two (and its cons_val_three sibling) to reduce — those are not simp lemmas.

Show solution

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

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)

With $n = 4$ and $n_{1} = n_{0} = 2$ , Neyman’s formula on the example reads $\frac{1}{2} + \frac{1}{2} - \frac{4 / 3}{4} = \frac{2}{3}$ .

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

theorem exP_var_DiM :
    (completeRandomization 4 2 (by norm_num)).var (DiM exP) = 2 / 3 := by
  sorry

Hint

obtain ⟨h1, h0, hτ⟩ := exP_variances, then rewrite with neyman_variance — it wants four arguments, exP itself and three side conditions each discharged by by norm_num — and with h1, h0, hτ; norm_num finishes. 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.

Show solution

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

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,522 lines of library preamble — the code above is at the bottom of the editor.

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]

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]

```
have hcompl : ({2, 3} : Assignment 4)ᶜ = {0, 1} := by decide
```
The control group of this particular assignment, computed rather than argued. Assignment 4 is Finset (Fin 4), equality of such sets is decidable, and decide asks the kernel to run that decision procedure. This works only because everything is a closed term of a finite type.
```
rw [hcompl, Finset.sum_pair (by decide : (2 : Fin 4) ≠ 3), Finset.sum_pair (by decide : (0 : Fin 4) ≠ 1), Finset.card_pair …, Finset.card_pair …]
```
Evaluate the four ingredients of DiM. Finset.sum_pair : a ≠ b → ∑ i ∈ {a, b}, f i = f a + f b and Finset.card_pair : a ≠ b → #{a, b} = 2 both need the two elements to be distinct, and that side condition is again decidable. Writing (by decide : (2 : Fin 4) ≠ 3) with the type ascription is how you tell Lean which disequality to prove.
```
norm_num [exP, Matrix.cons_val_two, Matrix.cons_val_three, Matrix.head_cons, Matrix.tail_cons]
```
Arithmetic. The goal is now $(2 + 2) / 2 - (0 + 0) / 2 = 2$ once the vector literals are evaluated, which is what the four Matrix lemmas do. exP_neymanVarEst_of_treated_two_three is the same four lines applied to neymanVarEst instead: both treated units have $y (1) = 2$ and both control units $y (0) = 0$ , so both sample variances are $0$ and so is the reported variance.

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,243 lines of library preamble — the code above is at the bottom of the editor.

noncomputable def sampleVar (v : Fin n → ℝ) (s : Finset (Fin n)) : ℝ :=
  (∑ i ∈ s, (v i - (∑ j ∈ s, v j) / (#s : ℝ)) ^ 2) / ((#s : ℝ) - 1)

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 = 1 - Z$ for the control group.

Theorem (sample variance of a simple random sample without replacement). Let $W_{1}, \dots, W_{n}$ be ${0, 1}$ -valued random variables with $E [W_{i}] = m / n$ for every $i$ and $E [W_{i} W_{j}] = \frac{m (m - 1)}{n (n - 1)}$ for every $i \neq j$ , and let $v$ be a fixed vector. Then

E [\sum_{i} W_{i} v_{i}^{2} - \frac{1}{m} (\sum_{i} W_{i} v_{i})^{2}] = (m - 1) S_{v}^{2},

where $S_{v}^{2} = \frac{1}{n - 1} \sum_{i} (v_{i} - \bar{v})^{2}$ . Dividing by $m - 1$ : the sample variance of a simple random sample of size $m$ drawn without replacement is unbiased for the finite-population variance, which is the classical reason the $n - 1$ is where it is.

Proof.

Write $T = \sum_{i} v_{i}$ and $Q = \sum_{i} v_{i}^{2}$ . The first term is immediate from linearity: $E [\sum_{i} W_{i} v_{i}^{2}] = \frac{m}{n} Q$ .

For the second, expand the square as a double sum, $(\sum_{i} W_{i} v_{i})^{2} = \sum_{i} \sum_{j} v_{i} v_{j} W_{i} W_{j}$ , and take expectations termwise. The diagonal has $E [W_{i}^{2}] = E [W_{i}] = m / n$ , because $W_{i}$ is ${0, 1}$ -valued, and the off-diagonal has $\frac{m (m - 1)}{n (n - 1)}$ ; since $\sum_{i} \sum_{j \neq i} v_{i} v_{j} = T^{2} - Q$ ,

E [(\sum_{i} W_{i} v_{i})^{2}] = \frac{m}{n} Q + \frac{m (m - 1)}{n (n - 1)} (T^{2} - Q) .

Subtracting one- $m$ -th of that from the first term,

\frac{m}{n} Q - \frac{1}{m} [\frac{m}{n} Q + \frac{m (m - 1)}{n (n - 1)} (T^{2} - Q)] = \frac{m - 1}{n} Q - \frac{m - 1}{n (n - 1)} (T^{2} - Q) = \frac{m - 1}{n (n - 1)} (n Q - T^{2}),

and $n Q - T^{2} = n (n - 1) S_{v}^{2}$ , so the whole thing is $(m - 1) S_{v}^{2}$ .∎

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

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
  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]
  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

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 the disjunction $W_{i} = 0$ or $W_{i} = 1$ is what keeps the proof inside ring.

```
have hA : D.expect (fun z => ∑ i, W z i * (v i) ^ 2) = (m / (n : ℝ)) * ∑ i, (v i) ^ 2
```
The first of the two moment computations, stated as a have because both of them are needed before anything can be combined. Its statement is the informal proof’s $E [\sum_{i} W_{i} v_{i}^{2}] = \frac{m}{n} Q$ ; the goal shown here is the theorem’s, before either have.
Goal
```
⊢ (D.expect fun z => ∑ i, W z i * v i ^ 2 - (∑ i, W z i * v i) ^ 2 / m) = (m - 1) * fpVar v
```
```
rw [D.expect_sum univ (fun i z => W z i * (v i) ^ 2), Finset.mul_sum]
```
Exchange $E$ and $\sum_{i}$ on the left, and push the constant $m / n$ into the sum on the right, so that both sides are sums over units and can be compared termwise.
Goal
```
⊢ (∑ i, D.expect fun z => W z i * v i ^ 2) = ∑ i, m / ↑n * v i ^ 2
```
```
exact Finset.sum_congr rfl fun i _ => by rw [D.expect_mul_const (fun z => W z i) ((v i) ^ 2), h1 i]
```
One unit at a time: expect_mul_const pulls the constant v i ^ 2 out of the expectation, and h1 i replaces the remaining E[W i] by m / n. That is the whole of hA.
```
have hsq : ∀ z, (∑ i, W z i * v i) ^ 2 = ∑ i, ∑ j, v i * v j * (W z i * W z j)
```
The algebraic identity behind the second moment, proved pointwise in z so that it can be pushed through the expectation by expect_congr. rw [sq, Finset.sum_mul_sum] turns the square into the double sum ∑ i, ∑ j, (W z i * v i) * (W z j * v j), and the goal shown is what remains: reordering each factor, which ring does term by term inside two nested Finset.sum_congrs.
Goal
```
z : Assignment n
⊢ ∑ i, ∑ j, W z i * v i * (W z j * v j) = ∑ i, ∑ j, v i * v j * (W z i * W z j)
```
```
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 …]
```
Move the expectation inside both sums. expect_congr hsq replaces the squared statistic by the double sum under the binder; the first expect_sum exchanges $E$ with the outer $\sum_{i}$ ; and the sum_congr applies expect_sum again, once per row, to exchange it with the inner $\sum_{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)

Pull the constant v i * v j out of every one of the $n^{2}$ expectations. What is left inside is exactly E[W i * W j], which is what the two moment hypotheses talk about.

Goal

⊢ (∑ i, ∑ j, v i * v j * D.expect fun z => W z i * W z j) =
    m / ↑n * ∑ i, v i ^ 2 + m * (m - 1) / (↑n * (↑n - 1)) * ((∑ i, v i) ^ 2 - ∑ i, v i ^ 2)

```
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)) … …
```
The same combinatorial identity used for exchangeable designs, applied now with the kernel E[W i * W j]: its diagonal value is $m / n$ and its off-diagonal value $\frac{m (m - 1)}{n (n - 1)}$ , and the identity supplies $\sum_{i} \sum_{j \neq i} v_{i} v_{j} = T^{2} - Q$ . The diagonal side condition is by rw [D.expect_congr (fun z => hidem z i), h1 i] — idempotence under the binder, then h1 — and the off-diagonal one is h2 i j hij directly. This finishes hB.
```
rw [D.expect_sub _ _]
```
Back to the main goal, and the last four rewrites are the informal proof’s final display. expect_sub splits the expectation of the difference into the two pieces hA and hB were computed for.
Goal
```
⊢ ((D.expect fun z => ∑ i, W z i * v i ^ 2) -
      D.expect fun z => (∑ i, W z i * v i) ^ 2 / m) =
    (m - 1) * fpVar v
```
```
hA, D.expect_div_const _ m, hB
```
Substitute the first moment; move the / m outside the second expectation with expect_div_const : E[f / c] = E[f] / c; substitute the second moment. The goal is now free of expectations altogether.

fpVar_eq_sum_sq hn

Rewrite the target’s fpVar v as $\frac{n Q - T^{2}}{n (n - 1)}$ , the same computational form used in the complete-randomization proof, so that both sides are rational functions of $T$ , $Q$ , $n$ and $m$ .

Goal

⊢ m / ↑n * ∑ i, v i ^ 2 -
      (m / ↑n * ∑ i, v i ^ 2 +
        m * (m - 1) / (↑n * (↑n - 1)) * ((∑ i, v i) ^ 2 - ∑ i, v i ^ 2)) / m =
    (m - 1) * ((↑n * ∑ i, v i ^ 2 - (∑ i, v i) ^ 2) / (↑n * (↑n - 1)))

```
field_simp; ring
```
Clear the denominators — hn, hn1 and hm are the three non-vanishing facts in the hypothesis list, and this is what they are for — and check the resulting polynomial identity. This step is the two lines of cancellation the informal proof performs between its last two displays.

The engine now applies twice. The treated group of a completely randomized experiment is a simple random sample of size $n_{1}$ , and the control group is one of size $n_{0}$ ; in both cases the two moments the theorem asks for are exactly the propensity and joint propensity computed in chapter 7.

Theorem (the sample variance is unbiased). Under complete randomization with $1 < n_{1} < n$ , the treated sample variance is unbiased for $S_{1}^{2}$ : $E [s_{1}^{2}] = S_{1}^{2}$ . Symmetrically, if the control group has at least two units, $E [s_{0}^{2}] = S_{0}^{2}$ .

Proof.

For the treated group, apply the previous theorem with $W_{i} = Z_{i}$ , $m = n_{1}$ and $v = y (1)$ . Its three hypotheses hold: $Z_{i}^{2} = Z_{i}$ ; $E [Z_{i}] = π_{i} = n_{1} / n$ ; and $E [Z_{i} Z_{j}] = π_{i j} = \frac{n_{1} (n_{1} - 1)}{n (n - 1)}$ for $i \neq j$ . On any assignment with $| z | = n_{1}$ the treated sample variance is $\frac{1}{n_{1} - 1} (\sum_{i} Z_{i} y_{i} (1)^{2} - \frac{1}{n_{1}} (\sum_{i} Z_{i} y_{i} (1))^{2})$ , by the “sum of squares minus square of sum” identity, and complete randomization gives every other assignment probability zero. So its expectation is $\frac{1}{n_{1} - 1} \cdot (n_{1} - 1) S_{1}^{2} = S_{1}^{2}$ , where the $n_{1} - 1$ in the denominator is why $n_{1} \geq 2$ is needed.

For the control group, take $W_{i} = 1 - Z_{i}$ and $m = n - n_{1}$ . The two moment conditions have to be recomputed: $E [1 - Z_{i}] = 1 - n_{1} / n = (n - n_{1}) / n$ , and for $i \neq j$

E [(1 - Z_{i}) (1 - Z_{j})] = 1 - π_{i} - π_{j} + π_{i j} = 1 - \frac{2 n_{1}}{n} + \frac{n_{1} (n_{1} - 1)}{n (n - 1)} = \frac{(n - n_{1}) (n - n_{1} - 1)}{n (n - 1)},

which is the same formula with $n - n_{1}$ in place of $n_{1}$ . In other words the control group of a completely randomized experiment is itself a simple random sample without replacement, of size $n - n_{1}$ , and the rest of the argument is identical.∎

E[s₁²] = S₁²New tabOpens with 1,299 lines of library preamble — the code above is at the bottom of the editor.

theorem completeRandomization_expect_sampleVar_treated {n₁ : ℕ} (hle : n₁ ≤ n)
    (P : Population n) (hn₁ : 1 < n₁) (hlt : n₁ < n) :
    (completeRandomization n n₁ hle).expect (fun z => sampleVar P.y1 z) = P.S1sq := by
  have hn : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr (by omega)
  have hnm1 : (n : ℝ) - 1 ≠ 0 := by
    have : (2 : ℝ) ≤ (n : ℝ) := by exact_mod_cast (by omega : 2 ≤ n)
    linarith
  have hm : (n₁ : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr (by omega)
  have hm1 : (n₁ : ℝ) - 1 ≠ 0 := by
    have : (2 : ℝ) ≤ (n₁ : ℝ) := by exact_mod_cast hn₁
    linarith
  -- on the support the group size is `n₁`, so the statistic has the exchangeable shape
  have hsupp : ∀ z, (completeRandomization n n₁ hle).prob z ≠ 0 →
      sampleVar P.y1 z
        = ((∑ i, Z z i * (P.y1 i) ^ 2) - (∑ i, Z z i * P.y1 i) ^ 2 / (n₁ : ℝ))
            / ((n₁ : ℝ) - 1) := by
    intro z hz
    have hcard : #z = n₁ := card_eq_of_prob_ne_zero hz
    rw [sampleVar_eq P.y1 z (by rw [hcard]; exact hm), hcard,
      sum_mem_eq_sum_Z_mul z (fun i => (P.y1 i) ^ 2), sum_mem_eq_sum_Z_mul z P.y1]
  have hss : (completeRandomization n n₁ hle).expect
      (fun z => (∑ i, Z z i * (P.y1 i) ^ 2) - (∑ i, Z z i * P.y1 i) ^ 2 / (n₁ : ℝ))
      = ((n₁ : ℝ) - 1) * fpVar P.y1 :=
    (completeRandomization n n₁ hle).expect_ss_of_exchangeable P.y1 Z (n₁ : ℝ) hn hnm1 hm
      (fun z i => Z_mul_self z i)
      (fun i => completeRandomization_propensity i)
      (fun i j hij => completeRandomization_jointPropensity hij)
  rw [(completeRandomization n n₁ hle).expect_congr_of_support hsupp,
    (completeRandomization n n₁ hle).expect_div_const _ ((n₁ : ℝ) - 1), hss,
    Population.S1sq]
  field_simp

```
have hn : (n : ℝ) ≠ 0 := …; have hnm1 : (n : ℝ) - 1 ≠ 0 := …; have hm : (n₁ : ℝ) ≠ 0 := …; have hm1 : (n₁ : ℝ) - 1 ≠ 0 := …
```
Four non-vanishing facts, in the order the engine lemma and the final field_simp will want them. hm1 is the one that needs the full strength of 1 < n₁ rather than 0 < n₁: it is the denominator of the sample variance, and at a single treated unit that denominator is zero.
```
rw [sampleVar_eq P.y1 z (by rw [hcard]; exact hm), hcard, sum_mem_eq_sum_Z_mul z (fun i => (P.y1 i) ^ 2), sum_mem_eq_sum_Z_mul z P.y1]
```
The hsupp step: rewrite the sample variance, on an assignment of positive probability, into the shape the engine lemma expects. sampleVar_eq is the “sum of squares minus square of sum” identity over the group z, and it turns the left-hand side into the state shown below — its side condition (#z : ℝ) ≠ 0 is discharged by rewriting with hcard : #z = n₁ first. Then hcard replaces the two remaining #z, and the two sum_mem_eq_sum_Z_mul rewrites turn the sums over the treated group into indicator-weighted sums over everybody, once for y1 ^ 2 and once for y1.
Goal
```
⊢ (∑ i ∈ z, P.y1 i ^ 2 - (∑ i ∈ z, P.y1 i) ^ 2 / ↑(#z)) / (↑(#z) - 1) =
    (∑ i, Z z i * P.y1 i ^ 2 - (∑ i, Z z i * P.y1 i) ^ 2 / ↑n₁) / (↑n₁ - 1)
```
```
have hss : … := (completeRandomization n n₁ hle).expect_ss_of_exchangeable P.y1 Z (n₁ : ℝ) hn hnm1 hm (fun z i => Z_mul_self z i) (fun i => completeRandomization_propensity i) (fun i j hij => completeRandomization_jointPropensity hij)
```
The engine, instantiated. W := Z, m := n₁, v := P.y1, and the three hypotheses are chapter 6’s Z_mul_self : Z z i * Z z i = Z z i and chapter 7’s two counting lemmas — the propensity is $n_{1} / n$ and the joint propensity is $n_{1} (n_{1} - 1) / (n (n - 1))$ , which is exactly what h1 and h2 demand. That the counting lemmas match the hypotheses syntactically is not an accident: expect_ss_of_exchangeable was stated in their shape.

rw [(completeRandomization n n₁ hle).expect_congr_of_support hsupp, (completeRandomization n n₁ hle).expect_div_const _ ((n₁ : ℝ) - 1)]

Transfer the expectation onto the rewritten statistic, then pull the constant / (n₁ - 1) outside it. The remaining expectation is literally the left-hand side of hss.

Goal

⊢ ((completeRandomization n n₁ hle).expect fun z =>
      ∑ i, Z z i * P.y1 i ^ 2 - (∑ i, Z z i * P.y1 i) ^ 2 / ↑n₁) / (↑n₁ - 1) =
    P.S1sq

```
hss, Population.S1sq
```
Substitute the engine’s answer and unfold S1sq to fpVar P.y1. The goal is (n₁ - 1) * fpVar P.y1 / (n₁ - 1) = fpVar P.y1.
Goal
```
⊢ (↑n₁ - 1) * fpVar P.y1 / (↑n₁ - 1) = fpVar P.y1
```
```
field_simp
```
Cancel, using hm1 : (n₁ : ℝ) - 1 ≠ 0. If the theorem had said 1 ≤ n₁ instead of 1 < n₁ this is the line that would fail — and it would fail for a real reason, since with a single treated unit the sample variance is $0 / 0$ .

E[s₀²] = S₀²New tabOpens with 1,336 lines of library preamble — the code above is at the bottom of the editor.

theorem completeRandomization_expect_sampleVar_control {n₁ : ℕ} (hle : n₁ ≤ n)
    (P : Population n) (hpos : 0 < n₁) (hn₀ : 1 < n - n₁) :
    (completeRandomization n n₁ hle).expect (fun z => sampleVar P.y0 zᶜ) = P.S0sq := by
  set D := completeRandomization n n₁ hle with hD
  have hlt : n₁ < n := by omega
  have hn : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr (by omega)
  have hnm1 : (n : ℝ) - 1 ≠ 0 := by
    have : (2 : ℝ) ≤ (n : ℝ) := by exact_mod_cast (by omega : 2 ≤ n)
    linarith
  have hltR : (n₁ : ℝ) < (n : ℝ) := by exact_mod_cast hlt
  have hm : (n : ℝ) - (n₁ : ℝ) ≠ 0 := by linarith
  have hm1 : ((n : ℝ) - (n₁ : ℝ)) - 1 ≠ 0 := by
    have h2 : ((n₁ : ℝ) + 2) ≤ (n : ℝ) := by exact_mod_cast (by omega : n₁ + 2 ≤ n)
    linarith
  -- the two moment hypotheses for the control indicator
  have hone : ∀ i, D.expect (fun z => 1 - Z z i) = ((n : ℝ) - (n₁ : ℝ)) / (n : ℝ) := by
    intro i
    rw [Design.expect_one_sub_Z, hD, completeRandomization_propensity]
    field_simp
  have hpair : ∀ i j, i ≠ j → D.expect (fun z => (1 - Z z i) * (1 - Z z j))
      = ((n : ℝ) - (n₁ : ℝ)) * (((n : ℝ) - (n₁ : ℝ)) - 1) / ((n : ℝ) * ((n : ℝ) - 1)) := by
    intro i j hij
    have hexp : D.expect (fun z => (1 - Z z i) * (1 - Z z j))
        = 1 - D.propensity i - D.propensity j + D.jointPropensity i j := by
      rw [D.expect_congr (f := fun z => (1 - Z z i) * (1 - Z z j))
          (g := fun z => (1 - Z z i - Z z j) + Z z i * Z z j) (fun z => by ring),
        D.expect_add (fun z => 1 - Z z i - Z z j) (fun z => Z z i * Z z j),
        D.expect_sub (fun z => 1 - Z z i) (fun z => Z z j),
        Design.expect_one_sub_Z, Design.expect_Z]
      rfl
    rw [hexp, hD, completeRandomization_propensity, completeRandomization_propensity,
      completeRandomization_jointPropensity hij]
    field_simp
    ring
  -- on the support the control group has size `n - n₁`
  have hsupp : ∀ z, D.prob z ≠ 0 →
      sampleVar P.y0 zᶜ
        = ((∑ i, (1 - Z z i) * (P.y0 i) ^ 2)
            - (∑ i, (1 - Z z i) * P.y0 i) ^ 2 / ((n : ℝ) - (n₁ : ℝ)))
          / (((n : ℝ) - (n₁ : ℝ)) - 1) := 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 [sampleVar_eq P.y0 zᶜ (by rw [hc]; exact hm), hc,
      sum_compl_eq_sum_one_sub_Z_mul z (fun i => (P.y0 i) ^ 2),
      sum_compl_eq_sum_one_sub_Z_mul z P.y0]
  have hss : D.expect
      (fun z => (∑ i, (1 - Z z i) * (P.y0 i) ^ 2)
        - (∑ i, (1 - Z z i) * P.y0 i) ^ 2 / ((n : ℝ) - (n₁ : ℝ)))
      = (((n : ℝ) - (n₁ : ℝ)) - 1) * fpVar P.y0 :=
    D.expect_ss_of_exchangeable P.y0 (fun z i => 1 - Z z i) ((n : ℝ) - (n₁ : ℝ)) hn hnm1 hm
      (fun z i => by rw [one_sub_Z]; exact Z_mul_self zᶜ i) hone hpair
  rw [D.expect_congr_of_support hsupp,
    D.expect_div_const _ (((n : ℝ) - (n₁ : ℝ)) - 1), hss, Population.S0sq]
  field_simp

The control proof has the same five moves; only the two moment conditions are new, because 1 - Z is not a primitive of the library the way Z is.

```
set D := completeRandomization n n₁ hle with hD
```
Abbreviate the design. set introduces a local definition D and replaces every occurrence in the goal, keeping the defining equation as hD so that design-specific lemmas can still be applied by rewriting with it. The treated proof did without this; here the design appears often enough to be worth naming.
```
have hm1 : ((n : ℝ) - (n₁ : ℝ)) - 1 ≠ 0 := by have h2 : ((n₁ : ℝ) + 2) ≤ (n : ℝ) := by exact_mod_cast (by omega : n₁ + 2 ≤ n); linarith
```
The control-side counterpart of hm1. hn₀ : 1 < n - n₁ is a statement about natural subtraction, so it cannot be cast directly; omega converts it into the equivalent n₁ + 2 ≤ n, which casts cleanly, and linarith finishes in $R$ . Mixing the two subtractions is the mistake this two-step dance avoids.
```
rw [Design.expect_one_sub_Z, hD, completeRandomization_propensity]; field_simp
```
The first moment condition, hone. expect_one_sub_Z : E[1 - Z i] = 1 - π i is chapter 7’s; the propensity lemma makes that 1 - n₁ / n; and field_simp checks that this equals (n - n₁) / n.
```
rw [D.expect_congr (f := fun z => (1 - Z z i) * (1 - Z z j)) (g := fun z => (1 - Z z i - Z z j) + Z z i * Z z j) (fun z => by ring), D.expect_add …, D.expect_sub …, Design.expect_one_sub_Z, Design.expect_Z]
```
The second moment condition, hpair, and its have hexp. Expanding the product pointwise with ring under expect_congr turns $(1 - Z_{i}) (1 - Z_{j})$ into $(1 - Z_{i} - Z_{j}) + Z_{i} Z_{j}$ ; then expect_add and expect_sub split it and the two indicator-expectation lemmas evaluate the first three pieces. The rfl that follows identifies the leftover E[Z i * Z j] with jointPropensity i j, which is its definition. The result is $E [(1 - Z_{i}) (1 - Z_{j})] = 1 - π_{i} - π_{j} + π_{i j}$ , the inclusion–exclusion line of the informal proof.

rw [hexp, hD, completeRandomization_propensity, completeRandomization_propensity, completeRandomization_jointPropensity hij]; field_simp; ring

Substitute the three chapter-7 values into that identity and verify the resulting rational-function identity. The goal shown is the state before field_simp: it is the middle display of the informal proof, and field_simp; ring is the “which is” that follows it.

Goal

⊢ 1 - ↑n₁ / ↑n - ↑n₁ / ↑n + ↑n₁ * (↑n₁ - 1) / (↑n * (↑n - 1)) =
    (↑n - ↑n₁) * (↑n - ↑n₁ - 1) / (↑n * (↑n - 1))

```
rw [sampleVar_eq P.y0 zᶜ (by rw [hc]; exact hm), hc, sum_compl_eq_sum_one_sub_Z_mul z (fun i => (P.y0 i) ^ 2), sum_compl_eq_sum_one_sub_Z_mul z P.y0]
```
The support step, hsupp, with zᶜ in place of z: the same “sum of squares minus square of sum” rewrite, the control-group size hc : (#zᶜ : ℝ) = n - n₁ obtained from card_compl and Nat.cast_sub as before, and sum_compl_eq_sum_one_sub_Z_mul in place of sum_mem_eq_sum_Z_mul.
```
D.expect_ss_of_exchangeable P.y0 (fun z i => 1 - Z z i) ((n : ℝ) - (n₁ : ℝ)) hn hnm1 hm (fun z i => by rw [one_sub_Z]; exact Z_mul_self zᶜ i) hone hpair
```
The engine again, with $W = 1 - Z$ and $m = n - n_{1}$ . Idempotence of the control indicator is proved by one_sub_Z : 1 - Z z i = Z zᶜ i — the control indicator is the treated indicator of the complementary assignment — followed by Z_mul_self at zᶜ. The last two arguments are the moment conditions just established.
```
rw [D.expect_congr_of_support hsupp, D.expect_div_const _ (((n : ℝ) - (n₁ : ℝ)) - 1), hss, Population.S0sq]; field_simp
```
Identical to the treated ending: transfer onto the support, pull the constant denominator out, substitute the engine’s answer, unfold S0sq, cancel.

Adding the two:

Proposition. Under complete randomization with $1 < n_{1}$ and $1 < n_{0}$ , $E [\hat{V}] = \frac{S_{1}^{2}}{n_{1}} + \frac{S_{0}^{2}}{n_{0}}$ , where $\hat{V} = \frac{s_{1}^{2}}{n_{1}} + \frac{s_{0}^{2}}{n_{0}}$ is the reported variance.

Proof.

On the support the two group sizes are the constants $n_{1}$ and $n_{0}$ , so $\hat{V}$ is a fixed linear combination of the two sample variances; by linearity of expectation and the previous theorem its expectation is $S_{1}^{2} / n_{1} + S_{0}^{2} / n_{0}$ .∎

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

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₀]

```
rw [neymanVarEst, hcard, hc]
```
The hsupp step. neymanVarEst divides by the realized group sizes #z and #zᶜ; on the support those are the constants n₁ and (n : ℝ) - (n₁ : ℝ), and hcard and hc put them in. Until this rewrite the estimator is not a linear combination of anything, because its coefficients are random.
```
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₁ : ℝ)))]
```
Transfer the expectation onto that constant-coefficient version, then split it with linearity. The two funs spell out the summands, which Lean cannot infer from the goal.

D.expect_div_const _ (n₁ : ℝ), D.expect_div_const _ ((n : ℝ) - (n₁ : ℝ))

Pull each constant denominator outside its expectation. What remains are exactly the two theorems just proved.

Goal

⊢ D.expect (sampleVar P.y1) / ↑n₁ + (D.expect fun z => sampleVar P.y0 zᶜ) / (↑n - ↑n₁) =
    P.S1sq / ↑n₁ + P.S0sq / (↑n - ↑n₁)

```
hD, completeRandomization_expect_sampleVar_treated hle P hn₁ hlt, completeRandomization_expect_sampleVar_control hle P (by omega) hn₀
```
Substitute them. hD first, to turn D back into the design those theorems are stated about — the price of the set above. The (by omega) supplies 0 < n₁ for the control theorem, which is derivable from hn₁ : 1 < n₁.

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.

Proposition. $E [\hat{V}] = Var ({\hat{τ}}_{DiM}) + \frac{S_{τ}^{2}}{n}$ .

Proof.

Subtract the two formulas: $\frac{S_{1}^{2}}{n_{1}} + \frac{S_{0}^{2}}{n_{0}}$ is what the estimator estimates on average, and Neyman’s theorem says the true variance is that same quantity minus $S_{τ}^{2} / n$ .∎

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

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

```
rw [completeRandomization_expect_neymanVarEst hle P hn₁ hn₀, neyman_variance hle P (by omega) (by omega)]
```
Both formulas at once: the left-hand side becomes $S_{1}^{2} / n_{1} + S_{0}^{2} / n_{0}$ and the var on the right becomes Neyman’s three terms. The two by omega arguments derive 0 < n₁ and n₁ < n from the estimator’s stronger hypotheses 1 < n₁ and 1 < n - n₁ — note that the statement keeps the stronger ones, because the left-hand side needs them.
Goal
```
⊢ P.S1sq / ↑n₁ + P.S0sq / (↑n - ↑n₁) =
    P.S1sq / ↑n₁ + P.S0sq / (↑n - ↑n₁) - P.Stausq / ↑n + P.Stausq / ↑n
```
```
ring
```
The goal is $X = X - c + c$ . Nothing is being proved here beyond the cancellation, which is the point: the bias is not estimated or bounded, it is the literal difference between two closed-form expressions.

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)

On average, the reported variance is at least the true variance of difference in means: $Var ({\hat{τ}}_{DiM}) \leq E [\hat{V}]$ . Note that no side condition is needed beyond the two the estimator already required.

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

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

Hint

One rewrite and one linarith: neymanVarEst_bias turns the goal into Var ≤ Var + Stausq / n, and fpVar_nonneg (which needs no hypothesis) with div_nonneg says the added term is a non-negative real divided by a non-negative real.

Show solution

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

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

Here “conservative” refers to the precise statement

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 1,550 lines of library preamble — the code above is at the bottom of the editor.

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

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

```
obtain ⟨h1, h0, _⟩ := exP_variances
```
exP_variances is a conjunction of three equations; obtain with an anonymous constructor pattern destructures it into named hypotheses, and the _ throws away the one about $S_{τ}^{2}$ , which this theorem does not need.
```
rw [completeRandomization_expect_neymanVarEst (by norm_num) exP (by norm_num) (by norm_num), h1, h0]; norm_num
```
Quote the general formula at $n = 4$ , $n_{1} = 2$ , discharging its three side conditions ( $2 \leq 4$ , $1 < 2$ , $1 < 4 - 2$ ) by norm_num, then substitute the two variances. What is left is $1 / 2 + 1 / 2 = 1$ .
```
rw [exP_expect_neymanVarEst, exP_var_DiM]; norm_num
```
The gap, as a subtraction of two numbers each already computed: $1 - 2 / 3 = 1 / 3$ . Nothing here is a new mathematical fact; the value of the theorem is that the three numbers are pinned by the kernel rather than by arithmetic in the margin.

$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.

Theorem. Under complete randomization with $1 < n_{1}$ and $1 < n_{0}$ ,

E [\hat{V}] = Var ({\hat{τ}}_{DiM}) ⟺ τ_{i} = τ_{j} for all i, j,

that is, if and only if the treatment effect is the same for every unit.

Proof.

By the previous proposition the two sides differ by exactly $S_{τ}^{2} / n$ , so they are equal if and only if $S_{τ}^{2} = 0$ (using $n > 0$ ). A finite-population variance vanishes if and only if its vector is constant, and the vector here is the vector of individual treatment effects.∎

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

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]

```
have hn : 1 < n := by omega; have hnR : (0 : ℝ) < (n : ℝ) := …
```
Two facts about $n$ : at least two units, from 1 < n₁ and 1 < n - n₁, and its positivity in $R$ . The first is the hypothesis of fpVar_eq_zero_iff, the second is what rules out dividing by zero in the forward direction.

rw [neymanVarEst_bias hle P hn₁ hn₀]

Replace the expectation by “the variance plus the gap” — the previous proposition, used as a rewrite. The left-hand side of the ↔ is now visibly of the form X + c = X, and no expectation appears in it any more.

Goal

⊢ (completeRandomization n n₁ hle).var (DiM P) + P.Stausq / ↑n =
      (completeRandomization n n₁ hle).var (DiM P) ↔
    ∀ (i j : Fin n), P.effect i = P.effect j

```
constructor
```
Split the ↔ into its two implications, mp and mpr. From here the two branches are separate proofs, marked by the two · bullets.
```
intro h; have hq : P.Stausq / (n : ℝ) = 0 := by linarith
```
Forward direction. h : Var + Stausq / n = Var is a linear equation between reals, so linarith extracts Stausq / n = 0 from it directly — no cancellation lemma needed, because linarith works over ordered fields and treats Stausq / n as an atom.
```
exact (div_eq_zero_iff.mp hq).resolve_right hnR.ne'
```
A quotient is zero when the numerator is zero or the denominator is; div_eq_zero_iff gives that disjunction and resolve_right discards the second case using hnR.ne' : (n : ℝ) ≠ 0. This is the step that needs $n > 0$ , and it is the reason the hypothesis was derived at the top.
```
exact (fpVar_eq_zero_iff hn P.effect).mp hzero
```
Conclude. fpVar_eq_zero_iff : 1 < n → (fpVar v = 0 ↔ ∀ i j, v i = v j) applied to the effect vector turns Stausq = 0 — which is fpVar P.effect = 0 by definition — into constant effects.
```
intro h; have hzero : P.Stausq = 0 := (fpVar_eq_zero_iff hn P.effect).mpr h; rw [hzero, zero_div, add_zero]
```
Backward direction, the same lemma in the other direction: constant effects make the variance of the effect vector zero, and then zero_div and add_zero collapse Var + 0 / n to Var. Both branches go through the same ↔, which is why the theorem is an ↔ rather than two separate results.

Constant effects is the assumption “the treatment effect is the same for every unit”. 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.

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

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

Hint

Get P.y1 = fun i => P.y0 i + c by funext and linarith from hc, then fpVar_add_const gives $S_{1}^{2} = S_{0}^{2}$ and fpVar_eq_zero_iff gives $S_{τ}^{2} = 0$ ; rewrite both into neyman_variance and field_simp; ring finishes. 1 < n has to be derived from 0 < n₁ and n₁ < n by omega before fpVar_eq_zero_iff will apply, and field_simp wants (n₁ : ℝ) ≠ 0 and (n : ℝ) - n₁ ≠ 0 as haves.

Show solution

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

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.

Proposition. For any design,

Var ({\hat{τ}}_{HT}) = \sum_{i} \sum_{j} b_{i} b_{j} (π_{i j} - π_{i} π_{j}), b_{i} = \frac{1}{n} (\frac{y_{i} (1)}{π_{i}} + \frac{y_{i} (0)}{1 - π_{i}}) .

Proof.

Regroup the estimator: collecting the terms in $Z_{i}$ ,

{\hat{τ}}_{HT} (z) = \sum_{i} b_{i} Z_{i} - \frac{1}{n} \sum_{i} \frac{y_{i} (0)}{1 - π_{i}},

with $b_{i}$ as displayed, and the subtracted quantity does not depend on the assignment. As before, adding a constant does not change a variance, so the result is the general linear-statistic formula applied to the coefficients $b$ . The regrouping is valid on every assignment, so no support argument is needed.∎

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

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))

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]

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]

```
have h1 : HT D P z = ∑ i, (1 / (n : ℝ)) * (…) := Finset.mul_sum _ _ _
```
HT_eq_linear first. The definition of HT has the 1 / n outside the sum; Finset.mul_sum moves it in, so that each summand can be regrouped independently. The have is stated separately only because the calc block below wants it as a named first step.
```
_ = ∑ 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
```
The regrouping, one summand at a time. Unfolding htCoef makes the goal a polynomial identity in Z z i and the four constants, and ring verifies it — the same move as in DiM_eq_linear_of_card, and again it is where the control outcome acquires its + sign.

_ = (∑ i, htCoef D P i * Z z i) + ∑ i, -(…) := Finset.sum_add_distrib

Split into the linear part and the constant part.

```
rw [Finset.sum_neg_distrib, ← Finset.mul_sum]
```
Tidy the constant part into the shape the statement gives it: pull the negation out of the sum, then pull the 1 / n back out. The calc chain now ends at the statement.

have h := D.var_congr_of_support (f := HT D P) (g := fun z => (∑ i, htCoef D P i * Z z i) + -(…)) (fun z _ => HT_eq_linear D P z)

var_HT proper. The support lemma is being used purely as a rewriting device: its third argument should be a proof valid only where D.prob z ≠ 0, and what is supplied ignores that hypothesis entirely. The goal after rw [h] is the variance of the linear-plus-constant form.

Goal

⊢ (D.var fun z => ∑ i, htCoef D P i * Z z i + -(1 / ↑n * ∑ i, P.y0 i / (1 - D.propensity i))) =
    ∑ i, ∑ j, htCoef D P i * htCoef D P j * (D.jointPropensity i j - D.propensity i * D.propensity j)

```
rw [h, D.var_add_const, D.var_linear]
```
Two lemmas and the theorem is done: var_add_const deletes the constant tail and var_linear — the one formula this whole chapter is built on — evaluates what is left.

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)

Under the Bernoulli design the Horvitz–Thompson variance collapses to a single sum, $p (1 - p) \sum_{i} b_{i}^{2}$ , over the same coefficients $b_{i}$ as the proposition above, because distinct units are assigned independently.

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

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

Hint

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

Show solution

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

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

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

noncomputable def varEstBad {n : ℕ} (P : Population n) (n₁ : ℕ)
    (_z : Assignment n) : ℝ :=
  P.S1sq / (n₁ : ℝ) + P.S0sq / ((n : ℝ) - (n₁ : ℝ))

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.