Unbiasedness — MrCLean

Everything so far has been scaffolding. The population is fixed, the potential outcomes are fixed, the design is an explicit pmf over $2^{n}$ assignments, and expectation is a finite sum. This chapter spends that machinery on the two results the whole exercise was for:

the Horvitz–Thompson estimator is unbiased for $τ$ under any design whose propensity scores are strictly inside $(0, 1)$ ;
the difference-in-means estimator is unbiased for $τ$ under complete randomization.

The second falls out of the first in about fifteen lines, by an argument (“these two estimators are the same function on the support of the design”) that is worth more than either theorem.

The two estimators

Horvitz–Thompson reweights each observed outcome by the reciprocal of the probability that it was observed:

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

The Horvitz–Thompson estimatorNew tabOpens with 849 lines of library preamble — the code above is at the bottom of the editor.

/-- The **Horvitz--Thompson** (inverse-probability-weighted) estimator of the ATE:

`HT = (1/n) ∑ i, ( Z i · y1 i / π i  -  (1 - Z i) · y0 i / (1 - π i) )`.

Each observed outcome is reweighted by the reciprocal of the probability that it was
observed, so that units the design was unlikely to reveal count for more.  Note that the
summand only involves the potential outcome that the assignment actually reveals: the
indicator `Z i` is `0` exactly when `y1 i` is unobserved, and vice versa, so `HT` is a
genuine statistic computable from data (see `HT_eq_observed`). -/
noncomputable def HT (D : Design n) (P : Population n) (z : Assignment n) : ℝ :=
  (1 / (n : ℝ)) * ∑ i, (Z z i * P.y1 i / D.propensity i
    - (1 - Z z i) * P.y0 i / (1 - D.propensity i))

Two things to notice in the Lean. First, HT takes the design D as an argument, because $π_{i}$ is a property of the design and not of the data; under a design-based view that is not a nuisance parameter to be estimated, it is a number the experimenter chose. Second, HT is a function of z, the assignment: an estimator, formally, is a map from assignments to reals, and “unbiased” will be a statement about the expectation of that map.

It is not obvious from the definition that HT is computable from data — P.y1 i appears for every i, including units that were never treated. It is, because the indicator kills the unobserved half of every summand:

HT only ever touches the revealed outcomeNew tabOpens with 862 lines of library preamble — the code above is at the bottom of the editor.

/-- `HT` really is computable from observed data: every potential outcome appearing in it with
a nonzero coefficient is the one the assignment reveals. -/
theorem HT_eq_observed (D : Design n) (P : Population n) (z : Assignment n) :
    HT D P z = (1 / (n : ℝ)) * ∑ i, (Z z i * Yobs P z i / D.propensity i
      - (1 - Z z i) * Yobs P z i / (1 - D.propensity i)) := by
  unfold HT
  congr 1
  refine Finset.sum_congr rfl fun i _ => ?_
  by_cases hi : i ∈ z <;> simp [Yobs_eq_ite, Z, hi]

Difference in means is the estimator everyone actually runs:

The difference-in-means estimatorNew tabOpens with 918 lines of library preamble — the code above is at the bottom of the editor.

/-- The **difference-in-means** estimator: the average observed outcome among the treated
minus the average observed outcome among the controls.

Because a treated unit reveals `y1` and a control unit reveals `y0`, this is
`(∑ i ∈ z, y1 i) / #z - (∑ i ∈ zᶜ, y0 i) / #zᶜ`, which is the form we take as the
definition.  `DiM_eq_observed` confirms it is what a data analyst would actually compute.

Lean's convention `x / 0 = 0` handles the degenerate assignments (everyone treated, or nobody
treated) without any side conditions; those assignments have probability `0` under the
designs we care about anyway. -/
noncomputable def DiM (P : Population n) (z : Assignment n) : ℝ :=
  (∑ i ∈ z, P.y1 i) / (#z : ℝ) - (∑ i ∈ zᶜ, P.y0 i) / (#zᶜ : ℝ)

#z is the number of treated units and zᶜ is the control group — recall from chapter 6 that an assignment is the treated set, so the control group is free.

LemmaLean's division convention

Division by zero is not an error, an exception, or an undefined term: it is zero. div_zero is the lemma that says so, and it is used, silently, every time a formal statement about a mean or a ratio omits a positivity hypothesis.

New tab

example (x : ℝ) : x / 0 = 0 := div_zero x

Horvitz–Thompson is unbiased

Informally. Fix a finite population and a design. If every unit has $0 < π_{i} < 1$ , then $E [{\hat{τ}}_{HT}] = τ$ , where the expectation is over the randomization and nothing else. No model for the outcomes, no asymptotics, no sampling from a superpopulation.

Formally.

HT_unbiasedNew tabOpens with 890 lines of library preamble — the code above is at the bottom of the editor.

/-- **The Horvitz--Thompson estimator is unbiased for the ATE.**

For *any* design in which every unit has a propensity strictly between `0` and `1`,
`E[HT] = τ`.  The proof is linearity of expectation plus the one-unit cancellation
`expect_HT_term`; no structure of the design is used beyond positivity, which is why the
result covers Bernoulli, complete randomisation, stratified and clustered designs alike. -/
theorem HT_unbiased (D : Design n) (P : Population n)
    (hlo : ∀ i, 0 < D.propensity i) (hhi : ∀ i, D.propensity i < 1) :
    D.expect (HT D P) = P.tau := by
  have hexpand : D.expect (HT D P)
      = (1 / (n : ℝ)) * D.expect (fun z => ∑ i, (Z z i * P.y1 i / D.propensity i
          - (1 - Z z i) * P.y0 i / (1 - D.propensity i))) :=
    D.expect_const_mul _ (1 / (n : ℝ))
  rw [hexpand, D.expect_sum Finset.univ
    (fun i z => Z z i * P.y1 i / D.propensity i - (1 - Z z i) * P.y0 i / (1 - D.propensity i)),
    Finset.sum_congr rfl fun i _ => expect_HT_term D P i (hlo i) (hhi i), P.tau_eq]
  ring

Read the statement slowly, because every piece of it is load-bearing.

D : Design n is universally quantified and otherwise unconstrained. There is no independence assumption, no “the units are exchangeable”, no fixed number of treated units. Bernoulli, complete randomization, stratified, clustered, and any hand-built pmf you like all satisfy this.
P : Population n is universally quantified too. The theorem holds for every configuration of potential outcomes, which is what “design-based” means: the randomness is entirely in D.
D.expect (HT D P) is $E [{\hat{τ}}_{HT}]$ , i.e. ∑ z, D.prob z * HT D P z. It is a sum over all $2^{n}$ assignments.
P.tau is the population ATE, mean P.effect. It is a fixed real number, not a random variable.
hlo and hhi are the positivity conditions, stated pointwise over units.

What `0 < π i` is doing there

Suppose you drop it. The theorem does not become vacuous — it becomes false, and the reason is div_zero. Showing that takes a population and a design, so here are both. The population is as plain as it can be and still have an effect to miss: two units, $y_{1 i} = 1$ , $y_{0 i} = 0$ , hence $τ = 1$ . Every counterexample in this chapter runs against it.

The population used by every counterexample belowNew tabOpens with 1,050 lines of library preamble — the code above is at the bottom of the editor.

/-- Two units, `y1 ≡ 1` and `y0 ≡ 0`.  Everything below is a counterexample against this
population, whose average treatment effect is `1`. -/
def twoUnits : Population 2 where
  y1 := fun _ => 1
  y0 := fun _ => 0

@[simp] theorem twoUnits_y1 : twoUnits.y1 = fun _ => (1 : ℝ) := rfl

@[simp] theorem twoUnits_y0 : twoUnits.y0 = fun _ => (0 : ℝ) := rfl

theorem twoUnits_mean_y1 : mean twoUnits.y1 = 1 := by
  rw [twoUnits_y1]; exact mean_const (by norm_num) 1

theorem twoUnits_mean_y0 : mean twoUnits.y0 = 0 := by
  rw [twoUnits_y0]; exact mean_const (by norm_num) 0

theorem twoUnits_tau : twoUnits.tau = 1 := by
  rw [Population.tau_eq_mean_sub_mean, twoUnits_mean_y1, twoUnits_mean_y0, sub_zero]

And here is a legitimate design in which nobody is ever treated. It is a point mass on the empty assignment, it satisfies both Design axioms — the probabilities are non-negative and they sum to one — and every propensity score is $0$ :

A design that never treats anybodyNew tabOpens with 1,069 lines of library preamble — the code above is at the bottom of the editor.

/-- The design that assigns nobody to treatment: a point mass on the empty assignment.

It is a perfectly legal `Design` — the probabilities are nonnegative and sum to one — and
every propensity score is `0`, so it is exactly the design `HT_unbiased` excludes. -/
noncomputable def neverTreat (n : ℕ) : Design n where
  prob z := if z = ∅ then 1 else 0
  nonneg z := by by_cases h : z = ∅ <;> simp [h]
  sum_one := by simp

/-- Under a point mass, expectation is evaluation. -/
theorem expect_neverTreat (f : Assignment n → ℝ) : (neverTreat n).expect f = f ∅ := by
  simp [Design.expect, neverTreat]

/-- Every propensity score of `neverTreat` is `0`. -/
@[simp] theorem neverTreat_propensity (i : Fin n) : (neverTreat n).propensity i = 0 := by
  rw [Design.propensity, expect_neverTreat]
  simp [Z]

Under this design the treated half of each Horvitz–Thompson summand is 0 * y1 i / 0. In ordinary mathematics that is undefined and the discussion stops; in Lean it is 0, so the estimator quietly evaluates to minus the control mean — and against twoUnits its expectation is $0$ while $τ = 1$ :

Without positivity the theorem is false, not vacuousNew tabOpens with 1,089 lines of library preamble — the code above is at the bottom of the editor.

/-- With every propensity equal to `0`, the treated half of every Horvitz--Thompson summand
is `0 * y1 i / 0`, which Lean evaluates to `0` rather than to nonsense.  What survives is the
control half, so `HT` collapses to `-mean y0`.

This is the precise sense in which `HT_unbiased` *fails* at `π i = 0` instead of being
vacuously true: both sides of the equation are perfectly well-defined real numbers, and they
are different. -/
theorem HT_neverTreat_empty (P : Population n) :
    HT (neverTreat n) P ∅ = -mean P.y0 := by
  have hz : ∀ i : Fin n, Z (∅ : Assignment n) i = 0 := fun i => by simp [Z]
  unfold HT
  simp only [hz, neverTreat_propensity, zero_mul, zero_div, sub_zero, div_one, one_mul,
    zero_sub, Finset.sum_neg_distrib]
  simp only [mean]
  ring

/-- **`HT_unbiased` is false without `0 < π i`.**  Against `twoUnits` the expectation of the
Horvitz--Thompson estimator under `neverTreat` is `0`, while `τ = 1`. -/
theorem HT_not_unbiased_at_zero_propensity :
    (neverTreat 2).expect (HT (neverTreat 2) twoUnits) ≠ twoUnits.tau := by
  rw [expect_neverTreat, HT_neverTreat_empty, twoUnits_mean_y0, twoUnits_tau]
  norm_num

≠ is how the falsity is recorded: HT_not_unbiased_at_zero_propensity is a proof that the two sides are different real numbers, which is a stronger and more useful thing to have on file than the absence of a proof that they are equal.

hhi : ∀ i, D.propensity i < 1 plays the mirror-image role. At $π_{i} = 1$ the unit is always treated, 1 - Z z i is always 0, and the control half of the summand is 0 * y0 i / 0 = 0: the estimator behaves as though $y_{0 i} = 0$ , which it generally is not. Note that hhi is not vacuous either — Design.propensity_le_one says $π_{i} \leq 1$ always, so < 1 excludes exactly the degenerate designs and nothing else. (An autoformalization that wrote 1 < D.propensity i would be excluding everything; see the last section.)

The proof, lemma by lemma

The whole argument is one line of statistics — linearity of expectation, then cancel the weight — and the Lean tracks it step for step.

Step 1: one unit at a time. The only real content is that a single summand has expectation equal to that unit’s individual treatment effect:

expect_HT_termNew tabOpens with 872 lines of library preamble — the code above is at the bottom of the editor.

/-- The contribution of a single unit to the expectation of `HT` is exactly that unit's
individual treatment effect `y1 i - y0 i`.

This is the heart of the Horvitz--Thompson argument: `E[Z i] = π i` cancels the weight
`1 / π i`, and `E[1 - Z i] = 1 - π i` cancels the weight `1 / (1 - π i)`.  Positivity
(`0 < π i < 1`) is exactly what makes both cancellations legal. -/
theorem expect_HT_term (D : Design n) (P : Population n) (i : Fin n)
    (hlo : 0 < D.propensity i) (hhi : D.propensity i < 1) :
    D.expect (fun z => Z z i * P.y1 i / D.propensity i
      - (1 - Z z i) * P.y0 i / (1 - D.propensity i)) = P.y1 i - P.y0 i := by
  have hne : D.propensity i ≠ 0 := ne_of_gt hlo
  have hne' : 1 - D.propensity i ≠ 0 := sub_ne_zero_of_ne (ne_of_gt hhi)
  rw [D.expect_congr (g := fun z => Z z i * (P.y1 i / D.propensity i)
      - (1 - Z z i) * (P.y0 i / (1 - D.propensity i)))
      fun z => by rw [mul_div_assoc, mul_div_assoc],
    D.expect_sub, D.expect_mul_const, D.expect_mul_const, Design.expect_Z,
    Design.expect_one_sub_Z, mul_div_cancel₀ _ hne, mul_div_cancel₀ _ hne']

The proof is a single rw (rewrite: replace occurrences of a lemma’s left-hand side by its right-hand side in the goal) whose bracketed list does five things, in order:

D.expect_congr … fun z => by rw [mul_div_assoc, mul_div_assoc] re-associates Z z i * y1 i / π i into Z z i * (y1 i / π i). expect_congr is how you rewrite under the binder: two statistics that agree at every z have the same expectation. You cannot rw inside a fun z => … directly, so this lemma is the standard workaround.
D.expect_sub splits the expectation of a difference. This is linearity, and in this library its entire proof is simp only [expect, mul_sub, Finset.sum_sub_distrib] — there is no measure theory underneath it, only the distributive law and a rearranged finite sum.
D.expect_mul_const, twice, pulls the constants y1 i / π i and y0 i / (1 - π i) out of the expectation.
Design.expect_Z and Design.expect_one_sub_Z are the indicator expectation step: $E [Z_{i}] = π_{i}$ and $E [1 - Z_{i}] = 1 - π_{i}$ . The first is rfl — true by definition, because propensity is defined as E[Z i]. That is a design decision worth copying: define the quantity you will need as the expectation you will meet.
mul_div_cancel₀ _ hne and mul_div_cancel₀ _ hne' are the cancellations $π_{i} \cdot (y_{1 i} / π_{i}) = y_{1 i}$ . This is the only place the positivity hypotheses are used, and they enter as hne : D.propensity i ≠ 0 and hne' : 1 - D.propensity i ≠ 0.

Step 2: sum over units. With the per-unit fact in hand, the theorem is bookkeeping. Reading the proof of HT_unbiased above:

have hexpand : … := D.expect_const_mul _ (1 / (n : ℝ)) pulls the leading $1 / n$ outside the expectation. have introduces a named intermediate fact; here its statement is written out in full so a reader can see the shape of the expression at that point.
D.expect_sum Finset.univ (fun i z => …) exchanges $E$ and $\sum_{i}$ . This is Fubini for finite sums, and it is proved by Finset.sum_comm, i.e. by swapping the order of a double sum. Nothing deeper is happening.
Finset.sum_congr rfl fun i _ => expect_HT_term D P i (hlo i) (hhi i) applies step 1 inside the sum, once per unit. rfl says the index set is unchanged; the function argument supplies, for each i, a proof that its summand equals y1 i - y0 i.
P.tau_eq unfolds $τ$ to (∑ i, (P.y1 i - P.y0 i)) / n.
ring closes the remaining goal (1/n) * ∑ i, (y1 i - y0 i) = (∑ i, (y1 i - y0 i)) / n. ring proves any identity that holds in a commutative ring, treating division as multiplication by an inverse, so it does not need to know that $n \neq 0$ .

That is the whole theorem. The reason it is short is that the work was done earlier: Design.expect_sum, Design.expect_const_mul and Design.expect_Z are the entire probabilistic content, and each is a two-line simp only.

For the Bernoulli design the propensities are constant and the corollary is immediate:

HT_unbiased_bernoulliNew tabOpens with 993 lines of library preamble — the code above is at the bottom of the editor.

/-- **The Horvitz--Thompson estimator is unbiased under the Bernoulli design** with any
`0 < p < 1`.  A direct corollary of `HT_unbiased` and `bernoulli_propensity`. -/
theorem HT_unbiased_bernoulli {p : ℝ} (hp0 : 0 ≤ p) (hp1 : p ≤ 1) (P : Population n)
    (hlo : 0 < p) (hhi : p < 1) :
    (bernoulliDesign n p hp0 hp1).expect (HT (bernoulliDesign n p hp0 hp1) P) = P.tau :=
  HT_unbiased _ P (fun i => by rw [bernoulli_propensity]; exact hlo)
    (fun i => by rw [bernoulli_propensity]; exact hhi)

Difference in means, under complete randomization

Informally. Under complete randomization with $n_{1}$ treated units, $0 < n_{1} < n$ , the difference in means is unbiased for $τ$ . This is Neyman (1923).

Formally.

DiM_unbiased_completeRandomizationNew tabOpens with 970 lines of library preamble — the code above is at the bottom of the editor.

/-- **The difference-in-means estimator is unbiased for the ATE under complete randomisation.**

For `0 < n₁ < n`, `E[DiM] = τ`.  This is the classical Neyman (1923) result, and it is proved
here by showing that difference-in-means and Horvitz--Thompson agree on every assignment the
design can produce, and then quoting `HT_unbiased`. -/
theorem DiM_unbiased_completeRandomization {n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n)
    (hpos : 0 < n₁) (hlt : n₁ < n) :
    (completeRandomization n n₁ hle).expect (DiM P) = P.tau := by
  have hnR : (0 : ℝ) < (n : ℝ) := by exact_mod_cast Nat.lt_of_lt_of_le hpos hle
  have hn1R : (0 : ℝ) < (n₁ : ℝ) := by exact_mod_cast hpos
  have hltR : (n₁ : ℝ) < (n : ℝ) := by exact_mod_cast hlt
  have hsupport : ∀ z, (completeRandomization n n₁ hle).prob z ≠ 0 →
      DiM P z = HT (completeRandomization n n₁ hle) P z :=
    fun z hz => DiM_eq_HT_of_card hle P hpos hlt (card_eq_of_prob_ne_zero hz)
  rw [(completeRandomization n n₁ hle).expect_congr_of_support hsupport]
  refine HT_unbiased _ P (fun i => ?_) (fun i => ?_)
  · rw [completeRandomization_propensity]
    exact div_pos hn1R hnR
  · rw [completeRandomization_propensity]
    exact (div_lt_one hnR).mpr hltR

0 < n₁ and n₁ < n are needed for the obvious reason. At $n_{1} = 0$ there is exactly one subset of size $0$ , so the design is a point mass on the empty assignment — the same pmf as neverTreat above, written differently — and DiM returns -mean y0. At $n_{1} = n$ it is a point mass on everybody, and DiM returns mean y1. Neither is $τ$ except by accident. hle : n₁ ≤ n is not a mathematical hypothesis at all: completeRandomization n n₁ hle cannot be written down without it, because a uniform distribution over the size- $n_{1}$ subsets of an $n$ -element set needs there to be some.

Equal on the support

The proof does not touch a single sum over assignments. It observes that under complete randomization DiM and HT are the same function wherever the design puts any mass, and then quotes HT_unbiased. The tool is:

Design.expect_congr_of_supportNew tabOpens with 274 lines of library preamble — the code above is at the bottom of the editor.

/-- Two statistics that agree on the **support** of the design have the same expectation.

This is how we transfer a computation from "all assignments" to "the assignments the design
can actually produce" -- for instance, difference-in-means only behaves well on assignments
with the right number of treated units, and under complete randomisation those are exactly
the assignments of positive probability. -/
theorem expect_congr_of_support {f g : Assignment n → ℝ}
    (h : ∀ z, D.prob z ≠ 0 → f z = g z) : D.expect f = D.expect g := by
  refine Finset.sum_congr rfl fun z _ => ?_
  by_cases hz : D.prob z = 0
  · rw [hz, zero_mul, zero_mul]
  · rw [h z hz]

In words: to compute $E [f]$ you may replace $f$ by any $g$ that agrees with it on the support. The proof is four lines: Finset.sum_congr reduces to one assignment at a time, by_cases splits on whether D.prob z = 0, and in the zero branch both terms of the sum are 0 * … and no hypothesis is needed.

This is a technique, not a one-off. Estimators in this subject are typically defined by formulas that misbehave on assignments the design never produces: DiM divides by #z, the Neyman variance estimator divides by #z - 1, a blocked estimator divides by a per-block count. Rather than carrying side conditions through every lemma, define the estimator totally, prove it equals a well-behaved companion on the support, and transfer the expectation. You will meet expect_congr_of_support again whenever an estimator has a denominator.

The support of complete randomization is described by one lemma, which is immediate from the definition of the pmf:

The support of complete randomizationNew tabOpens with 619 lines of library preamble — the code above is at the bottom of the editor.

/-- Only assignments treating exactly `n₁` units can occur: this describes the **support**
of the design, and is what lets us replace an estimator by an equal one on the support. -/
theorem card_eq_of_prob_ne_zero {z : Assignment n}
    (hz : (completeRandomization n n₁ h).prob z ≠ 0) : #z = n₁ := by
  by_contra hc
  exact hz (by simp [completeRandomization_prob, hc])

DiM equals HT when the design is balanced

Here is the identity that makes the whole thing work.

DiM_eq_HT_of_cardNew tabOpens with 941 lines of library preamble — the code above is at the bottom of the editor.

/-- **On the support of complete randomisation, difference-in-means *is* Horvitz--Thompson.**

Once `#z = n₁` is known, `#z = n₁` and `#zᶜ = n - n₁`, while every propensity equals `n₁/n`.
The factor `1/n` in front of `HT` then cancels against the constant weights, turning the two
inverse-probability sums into the treated and control sample means. -/
theorem DiM_eq_HT_of_card {n₁ : ℕ} (hle : n₁ ≤ n) (P : Population n)
    (hpos : 0 < n₁) (hlt : n₁ < n) {z : Assignment n} (hz : #z = n₁) :
    DiM P z = HT (completeRandomization n n₁ hle) P z := by
  have hnR : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr (by omega)
  have hn1R : (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 hcompl : ((#zᶜ : ℕ) : ℝ) = (n : ℝ) - (n₁ : ℝ) := by
    rw [card_compl, hz, Nat.cast_sub hle]
  -- rewrite each summand of `HT` into a difference of constant-weight terms
  have hterm : ∀ i : Fin n,
      (1 / (n : ℝ)) * (Z z i * P.y1 i / ((n₁ : ℝ) / (n : ℝ))
        - (1 - Z z i) * P.y0 i / (1 - (n₁ : ℝ) / (n : ℝ)))
      = Z z i * P.y1 i / (n₁ : ℝ) - (1 - Z z i) * P.y0 i / ((n : ℝ) - (n₁ : ℝ)) := by
    intro i
    rw [div_div_eq_mul_div, one_sub_div hnR]
    field_simp
  unfold DiM HT
  simp only [completeRandomization_propensity]
  rw [Finset.mul_sum, Finset.sum_congr rfl fun i _ => hterm i, Finset.sum_sub_distrib,
    ← Finset.sum_div, ← Finset.sum_div, sum_mem_eq_sum_Z_mul, sum_compl_eq_sum_one_sub_Z_mul,
    hz, hcompl]

The statistical content is one sentence: when every propensity is the same constant, inverse-probability weighting is just averaging. With $π_{i} = n_{1} / n$ for all $i$ , the treated half of HT is

\frac{1}{n} \sum_{i} \frac{Z_{i} y_{1 i}}{n_{1} / n} = \frac{1}{n_{1}} \sum_{i \in z} y_{1 i},

the treated sample mean, and the control half is the control sample mean. The $1 / n$ out front cancels against the $n$ in the denominator of the weight; that cancellation is the whole proof.

In Lean it takes twenty lines, and almost all of them are casts. The interesting ones:

hcompl : ((#zᶜ : ℕ) : ℝ) = (n : ℝ) - (n₁ : ℝ) is proved by rw [card_compl, hz, Nat.cast_sub hle]. card_compl lives in ℕ, where subtraction is truncated (3 - 5 = 0), so moving it to ℝ requires Nat.cast_sub, which requires the hypothesis n₁ ≤ n. This is the single most common source of wrong autoformalized statements about finite populations, and it is why the hypothesis hle is threaded everywhere.
hterm does the algebra on one summand: div_div_eq_mul_div and one_sub_div put the two weights over a common denominator, and field_simp (clear denominators, given that they are nonzero) finishes. The haves hnR, hn1R, hn0R at the top of the proof exist only to supply one_sub_div and field_simp with the three non-vanishing facts they need.
sum_mem_eq_sum_Z_mul and sum_compl_eq_sum_one_sub_Z_mul convert between “sum over the treated set” and “sum over everybody, weighted by the indicator”. These are the two directions the algebra needs, and having both is why the treated and control sides are symmetric.

With that, DiM_unbiased_completeRandomization has nothing left to do: rewrite the expectation of DiM into the expectation of HT using expect_congr_of_support and card_eq_of_prob_ne_zero, then apply HT_unbiased, discharging its two hypotheses with completeRandomization_propensity (which says $π_{i} = n_{1} / n$ ), div_pos and div_lt_one.

How wrong can the weights be?

HT_unbiased needs the true propensity scores. Suppose an analyst uses some other number $q$ — a misremembered design, a protocol deviation, a “probability of treatment” copied from a different arm of the study. Nothing in the definition of the estimator prevents it:

Inverse-probability weighting with the wrong probabilityNew tabOpens with 1,114 lines of library preamble — the code above is at the bottom of the editor.

/-- The Horvitz--Thompson estimator with the *wrong* weight: the analyst divides by `q`
instead of by the true propensity score.  Nothing in the definition knows what the design is,
so this is a statistic like any other. -/
noncomputable def HTq (P : Population n) (q : ℝ) (z : Assignment n) : ℝ :=
  (1 / (n : ℝ)) * ∑ i, (Z z i * P.y1 i / q - (1 - Z z i) * P.y0 i / (1 - q))

/-- One unit's contribution, with no hypothesis at all: expectation is linear, so the
indicator is replaced by its propensity and nothing else happens.  Compare `expect_HT_term`,
which additionally cancels — and needs `0 < π i < 1` to do so. -/
theorem expect_HTq_term (D : Design n) (P : Population n) (q : ℝ) (i : Fin n) :
    D.expect (fun z => Z z i * P.y1 i / q - (1 - Z z i) * P.y0 i / (1 - q))
      = D.propensity i * (P.y1 i / q) - (1 - D.propensity i) * (P.y0 i / (1 - q)) := by
  rw [D.expect_congr (g := fun z => Z z i * (P.y1 i / q) - (1 - Z z i) * (P.y0 i / (1 - q)))
      fun z => by rw [mul_div_assoc, mul_div_assoc],
    D.expect_sub, D.expect_mul_const, D.expect_mul_const, Design.expect_Z,
    Design.expect_one_sub_Z]

/-- **The bias of the mis-weighted estimator, in closed form.**  Under any design with constant
propensity `p`,

`E[HTq] = (p / q) · mean y1 - ((1 - p) / (1 - q)) · mean y0`,

so the two halves are inflated by `p / q` and `(1 - p) / (1 - q)` respectively.  At `q = p`
both factors are `1` and this is `τ`; at any other `q` it is not. -/
theorem expect_HTq (D : Design n) (P : Population n) {p : ℝ} (q : ℝ)
    (hp : ∀ i, D.propensity i = p) :
    D.expect (HTq P q) = (p / q) * mean P.y1 - ((1 - p) / (1 - q)) * mean P.y0 := by
  have hexpand : D.expect (HTq P q)
      = (1 / (n : ℝ)) * D.expect (fun z => ∑ i, (Z z i * P.y1 i / q
          - (1 - Z z i) * P.y0 i / (1 - q))) :=
    D.expect_const_mul _ (1 / (n : ℝ))
  rw [hexpand, D.expect_sum Finset.univ
    (fun i z => Z z i * P.y1 i / q - (1 - Z z i) * P.y0 i / (1 - q)),
    Finset.sum_congr rfl fun i _ => expect_HTq_term D P q i]
  simp only [hp]
  rw [Finset.sum_congr rfl fun i (_ : i ∈ Finset.univ) =>
      show p * (P.y1 i / q) - (1 - p) * (P.y0 i / (1 - q))
        = (p / q) * P.y1 i - ((1 - p) / (1 - q)) * P.y0 i from by ring,
    Finset.sum_sub_distrib, ← Finset.mul_sum, ← Finset.mul_sum]
  simp only [mean]
  ring

expect_HTq_term is worth comparing with expect_HT_term. It has no hypotheses, because linearity of expectation never needed any: the indicator is replaced by its propensity and the algebra stops there. The cancellation is what required $0 < π_{i} < 1$ , and with the wrong denominator there is nothing to cancel. The bias is exactly the pair of factors $π / q$ and $(1 - π) / (1 - q)$ .

A two-unit instance, with fair coins and $q = 1 / 4$ :

A concrete biased estimatorNew tabOpens with 1,156 lines of library preamble — the code above is at the bottom of the editor.

/-- Fair coins on two units. -/
noncomputable def halfBernoulli : Design 2 :=
  bernoulliDesign 2 (1 / 2) (by norm_num) (by norm_num)

@[simp] theorem halfBernoulli_propensity (i : Fin 2) : halfBernoulli.propensity i = 1 / 2 :=
  bernoulli_propensity i

/-- **A concrete two-unit counterexample.**  Coin flips (`p = 1/2`), an analyst who divides by
`q = 1/4`, and the population `twoUnits`.  The expectation of the estimator is `2`; the
estimand is `1`.  The bias is not asymptotic, not small, and not a modelling assumption: it
is the ratio `p / q = 2`. -/
theorem expect_HTq_halfBernoulli : halfBernoulli.expect (HTq twoUnits (1 / 4)) = 2 := by
  rw [expect_HTq halfBernoulli twoUnits (1 / 4) halfBernoulli_propensity,
    twoUnits_mean_y1, twoUnits_mean_y0]
  norm_num

theorem HTq_biased : halfBernoulli.expect (HTq twoUnits (1 / 4)) ≠ twoUnits.tau := by
  rw [expect_HTq_halfBernoulli, twoUnits_tau]
  norm_num

The expectation is $2$ ; the estimand is $1$ . This is not an asymptotic statement or a simulation: it is an equation between two real numbers, checked by the compiler.

Why difference in means under a Bernoulli design is not here

The library proves HT_unbiased for the Bernoulli design and DiM_unbiased_completeRandomization for complete randomization, and it does not prove that DiM is unbiased under a Bernoulli design. That is not an oversight and not laziness. The statement is false.

Under a Bernoulli design with probability $p$ , the number of treated units $n_{1} = | z |$ is itself random, and it can be $0$ or $n$ . Both degenerate assignments have strictly positive probability:

Pr (z = \emptyset) = (1 - p)^{n} > 0, Pr (z = {1, \dots, n}) = p^{n} > 0 .

On the first, DiM returns -mean y0; on the second, mean y1. Neither has any reason to equal $τ$ , and neither can be argued away, because the “equal on the support” move requires the bad assignments to be off the support. They are not. The positive probability of the empty assignment is an exercise below.

Two units, $p = 1 / 2$ , and twoUnits again make the failure concrete, and at that size the whole expectation is a sum of four terms, so it can simply be computed. All four assignments have probability $1 / 4$ ; DiM is $0$ on $\emptyset$ and $1$ on the other three:

Difference in means is biased under a Bernoulli designNew tabOpens with 1,178 lines of library preamble — the code above is at the bottom of the editor.

/-- Two units have four assignments, so a sum over `Assignment 2` is a sum of four terms.
`decide` checks the enumeration of `univ` by evaluation, and `Finset.sum_insert` peels off one
term at a time (its side condition — the head is not in the tail — is also decidable). -/
theorem sum_assignments_two (f : Assignment 2 → ℝ) :
    ∑ z, f z = f ∅ + f {0} + f {1} + f {0, 1} := by
  rw [show (univ : Finset (Assignment 2)) = {∅, {0}, {1}, {0, 1}} from by decide]
  repeat rw [Finset.sum_insert (by decide)]
  rw [Finset.sum_singleton]
  ring

/-- Difference in means on `twoUnits`: the treated mean is `#z / #z` and the control mean is
`0`, so the estimator is `1` on every assignment except the empty one, where Lean's `x / 0 = 0`
convention makes it `0`.  Note that this is a statement about *all* four assignments, including
the two degenerate ones — nothing here is conditioned away. -/
theorem DiM_twoUnits (z : Assignment 2) :
    DiM twoUnits z = if z = ∅ then 0 else 1 := by
  have hz : DiM twoUnits z = (#z : ℝ) / (#z : ℝ) := by
    rw [DiM]
    simp only [twoUnits_y1, twoUnits_y0, Finset.sum_const, nsmul_eq_mul, mul_one,
      mul_zero, zero_div, sub_zero]
  rw [hz]
  split
  · next h => simp only [h, Finset.card_empty, Nat.cast_zero, div_zero]
  · next h =>
      exact div_self (Nat.cast_ne_zero.mpr
        (Finset.card_ne_zero.mpr (Finset.nonempty_iff_ne_empty.mpr h)))

/-- **The bias of difference in means under a Bernoulli design, computed.**  Fair coins on two
units put mass `1/4` on each of the four assignments; `DiM` is `0` on one of them and `1` on
the other three. -/
theorem expect_DiM_halfBernoulli : halfBernoulli.expect (DiM twoUnits) = 3 / 4 := by
  rw [Design.expect, sum_assignments_two (fun z => halfBernoulli.prob z * DiM twoUnits z)]
  simp only [DiM_twoUnits, halfBernoulli, bernoulliDesign_prob]
  norm_num

/-- **`DiM` is biased under a Bernoulli design.**  The estimand is `1`; the expectation is
`3/4`.  The whole of the gap comes from the single assignment on which nobody is treated. -/
theorem DiM_biased_halfBernoulli :
    halfBernoulli.expect (DiM twoUnits) ≠ twoUnits.tau := by
  rw [expect_DiM_halfBernoulli, twoUnits_tau]
  norm_num

sum_assignments_two is the whole trick. Assignment 2 is Finset (Fin 2), which has four elements, and decide proves univ = {∅, {0}, {1}, {0,1}} by evaluating both sides — there is no cleverness to it, just a finite check the kernel is willing to run. Finset.sum_insert then peels the four terms off one at a time; its side condition (the head is not already in the tail) is decidable too, which is why by decide discharges it. After that the goal is arithmetic in $Q$ and norm_num closes it.

So $E [{\hat{τ}}_{DiM}] = 3 / 4 \neq 1 = τ$ , with the $3 / 4$ checked by the kernel rather than by me. Both degenerate assignments have a $0 / 0$ in them, but only $\emptyset$ costs anything: on the all-treated assignment the missing control mean would have been $\bar{y} (0) = 0$ anyway, so div_zero guesses right by luck. That is what a bias of exactly $1 / 4$ looks like — one assignment in four, and the estimator off by one on it.

Three honest repairs, in increasing order of usefulness:

State a different estimand. $E [{\hat{τ}}_{DiM}]$ exists and is a perfectly well-defined number under a Bernoulli design; it is just not $τ$ . You could compute it, and it would be an ugly expression involving $Pr (| z | = k)$ for every $k$ .
Condition on the realized count. The Bernoulli pmf is constant on ${z : | z | = n_{1}}$ — every such assignment has probability $p^{n_{1}} (1 - p)^{n - n_{1}}$ — so the Bernoulli design conditioned on $| z | = n_{1}$ is exactly complete randomization with $n_{1}$ treated units. Hence $E [{\hat{τ}}_{DiM} ∣ | z | = n_{1}] = τ$ for every $0 < n_{1} < n$ , and by the law of total expectation $E [{\hat{τ}}_{DiM} ∣ 0 < | z | < n] = τ$ as well — a weighted average of unbiased estimators (see the core exercise below). Conditioning is the standard fix, and it is honest about what it costs: the estimator is unbiased for $τ$ under a different probability measure than the one the experiment ran under.
Use HT instead. HT_unbiased_bernoulli holds unconditionally for $0 < p < 1$ , because Horvitz–Thompson divides by the design’s $π_{i} = p$ , which is not random, rather than by the realized $| z |$ , which is.

What the theorem actually rests on

A Lean proof is a term, and the term is checked by a kernel of a few thousand lines. What the kernel cannot check is whether you smuggled in an axiom. #print axioms walks the finished proof term transitively and lists every axiom reached:

The axiom auditNew tabOpens with 1,038 lines of library preamble — the code above is at the bottom of the editor.

/- `#print axioms` walks the whole proof term, transitively, and reports every axiom it
rests on.  It is the only way to be sure that a theorem is not propped up by a `sorry`
(which is itself an axiom, `sorryAx`) somewhere three lemmas down. -/
#print axioms HT_unbiased
-- 'MrCLean.HT_unbiased' depends on axioms: [propext, Classical.choice, Quot.sound]

#print axioms DiM_unbiased_completeRandomization
-- 'MrCLean.DiM_unbiased_completeRandomization' depends on axioms:
--   [propext, Classical.choice, Quot.sound]

Three names come back, and they are the same three that come back for almost every theorem in Mathlib:

propext — propositional extensionality: if p ↔ q then p = q. Propositions that imply each other are equal, so you may rewrite one for the other anywhere, including inside other types.
Classical.choice — a choice function on every non-empty type. This is the axiom of choice, and via Diaconescu’s theorem it delivers excluded middle, which is why by_cases and Classical.em are available. Every statement about the real numbers uses it, since ℝ itself is built with it.
Quot.sound — quotients respect their relation. ℝ is a quotient of Cauchy sequences, ℤ and ℚ are quotients, Finset is built on Multiset which is a quotient of List. Nothing here is optional.

What matters is the name that is absent. sorry is not a hole in the compiler’s checking; it is an axiom, sorryAx, and it shows up in this list the moment any lemma in the transitive dependency graph uses it. A proof three levels down that says sorry cannot hide from #print axioms, which is exactly why this is the check to run on any formalization you did not write yourself. native_decide also announces itself, as Lean.ofReduceBool — worth knowing, because it trusts the compiler and the hardware in a way the kernel does not.

Exercises

ExerciseOne weight cancels (warm-up)

The heart of Horvitz–Thompson in miniature. Design.expect_congr lets you rewrite under the fun z => … binder; Design.expect_mul_const pulls the constant out; Design.expect_Z turns E[Z i] into π i; and mul_div_cancel₀ _ hi cancels π i * (y1 i / π i). The goal after the first three steps is D.propensity i * (P.y1 i / D.propensity i) = P.y1 i.

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

/-- Exercise (warm-up).  **One inverse-probability weight cancels.**  The treated half of a
single Horvitz--Thompson summand has expectation `y1 i` exactly.

Sketch: `Z z i * (y1 i / π i)` is the same function as `Z z i * y1 i / π i`
(`mul_div_assoc`), so `Design.expect_mul_const` pulls the constant out, `Design.expect_Z`
turns `E[Z i]` into `π i`, and `mul_div_cancel₀ _ hi` cancels.  `Design.expect_congr` is how
you rewrite under the binder. -/
theorem ex_expect_one_weight (D : Design n) (P : Population n) (i : Fin n)
    (hi : D.propensity i ≠ 0) :
    D.expect (fun z => Z z i * P.y1 i / D.propensity i) = P.y1 i := by
  sorry

Show solution

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

/-- Exercise (warm-up).  **One inverse-probability weight cancels.**  The treated half of a
single Horvitz--Thompson summand has expectation `y1 i` exactly.

Sketch: `Z z i * (y1 i / π i)` is the same function as `Z z i * y1 i / π i`
(`mul_div_assoc`), so `Design.expect_mul_const` pulls the constant out, `Design.expect_Z`
turns `E[Z i]` into `π i`, and `mul_div_cancel₀ _ hi` cancels.  `Design.expect_congr` is how
you rewrite under the binder. -/
theorem ex_expect_one_weight (D : Design n) (P : Population n) (i : Fin n)
    (hi : D.propensity i ≠ 0) :
    D.expect (fun z => Z z i * P.y1 i / D.propensity i) = P.y1 i := by
  rw [D.expect_congr (g := fun z => Z z i * (P.y1 i / D.propensity i))
      fun z => by rw [mul_div_assoc],
    D.expect_mul_const, Design.expect_Z, mul_div_cancel₀ _ hi]

ExerciseDifference in means with nobody treated (warm-up)

Unfold both definitions and let simp only do the arithmetic. The lemmas you need are Finset.sum_empty, Finset.card_empty, Nat.cast_zero, div_zero (this is where 0 / 0 = 0 happens), Finset.compl_empty, Finset.card_univ, Fintype.card_fin and zero_sub. No hypothesis is required, which is the point: the statement is about every assignment.

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

/-- Exercise (warm-up).  **Difference in means on the degenerate assignment.**  When nobody is
treated, the treated mean is `0 / 0 = 0` and the control mean is the population mean of `y0`,
so `DiM P ∅ = -mean y0`.  No hypothesis is needed, and that is the point: Lean's division
convention makes the estimator total, so a statement about `DiM` is a statement about *every*
assignment unless you say otherwise.

Sketch: `unfold DiM mean`, then `simp only` with `Finset.sum_empty`, `Finset.card_empty`,
`Nat.cast_zero`, `div_zero`, `Finset.compl_empty`, `Finset.card_univ`, `Fintype.card_fin`,
`zero_sub`. -/
theorem ex_DiM_empty (P : Population n) : DiM P ∅ = -mean P.y0 := by
  sorry

Show solution

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

/-- Exercise (warm-up).  **Difference in means on the degenerate assignment.**  When nobody is
treated, the treated mean is `0 / 0 = 0` and the control mean is the population mean of `y0`,
so `DiM P ∅ = -mean y0`.  No hypothesis is needed, and that is the point: Lean's division
convention makes the estimator total, so a statement about `DiM` is a statement about *every*
assignment unless you say otherwise.

Sketch: `unfold DiM mean`, then `simp only` with `Finset.sum_empty`, `Finset.card_empty`,
`Nat.cast_zero`, `div_zero`, `Finset.compl_empty`, `Finset.card_univ`, `Fintype.card_fin`,
`zero_sub`. -/
theorem ex_DiM_empty (P : Population n) : DiM P ∅ = -mean P.y0 := by
  unfold DiM mean
  simp only [Finset.sum_empty, Finset.card_empty, Nat.cast_zero, div_zero,
    Finset.compl_empty, Finset.card_univ, Fintype.card_fin, zero_sub]

ExerciseWeighted averages of unbiased estimators (core)

Nothing about the design is used beyond linearity, so this is a template you will reuse: it is why blocked and stratified estimators inherit unbiasedness from their per-block pieces in chapter 10, and it is repair (2) of the Bernoulli discussion above. Exchange E and ∑ with Design.expect_sum, pull each w k out with Design.expect_const_mul, replace each E[T k] by θ inside the sum with Finset.sum_congr, then use Finset.sum_mul backwards and hw.

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

/-- Exercise (core).  **A weighted average of unbiased estimators is unbiased**, provided the
weights sum to one.  This is the abstract reason blocked and stratified estimators inherit
unbiasedness from their per-block pieces (chapter 10), and it uses nothing about the design
beyond linearity.

Sketch: `Design.expect_sum` exchanges `E` and `∑`, `Design.expect_const_mul` pulls each weight
out, `Finset.sum_congr` replaces each `E[T k]` by `θ`, and `Finset.sum_mul` (backwards) plus
`hw` finishes. -/
theorem ex_expect_weighted_average {ι : Type*} (D : Design n) (s : Finset ι)
    (w : ι → ℝ) (T : ι → Assignment n → ℝ) (θ : ℝ)
    (hw : ∑ k ∈ s, w k = 1) (hT : ∀ k ∈ s, D.expect (T k) = θ) :
    D.expect (fun z => ∑ k ∈ s, w k * T k z) = θ := by
  sorry

Show solution

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

/-- Exercise (core).  **A weighted average of unbiased estimators is unbiased**, provided the
weights sum to one.  This is the abstract reason blocked and stratified estimators inherit
unbiasedness from their per-block pieces (chapter 10), and it uses nothing about the design
beyond linearity.

Sketch: `Design.expect_sum` exchanges `E` and `∑`, `Design.expect_const_mul` pulls each weight
out, `Finset.sum_congr` replaces each `E[T k]` by `θ`, and `Finset.sum_mul` (backwards) plus
`hw` finishes. -/
theorem ex_expect_weighted_average {ι : Type*} (D : Design n) (s : Finset ι)
    (w : ι → ℝ) (T : ι → Assignment n → ℝ) (θ : ℝ)
    (hw : ∑ k ∈ s, w k = 1) (hT : ∀ k ∈ s, D.expect (T k) = θ) :
    D.expect (fun z => ∑ k ∈ s, w k * T k z) = θ := by
  rw [D.expect_sum s (fun k z => w k * T k z)]
  have h : ∀ k ∈ s, D.expect (fun z => w k * T k z) = w k * θ := by
    intro k hk
    rw [D.expect_const_mul (T k) (w k), hT k hk]
  rw [Finset.sum_congr rfl h, ← Finset.sum_mul, hw, one_mul]

ExerciseThe degenerate assignment really happens (core)

This is the formal content of “Bernoulli difference-in-means is not covered”. bernoulliDesign_prob is a @[simp] lemma unfolding the pmf; after Finset.card_empty, pow_zero, one_mul and Nat.sub_zero you are looking at 0 < (1 - p) ^ n, which is pow_pos applied to 0 < 1 - p, which is linarith from hlt.

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

/-- Exercise (core).  **Under a Bernoulli design the degenerate assignment really happens.**
If `p < 1` then `P(nobody treated) = (1 - p)^n > 0`, so the set of assignments on which
`DiM` is degenerate is not a null set and the "agree on the support" argument that proves
`DiM_unbiased_completeRandomization` has nothing to stand on.

Sketch: `bernoulliDesign_prob` unfolds the pmf; `Finset.card_empty`, `pow_zero`, `one_mul`
and `Nat.sub_zero` reduce it to `(1 - p) ^ n`; then `pow_pos`. -/
theorem ex_bernoulli_prob_empty_pos {p : ℝ} (hp0 : 0 ≤ p) (hp1 : p ≤ 1) (hlt : p < 1) :
    0 < (bernoulliDesign n p hp0 hp1).prob ∅ := by
  sorry

Show solution

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

/-- Exercise (core).  **Under a Bernoulli design the degenerate assignment really happens.**
If `p < 1` then `P(nobody treated) = (1 - p)^n > 0`, so the set of assignments on which
`DiM` is degenerate is not a null set and the "agree on the support" argument that proves
`DiM_unbiased_completeRandomization` has nothing to stand on.

Sketch: `bernoulliDesign_prob` unfolds the pmf; `Finset.card_empty`, `pow_zero`, `one_mul`
and `Nat.sub_zero` reduce it to `(1 - p) ^ n`; then `pow_pos`. -/
theorem ex_bernoulli_prob_empty_pos {p : ℝ} (hp0 : 0 ≤ p) (hp1 : p ≤ 1) (hlt : p < 1) :
    0 < (bernoulliDesign n p hp0 hp1).prob ∅ := by
  rw [bernoulliDesign_prob]
  simp only [Finset.card_empty, pow_zero, one_mul, Nat.sub_zero]
  exact pow_pos (by linarith) n

ExerciseThe treated half is unbiased on its own (stretch)

Reconstruct HT_unbiased’s proof skeleton on a simpler estimator: pull the 1 / n out with Design.expect_const_mul, exchange E and ∑ with Design.expect_sum, apply the warm-up exercise to every term with Finset.sum_congr, then simp only [mean] and ring. Note which hypothesis you do not need: nothing here mentions π i < 1, because the control weight never appears.

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

/-- Exercise (stretch).  **The treated half of Horvitz--Thompson, on its own, is unbiased for
`mean y1`.**  Only the lower positivity hypothesis is needed: the control weight never
appears, so `π i < 1` is irrelevant here.

This is `HT_unbiased`'s proof skeleton with half the algebra removed, so it is the exercise to
do if you want to be able to reconstruct the real theorem from memory.  The three steps are:
pull the `1 / n` out (`Design.expect_const_mul`), exchange `E` and `∑` (`Design.expect_sum`),
and apply the warm-up exercise above to every term (`Finset.sum_congr`).  Then `simp only
[mean]` and `ring`. -/
theorem ex_HT_treated_mean_unbiased (D : Design n) (P : Population n)
    (hlo : ∀ i, D.propensity i ≠ 0) :
    D.expect (fun z => (1 / (n : ℝ)) * ∑ i, Z z i * P.y1 i / D.propensity i) = mean P.y1 := by
  sorry

Show solution

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

/-- Exercise (stretch).  **The treated half of Horvitz--Thompson, on its own, is unbiased for
`mean y1`.**  Only the lower positivity hypothesis is needed: the control weight never
appears, so `π i < 1` is irrelevant here.

This is `HT_unbiased`'s proof skeleton with half the algebra removed, so it is the exercise to
do if you want to be able to reconstruct the real theorem from memory.  The three steps are:
pull the `1 / n` out (`Design.expect_const_mul`), exchange `E` and `∑` (`Design.expect_sum`),
and apply the warm-up exercise above to every term (`Finset.sum_congr`).  Then `simp only
[mean]` and `ring`. -/
theorem ex_HT_treated_mean_unbiased (D : Design n) (P : Population n)
    (hlo : ∀ i, D.propensity i ≠ 0) :
    D.expect (fun z => (1 / (n : ℝ)) * ∑ i, Z z i * P.y1 i / D.propensity i) = mean P.y1 := by
  rw [D.expect_const_mul (fun z => ∑ i, Z z i * P.y1 i / D.propensity i) (1 / (n : ℝ)),
    D.expect_sum Finset.univ (fun i z => Z z i * P.y1 i / D.propensity i),
    Finset.sum_congr rfl fun i _ => ex_expect_one_weight D P i (hlo i)]
  simp only [mean]
  ring

Reading autoformalized Lean

Both of the following compile, both look like this chapter’s theorem, and neither says what it appears to say.

A hypothesis nothing satisfies

-- WRONG: do not trust this statement
theorem HT_unbiased_bad (hn : 0 < n) (D : Design n) (P : Population n)
    (hlo : ∀ i, 0 < D.propensity i) (hhi : ∀ i, 1 < D.propensity i) :
    D.expect (HT D P) = P.tau :=
  absurd (hhi ⟨0, hn⟩) (not_lt.mpr (D.propensity_le_one _))

hhi has drifted from D.propensity i < 1 to 1 < D.propensity i. Since Design.propensity_le_one proves $π_{i} \leq 1$ for every design, no D satisfies hhi once there is at least one unit, and the theorem is vacuously true: the “proof” never mentions HT at all, it just derives a contradiction from the hypotheses. The tell here is the shape of the proof — a real unbiasedness proof has to touch the estimator — but the reliable check is to ask whether the hypotheses are satisfiable, by exhibiting one object that meets them. Here you cannot, and #check @HT_unbiased_bad will not warn you.

An estimator that is not a statistic

-- WRONG: this is not an estimator
noncomputable def HT_bad (D : Design n) (P : Population n) (z : Assignment n) : ℝ :=
  (1 / (n : ℝ)) * ∑ i, (P.y1 i - P.y0 i)

theorem HT_bad_unbiased (D : Design n) (P : Population n) :
    D.expect (HT_bad D P) = P.tau := by
  unfold HT_bad
  rw [Design.expect_const, Population.tau_eq]
  ring

Perfectly true, and worthless: z is bound and never used, so HT_bad is the constant function equal to $τ$ , and Design.expect_const — the lemma that says the expectation of a constant is that constant, the one place sum_one gets spent — closes the goal without the design ever being looked at. An “estimator” that ignores the assignment cannot be computed from data, because the data is what the assignment reveals.

Two habits catch this. First, check that the argument standing for the data actually occurs in the body. Lean accepts an unused binder, but its unusedVariables linter does say so — the definition above compiles with Variable name z is not explicitly referenced — so the tell is in the build output if you are reading it. It is also easy to suppress: renaming the binder _z silences the linter and the definition stays exactly as useless. Second, look for the companion lemma: this library states HT_eq_observed and DiM_eq_observed precisely so that “is a statistic” is a checked claim rather than a reader’s assumption. If an autoformalization of an estimator does not come with such a lemma, write one before you believe the unbiasedness theorem.

Two further checks worth making a habit of, both illustrated by this chapter: does every / in the statement have a nonzero denominator, or is a div_zero silently doing work; and is a subtraction happening in ℕ, where n₁ - 2 = 0 when n₁ = 1, rather than in ℝ? The complete-randomization propensity formulas in chapter 7 are stated with real subtraction on the right-hand side for exactly that reason.

The two estimators

Horvitz–Thompson is unbiased

What 0 < π i is doing there

The proof, lemma by lemma

Difference in means, under complete randomization

Equal on the support

DiM equals HT when the design is balanced

How wrong can the weights be?

Why difference in means under a Bernoulli design is not here

What the theorem actually rests on

Exercises

Reading autoformalized Lean

What `0 < π i` is doing there