Causal inferenceChapter 08

Unbiasedness

Horvitz–Thompson is unbiased under any design with strictly interior propensities; difference in means is unbiased under complete randomization.

Contents
  1. The two estimators
  2. Horvitz–Thompson is unbiased
    1. What 0 < π i is doing there
    2. The proof
  3. Difference in means, under complete randomization
    1. Equal on the support
    2. DiM equals HT when the design is balanced
  4. How wrong can the weights be?
  5. Why difference in means under a Bernoulli design is not here
  6. What the theorem actually rests on
  7. Exercises
  8. Reading autoformalized Lean

This chapter proves two results using the finite-population definitions from the preceding chapters:

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

The second result follows from the first by showing that the two estimators agree on the support of complete randomization.

The two estimators

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

τ^HT=1n∑i=1n(Ziy1iπi−(1−Zi)y0i1−πi).
The Horvitz–Thompson estimatorNew tabOpens with 658 lines of library preamble — the code above is at the bottom of the editor.
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 — y1i appears for every unit i, including units that were never treated. It is, and the reason is that the indicator kills the unobserved half of every summand.

Proposition (Horvitz–Thompson is a statistic). For every assignment,

τ^HT=1n∑i=1n(ZiYiobsπi−(1−Zi)Yiobs1−πi),

where Yiobs=Ziy1i+(1−Zi)y0i is the outcome the experiment reveals for unit i.

Proof.

Compare the two expressions summand by summand and split on whether unit i is treated. If i∈z then Zi=1 and Yiobs=y1i: the treated half of the summand is y1i/πi on both sides, and the control half is 0 on both sides. If i∉z then Zi=0 and Yiobs=y0i: the treated half is 0 on both sides and the control half is y0i/(1−πi) on both. Either way the summands agree, and so do the sums.∎

HT only ever touches the revealed outcomeNew tabOpens with 662 lines of library preamble — the code above is at the bottom of the editor.
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]
  1. unfold HT

    Replace HT D P z on the left by its definition, so that both sides of the goal are now 1 / n * ∑ i, … with visibly different summands. unfold is rw with a definition instead of a lemma.

  2. congr 1

    Both sides are products of the same factor 1 / ↑n with a sum. congr 1 peels one layer of application off an equation between two applications of the same head symbol, leaving one goal per differing argument: here the 1 / ↑n factors match syntactically and the only goal left, labelled e_a, is the equality of the two sums. Using congr rather than field_simp keeps n out of the rest of the proof — nothing below needs n≠0.

    Goal
    case e_a
    ⊢ ∑ i, (Z z i * P.y1 i / D.propensity i - (1 - Z z i) * P.y0 i / (1 - D.propensity i)) =
        ∑ i, (Z z i * Yobs P z i / D.propensity i - (1 - Z z i) * Yobs P z i / (1 - D.propensity i))
  3. refine Finset.sum_congr rfl fun i _ => ?_

    Finset.sum_congr states that two sums are equal when their index sets are equal and their summands agree on that set. The rfl discharges “same index set” (both are univ), the fun i _ => ?_ leaves the summand goal for a general unit i, and the _ swallows the unused proof that i ∈ univ. This is the standard way to get from a statement about a sum to a statement about one term.

    Goal
    i : Fin n
    ⊢ Z z i * P.y1 i / D.propensity i - (1 - Z z i) * P.y0 i / (1 - D.propensity i) =
        Z z i * Yobs P z i / D.propensity i - (1 - Z z i) * Yobs P z i / (1 - D.propensity i)
  4. by_cases hi : i ∈ z <;> simp [Yobs_eq_ite, Z, hi]

    The two cases of the informal proof. by_cases adds hi : i ∈ z in the first branch and hi : i ∉ z in the second, and <;> runs the same simp on both. Yobs_eq_ite rewrites Yobs P z i to if i ∈ z then P.y1 i else P.y0 i, unfolding Z exposes the same if i ∈ z then 1 else 0, and hi decides both conditionals. In the treated branch what is left is 1 * y1 i / π i - 0 * y0 i / (1 - π i) on each side, which simp finishes with one_mul, sub_self, zero_mul, zero_div and sub_zero; the control branch is the mirror image, ending in zero_sub.

Difference in means is the estimator everyone actually runs: the average observed outcome among the treated units minus the average among the controls,

τ^DiM(z)=1n1∑i∈zy1i−1n0∑i∉zy0i,n1=|z|,n0=n−n1.
The difference-in-means estimatorNew tabOpens with 702 lines of library preamble — the code above is at the bottom of the editor.
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[τ^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 682 lines of library preamble — the code above is at the bottom of the editor.
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

Each part of the statement has a specific role.

  • 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[τ^HT], i.e. ∑ z, D.prob z * HT D P z. It is a sum over all 2n 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, y1i=1, y0i=0, hence τ=1. Every counterexample in this chapter runs against it.

The population used by every counterexample belowNew tabOpens with 775 lines of library preamble — the code above is at the bottom of the editor.
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]
  1. theorem twoUnits_y1 : twoUnits.y1 = fun _ => (1 : ℝ) := rfl

    A structure’s field is by definition what the anonymous constructor put there, so the projection lemma is rfl. Marking the two of them @[simp] is what lets later proofs about twoUnits be norm_num calls rather than unfold gymnastics: simp will replace twoUnits.y1 by the constant function without being asked.

  2. rw [twoUnits_y1]; exact mean_const (by norm_num) 1

    The mean of the constant vector. rw turns the goal mean twoUnits.y1 = 1 into mean (fun _ => (1 : ℝ)) = 1, and mean_const : 0 < n → mean (fun _ => c) = c closes it. mean_const needs 0<n because mean divides by n and 1 / 0 = 0 would make the claim false at n=0; here n=2, so by norm_num supplies it. twoUnits_mean_y0 is the same two lines with 0 in place of 1.

  3. rw [Population.tau_eq_mean_sub_mean, twoUnits_mean_y1, twoUnits_mean_y0, sub_zero]

    Four rewrites, left to right. tau_eq_mean_sub_mean is chapter 6’s τ=y¯1−y¯0, and it turns the goal into mean twoUnits.y1 - mean twoUnits.y0 = 1. The next two rewrites replace those means by 1 and 0 using the lemmas just proved, leaving 1 - 0 = 1, and sub_zero finishes by rewriting 1 - 0 to 1 — after which the goal is 1 = 1 and rw closes it with rfl automatically.

    Goal
    ⊢ 1 - 0 = 1

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

theorem expect_neverTreat (f : Assignment n → ℝ) : (neverTreat n).expect f = f ∅ := by
  simp [Design.expect, neverTreat]

@[simp] theorem neverTreat_propensity (i : Fin n) : (neverTreat n).propensity i = 0 := by
  rw [Design.propensity, expect_neverTreat]
  simp [Z]
  1. nonneg z := by by_cases h : z = ∅ <;> simp [h]

    The first Design axiom. The pmf is an if, so split on its condition: with h : z = ∅ the value is 1 and with h : z ≠ ∅ it is 0, and simp [h] reduces the conditional and checks 0 ≤ 1 and 0 ≤ 0. A structure instance is a proof obligation per field, and these two by blocks are ordinary tactic proofs that happen to be written inside a definition.

  2. sum_one := by simp

    The second axiom, ∑ z, (if z = ∅ then 1 else 0) = 1. Exactly one term of the sum is non-zero, and simp knows the lemma that collapses such a sum — Finset.sum_ite_eq', which evaluates ∑ z ∈ s, if z = a then f z else 0 to f a when a ∈ s — leaving 1 = 1 after mem_univ discharges the membership.

  3. simp [Design.expect, neverTreat]

    The proof of expect_neverTreat, that expectation under this design is evaluation at ∅. Unfolding expect gives ∑ z, (if z = ∅ then 1 else 0) * f z; ite_mul, one_mul and zero_mul push the conditional outward to if z = ∅ then f z else 0, and Finset.sum_ite_eq' again collapses the sum to f ∅. A point-mass design is the one case where “expectation” and “the value there” are the same thing, which is why every counterexample below is a one-liner.

  4. rw [Design.propensity, expect_neverTreat]; simp [Z]

    The propensity. propensity i is defined as E[Z i], so unfolding it and applying the previous lemma leaves Z ∅ i = 0; unfolding Z makes that (if i ∈ (∅ : Assignment n) then 1 else 0) = 0, and simp closes it with Finset.notMem_empty. Nothing is treated, so every propensity is zero — and neverTreat therefore satisfies the hypothesis of no theorem in this chapter.

Under this design the treated half of each Horvitz–Thompson summand is 0⋅y1i/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.

Proposition. Under the design that never treats anybody, τ^HT=−y¯0 on the only assignment that occurs. Against the two-unit population its expectation is therefore 0, while τ=1.

Proof.

Every unit has Zi=0 and πi=0. The treated half of the i-th summand is 0⋅y1i/0, which the division convention makes 0; the control half is 1⋅y0i/(1−0)=y0i. So the estimator equals 1n∑i(0−y0i)=−y¯0. The design is a point mass, so its expectation is that same number. For the two-unit population y¯0=0 and τ=1, and 0≠1.∎

Without positivity the theorem is false, not vacuousNew tabOpens with 805 lines of library preamble — the code above is at the bottom of the editor.
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

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
  1. have hz : ∀ i : Fin n, Z (∅ : Assignment n) i = 0 := fun i => by simp [Z]

    The pointwise fact the rest of the proof rewrites with: nobody is in the empty set, so the indicator vanishes at every unit. Stating it as a have over all i is what makes the following simp only able to fire under the sum’s binder.

  2. simp only [hz, neverTreat_propensity, zero_mul, zero_div, sub_zero, div_one, one_mul, zero_sub, Finset.sum_neg_distrib]

    The arithmetic of the informal proof, one lemma per move. hz replaces every Z ∅ i by 0 and neverTreat_propensity every π i by 0. Then zero_mul and zero_div kill the treated half (0 * y1 i / 0 = 0); sub_zero turns 1 - 0 into 1 in the control weight; div_one removes the resulting / 1; one_mul removes the 1 *; zero_sub turns each summand 0 - y0 i into -(y0 i); and Finset.sum_neg_distrib pulls the negation out of the sum. Note that this is a simp only, so nothing beyond these nine lemmas is tried.

    Goal
    ⊢ 1 / ↑n * -∑ x, P.y0 x = -mean P.y0
  3. simp only [mean]; ring

    Unfold mean on the right, and the goal is an identity of the commutative field R with ∑ x, P.y0 x as an opaque atom. ring treats division as multiplication by an inverse, so it does not need n≠0 — which matters, because there is no such hypothesis here.

    Goal
    ⊢ 1 / ↑n * -∑ x, P.y0 x = -((∑ x, P.y0 x) / ↑n)
  4. rw [expect_neverTreat, HT_neverTreat_empty, twoUnits_mean_y0, twoUnits_tau]; norm_num

    The counterexample itself. The three rewrites turn (neverTreat 2).expect (HT (neverTreat 2) twoUnits) ≠ twoUnits.tau into -0 ≠ 1: expect_neverTreat evaluates the expectation at ∅, HT_neverTreat_empty replaces the estimator by -mean y0, and the two twoUnits lemmas put in the numbers. norm_num then decides the numeric disequality. Rewriting inside a ≠ is legitimate because a ≠ b is notation for a = b → False, and rw rewrites in the hypothesis part just as happily.

≠ 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 y0i=0, which it generally is not. Note that hhi is not vacuous either — Design.propensity_le_one says πi≤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

Proof.

By linearity of expectation it is enough to compute the expectation of a single summand, since 1/n and the sum over units are both outside the randomization. Fix a unit i. The numbers y1i, y0i, πi are constants — nothing about them depends on the assignment — so

E[Ziy1iπi−(1−Zi)y0i1−πi]=y1iπiE[Zi]−y0i1−πiE[1−Zi].

Now E[Zi]=πi and hence E[1−Zi]=1−πi, and since 0<πi<1 neither denominator is zero, so both weights cancel and the right-hand side is y1i−y0i. Summing over the n units and dividing by n,

E[τ^HT]=1n∑i=1n(y1i−y0i)=τ.

∎

The Lean proof follows that decomposition exactly: a lemma for the single summand, then the bookkeeping.

Step 1: one unit at a time.

Lemma. Fix a unit i with 0<πi<1. Then

E[Ziy1iπi−(1−Zi)y0i1−πi]=y1i−y0i,

the individual treatment effect of that unit.

expect_HT_termNew tabOpens with 670 lines of library preamble — the code above is at the bottom of the editor.
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, preceded by the two non-vanishing facts the last of them needs.

  1. have hne : D.propensity i ≠ 0 := ne_of_gt hlo

    Convert the hypothesis 0<πi into the form the cancellation lemma wants. ne_of_gt : a < b → b ≠ a — the direction matters, and it is the reason this is hne : D.propensity i ≠ 0 rather than 0 ≠ D.propensity i.

  2. have hne' : 1 - D.propensity i ≠ 0 := sub_ne_zero_of_ne (ne_of_gt hhi)

    The same for the control weight. ne_of_gt hhi gives 1 ≠ D.propensity i from πi<1, and sub_ne_zero_of_ne : a ≠ b → a - b ≠ 0 turns that into 1 - π i ≠ 0. These two haves are the only place hlo and hhi are used in this proof: everything before the final cancellation holds for any propensity at all.

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

    First entry of the rw list: re-associate Z z i * y1 i / π i into Z z i * (y1 i / π i), so that everything except the indicator is a single constant factor. Design.expect_congr says two statistics that agree at every z have the same expectation, and it is how you rewrite under the binder — you cannot rw inside a fun z => … directly, so this lemma is the standard workaround. Its argument is the pointwise proof, discharged by mul_div_assoc : a * b / c = a * (b / c) applied once per half.

    Goal
    ⊢ (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
  4. D.expect_sub

    Linearity: the expectation of a difference is the difference of expectations. 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.

  5. D.expect_mul_const, D.expect_mul_const

    Twice, once per half: E[f z * c] = E[f] * c pulls the constants y1 i / π i and y0 i / (1 - π i) out of the expectations, leaving only E[Z i] and E[1 - Z i] inside.

    Goal
    ⊢ (D.expect fun z => Z z i) * (P.y1 i / D.propensity i) -
          (D.expect fun z => 1 - Z z i) * (P.y0 i / (1 - D.propensity i)) =
        P.y1 i - P.y0 i
  6. Design.expect_Z, Design.expect_one_sub_Z

    The indicator expectation step: E[Zi]=πi and E[1−Zi]=1−πi. The first is rfl — true by definition, because propensity is defined as E[Z i]. That is a design decision that pays off here: the propensity is defined to be exactly the expectation this proof needs.

    Goal
    ⊢ D.propensity i * (P.y1 i / D.propensity i) -
          (1 - D.propensity i) * (P.y0 i / (1 - D.propensity i)) =
        P.y1 i - P.y0 i
  7. mul_div_cancel₀ _ hne, mul_div_cancel₀ _ hne'

    The cancellations πi⋅(y1i/πi)=y1i and (1−πi)⋅(y0i/(1−πi))=y0i. mul_div_cancel₀ takes the proof that the denominator is non-zero as its second argument, which is what hne and hne' are for. After both, the goal reads y1 i - y0 i = y1 i - y0 i and the rw closes it.

Step 2: sum over units. With the per-unit fact in hand, the theorem is bookkeeping.

  1. have hexpand : … := D.expect_const_mul _ (1 / (n : ℝ))

    Pull the leading 1/n outside the expectation, which is the ”1/n is not random” half of the informal first sentence. have introduces a named intermediate fact; here its statement is written out in full in the source so that a reader can see the shape of the expression at that point, but the proof is the one lemma expect_const_mul : E[c * f] = c * E[f].

  2. D.expect_sum Finset.univ (fun i z => …)

    Exchange E and ∑i — the other half of the first sentence. 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. The explicit fun i z => … argument is the family of statistics being summed; Lean cannot infer it from the goal, because there are many ways to see an expression as a sum of functions of z.

    Goal
    ⊢ (1 / ↑n * ∑ i, 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.tau
  3. Finset.sum_congr rfl fun i _ => expect_HT_term D P i (hlo i) (hhi i)

    Apply step 1 inside the sum, once per unit. rfl says the index set is unchanged; the function argument supplies, for each i, the proof that its summand equals y1 i - y0 i, with the two pointwise hypotheses obtained by instantiating hlo and hhi at that unit.

    Goal
    ⊢ 1 / ↑n * ∑ i, (P.y1 i - P.y0 i) = P.tau
  4. P.tau_eq

    Unfold the estimand: tau_eq : P.tau = (∑ i, (P.y1 i - P.y0 i)) / n. Only now does τ appear as a formula rather than a name, and the goal has become pure algebra.

    Goal
    ⊢ 1 / ↑n * ∑ i, (P.y1 i - P.y0 i) = (∑ i, (P.y1 i - P.y0 i)) / ↑n
  5. ring

    Close it. 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≠0: 1n⋅S and S/n are the same term after normalization whatever n is.

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; the first two are two-line simp onlys and the third is rfl.

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

HT_unbiased_bernoulliNew tabOpens with 754 lines of library preamble — the code above is at the bottom of the editor.
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)
  1. HT_unbiased _ P (fun i => …) (fun i => …)

    No by: this is a term-mode proof, the general theorem applied to four arguments. The _ is the design, which Lean reads off the goal (bernoulliDesign n p hp0 hp1), P is the population, and the two funs are the pointwise hypotheses ∀ i, 0 < π i and ∀ i, π i < 1 that HT_unbiased demands — supplied as functions of the unit, because that is what a ∀ i hypothesis is.

  2. fun i => by rw [bernoulli_propensity]; exact hlo

    The first hypothesis, one unit at a time. The goal inside the by is 0 < (bernoulliDesign n p hp0 hp1).propensity i; bernoulli_propensity : π i = p rewrites it to 0 < p, which is exactly the assumption hlo. A by block inside a term is an ordinary tactic proof of whatever that position requires — here, of a proposition about one unit.

  3. fun i => by rw [bernoulli_propensity]; exact hhi

    The second, identically: π i < 1 becomes p < 1, which is hhi. The whole corollary is therefore two rewrites, because under a Bernoulli design the propensity does not depend on the unit — and that is the only fact about the design being used.

Difference in means, under complete randomization

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

Formally.

DiM_unbiased_completeRandomizationNew tabOpens with 736 lines of library preamble — the code above is at the bottom of the editor.
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 n1=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 n1=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-n1 subsets of an n-element set needs there to be some.

Proof.

Under complete randomization every unit has the same propensity πi=n1/n, and 0<n1<n puts that number strictly inside (0,1), so Horvitz–Thompson is unbiased for τ by the previous theorem. It is therefore enough to show that the two estimators agree wherever the design puts mass.

Fix an assignment z with exactly n1 treated units. Substituting πi=n1/n into the treated half of Horvitz–Thompson,

1n∑i=1nZiy1in1/n=1n1∑i∈zy1i,

the treated sample mean, and likewise the control half is 1n0∑i∉zy0i with n0=n−n1, because 1−n1/n=n0/n and |zc|=n0. Their difference is τ^DiM(z).

Complete randomization assigns probability zero to every assignment with |z|≠n1, so the two estimators agree on the support, and expectations do not see off-support values. Hence E[τ^DiM]=E[τ^HT]=τ.∎

The Lean proof is that argument in the same order, and it does not touch a single sum over assignments. It needs three tools: transferring an expectation across a change on the support, the description of that support, and the identity between the two estimators.

Equal on the support

Design.expect_congr_of_supportNew tabOpens with 198 lines of library preamble — the code above is at the bottom of the editor.
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.

  1. refine Finset.sum_congr rfl fun z _ => ?_

    Both expectations are sums over all assignments, so compare them one assignment at a time. rfl is “same index set”, and the remaining goal is D.prob z * f z = D.prob z * g z for a fixed z. Note what this step has bought: the quantifier is gone and z is now a constant the hypothesis h can be applied to.

  2. by_cases hz : D.prob z = 0

    The case split the lemma exists for. The hypothesis h says nothing when the probability is zero, so that case has to be handled by arithmetic instead.

  3. rw [hz, zero_mul, zero_mul]

    The zero branch. Rewriting with hz makes both sides 0 * …, and zero_mul turns each into 0. No hypothesis about f or g is needed, which is the entire point of the lemma: off the support the two statistics may disagree wildly and the expectation cannot tell.

  4. rw [h z hz]

    The other branch, whose state on entry is the one shown below. Now hz : ¬D.prob z = 0 is exactly what h wants, so h z hz : f z = g z rewrites the left-hand side into the right and the goal closes by reflexivity.

    Goal
    z : Assignment n
    hz : ¬D.prob z = 0
    ⊢ D.prob z * f z = D.prob z * g z

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.

Lemma. Under complete randomization with n1 treated units, every assignment of positive probability has exactly n1 treated units.

Proof.

The pmf is 1/(nn1) on assignments of size n1 and 0 on all others, so an assignment of a different size has probability zero.∎

The support of complete randomizationNew tabOpens with 465 lines of library preamble — the code above is at the bottom of the editor.
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])
  1. by_contra hc

    The statement is an equation between two naturals and the hypothesis is a disequality, so the contrapositive is the convenient direction: assume hc : ¬#z = n₁ and derive a contradiction. After by_contra the goal is False and both hz and hc are available.

    Goal
    hc : ¬#z = n₁
    ⊢ False
  2. exact hz (by simp [completeRandomization_prob, hc])

    hz : D.prob z ≠ 0 is, unfolded, a function from D.prob z = 0 to False, so applying it to a proof that the probability is zero closes the goal. That proof is the inner by simp: completeRandomization_prob rewrites D.prob z to if #z = n₁ then (n.choose n₁ : ℝ)⁻¹ else 0, and hc sends the conditional down its else-branch.

DiM equals HT when the design is balanced

Here is the identity that makes the whole thing work.

Proposition. Let 0<n1<n and let z be an assignment with exactly n1 treated units. Then, with πi=n1/n for every unit,

τ^DiM(z)=τ^HT(z).
Proof.

The statistical content is one sentence: when every propensity is the same constant, inverse-probability weighting is just averaging. With πi=n1/n the treated half of Horvitz–Thompson is

1n∑iZiy1in1/n=1n1∑i∈zy1i,

the treated sample mean, since the 1/n out front cancels against the n in the numerator of the weight and ∑iZiy1i=∑i∈zy1i. For the control half, 1−n1/n=(n−n1)/n and the control group has exactly n−n1 units, so the same cancellation turns it into ∑i∉zy0i divided by n−n1. Subtracting gives difference in means. That cancellation is the whole proof.∎

DiM_eq_HT_of_cardNew tabOpens with 713 lines of library preamble — the code above is at the bottom of the editor.
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]
  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]

In Lean it takes twenty lines, and almost all of them are casts.

  1. have hnR : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr (by omega)

    The first of three non-vanishing denominators. Nat.cast_ne_zero : (↑m ≠ 0) ↔ (m ≠ 0) moves the claim from R to N, where omega — a decision procedure for linear arithmetic over ℕ and ℤ — gets n ≠ 0 out of hpos : 0 < n₁ and hle : n₁ ≤ n. hn1R : (n₁ : ℝ) ≠ 0 is the same line.

  2. have hn0R : (n : ℝ) - (n₁ : ℝ) ≠ 0 := by …

    The third, and the one that has to be done in R: cast hlt to (n₁ : ℝ) < (n : ℝ) with exact_mod_cast and let linarith conclude (n : ℝ) - (n₁ : ℝ) ≠ 0. Doing it as ((n - n₁ : ℕ) : ℝ) ≠ 0 instead would be a different, weaker statement, because natural subtraction truncates.

  3. have hcompl : ((#zᶜ : ℕ) : ℝ) = (n : ℝ) - (n₁ : ℝ) := by rw [card_compl, hz, Nat.cast_sub hle]

    The size of the control group, as a real number. card_compl is chapter 6’s #zᶜ = n - #z, which lives in ℕ where subtraction is truncated (3 - 5 = 0); hz replaces #z by n₁; and Nat.cast_sub hle is what licenses moving the now-ℕ subtraction across the cast, which it does only because 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.

  4. intro i

    Now the algebra on one summand, stated as hterm for all i so that it can later be applied under the sum’s binder. intro i fixes the unit; the goal is the i-th summand of Horvitz–Thompson with πi already substituted, against the i-th summand of the difference-in-means form.

    Goal
    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₁)
  5. rw [div_div_eq_mul_div, one_sub_div hnR]

    Put the two weights over a common denominator. div_div_eq_mul_div : a / (b / c) = a * c / b turns dividing by n1/n into multiplying by n and dividing by n1, and one_sub_div, whose hypothesis is hnR, rewrites 1 - n₁ / n to (n - n₁) / n. Both moves are exactly the two sentences of the informal proof.

    Goal
    i : Fin n
    ⊢ 1 / ↑n * (Z z i * P.y1 i * ↑n / ↑n₁ - (1 - Z z i) * P.y0 i / ((↑n - ↑n₁) / ↑n)) =
        Z z i * P.y1 i / ↑n₁ - (1 - Z z i) * P.y0 i / (↑n - ↑n₁)
  6. field_simp

    Clear the denominators and finish the summand. field_simp needs to know that each denominator is non-zero, and it finds hnR, hn1R and hn0R in the context — which is the only reason those three haves are at the top of the proof.

  7. unfold DiM HT

    Back to the main goal. Both estimators are replaced by their definitions, so the goal is an equation between (∑ i ∈ z, y1 i)/#z - (∑ i ∈ zᶜ, y0 i)/#zᶜ and 1 / n * ∑ i, (…).

  8. simp only [completeRandomization_propensity]

    Substitute πi=n1/n. It has to be a simp only rather than a rw, because the propensity occurs under the sum’s binder, where rw cannot reach; simp only rewrites everywhere, including inside funs and sums.

    Goal
    ⊢ (∑ i ∈ z, P.y1 i) / ↑(#z) - (∑ i ∈ zᶜ, P.y0 i) / ↑(#zᶜ) =
        1 / ↑n * ∑ x, (Z z x * P.y1 x / (↑n₁ / ↑n) - (1 - Z z x) * P.y0 x / (1 - ↑n₁ / ↑n))
  9. Finset.mul_sum

    First entry of the closing rw list: c * ∑ i, f i = ∑ i, c * f i moves the 1 / n inside the sum, so that hterm — which is a statement about one summand, with the 1 / n included — becomes applicable.

  10. Finset.sum_congr rfl fun i _ => hterm i

    Apply the summand identity to every unit. The right-hand side is now the sum of the difference-in-means summands, with the constants n1 and n−n1 as denominators and no π anywhere.

    Goal
    ⊢ (∑ i ∈ z, P.y1 i) / ↑(#z) - (∑ i ∈ zᶜ, P.y0 i) / ↑(#zᶜ) =
        ∑ i, (Z z i * P.y1 i / ↑n₁ - (1 - Z z i) * P.y0 i / (↑n - ↑n₁))
  11. Finset.sum_sub_distrib, ← Finset.sum_div, ← Finset.sum_div

    Split the sum of differences into a difference of sums, then pull each constant denominator back out of its sum — the ← runs Finset.sum_div : (∑ i, f i) / c = ∑ i, f i / c backwards, from the sum of quotients to the quotient of the sum. The right-hand side now has the shape of the left: one indicator-weighted sum over n1, minus another over n−n1.

    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₁)
  12. sum_mem_eq_sum_Z_mul, sum_compl_eq_sum_one_sub_Z_mul, hz, hcompl

    The last four rewrites convert the left side. sum_mem_eq_sum_Z_mul is ∑ i ∈ z, v i = ∑ i, Z z i * v i and sum_compl_eq_sum_one_sub_Z_mul is its control-side counterpart ∑ i ∈ zᶜ, v i = ∑ i, (1 - Z z i) * v i: they trade “sum over a group” for “sum over everybody, weighted by the indicator”, which is what makes the two sides comparable at all. Then hz and hcompl replace the two group sizes by n₁ and (n : ℝ) - (n₁ : ℝ), and the two sides are syntactically equal.

With that, DiM_unbiased_completeRandomization has nothing left to do.

  1. have hnR : (0 : ℝ) < (n : ℝ) := by exact_mod_cast Nat.lt_of_lt_of_le hpos hle

    The three numeric facts the last two lines will need, all obtained by casting the natural-number hypotheses into R: 0<n (from 0<n1≤n), then hn1R : 0 < (n₁ : ℝ) and hltR : (n₁ : ℝ) < (n : ℝ) by exact_mod_cast from hpos and hlt. exact_mod_cast closes a goal that differs from a hypothesis only by casts.

  2. have hsupport : ∀ z, … → DiM P z = HT … P z := fun z hz => DiM_eq_HT_of_card hle P hpos hlt (card_eq_of_prob_ne_zero hz)

    The two previous results, composed. Given an assignment of positive probability, card_eq_of_prob_ne_zero extracts #z = n₁ and feeds it to DiM_eq_HT_of_card, which returns the identity at that assignment. This one have is the whole “agree on the support” step.

  3. rw [(completeRandomization n n₁ hle).expect_congr_of_support hsupport]

    Transfer the expectation. The goal changes from expect (DiM P) = P.tau to expect (HT … P) = P.tau, and no sum over assignments has been touched.

  4. refine HT_unbiased _ P (fun i => ?_) (fun i => ?_)

    Quote the theorem. refine applies it but leaves the two pointwise hypotheses as goals — that is what the ?_ holes mean — one for 0 < π i and one for π i < 1, each for a general unit i.

  5. rw [completeRandomization_propensity]; exact div_pos hn1R hnR

    The first hole. Rewriting the propensity to n₁ / n reduces it to 0 < (n₁ : ℝ) / (n : ℝ), which is div_pos applied to the two positivity facts from the first step.

  6. rw [completeRandomization_propensity]; exact (div_lt_one hnR).mpr hltR

    The second hole. div_lt_one hnR : a / b < 1 ↔ a < b (for 0 < b), so the goal n₁ / n < 1 follows from hltR : (n₁ : ℝ) < (n : ℝ). This is the step that would fail at n1=n, and it is the only place hlt is used in this theorem.

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

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]

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

Proposition. Let the design have constant propensity πi=p for every unit, and let the analyst weight by q instead. Then

E[1n∑i(Ziy1iq−(1−Zi)y0i1−q)]=pqy¯1−1−p1−qy¯0.
Proof.

By linearity, and because q is a constant, the expectation of the i-th summand is p(y1i/q)−(1−p)(y0i/(1−q)); no hypothesis on p or q is needed for this, since nothing is being cancelled. Rewriting each summand as (p/q)y1i−((1−p)/(1−q))y0i, pulling the two constant factors out of the sum and dividing by n gives the claim. At q=p the two factors are 1 and the estimand τ=y¯1−y¯0 comes back; otherwise the estimator is biased by exactly those two factors.∎

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

    The proof of expect_HTq_term is expect_HT_term’s proof with the two haves and the last two rewrites deleted, so it opens with the same move: re-associate under the binder so that the analyst’s weights become constant factors.

  2. D.expect_sub, D.expect_mul_const, D.expect_mul_const, Design.expect_Z, Design.expect_one_sub_Z

    Linearity, then the two constants out, then the two indicator expectations — and stop. The goal is now closed, because the statement of expect_HTq_term keeps π i * (y1 i / q) unevaluated instead of claiming a cancellation. That is why the lemma has no hypotheses at all: linearity of expectation never needed any, and the cancellation is what required 0<πi<1.

  3. rw [hexpand, D.expect_sum Finset.univ (fun i z => …), Finset.sum_congr rfl fun i _ => expect_HTq_term D P q i]

    expect_HTq now repeats step 2 of HT_unbiased verbatim: hexpand pulls the 1 / n out with expect_const_mul, expect_sum exchanges E and ∑i, and sum_congr applies the per-unit lemma inside the sum.

    Goal
    ⊢ 1 / ↑n * ∑ i, (D.propensity i * (P.y1 i / q) - (1 - D.propensity i) * (P.y0 i / (1 - q))) =
        p / q * mean P.y1 - (1 - p) / (1 - q) * mean P.y0
  4. simp only [hp]

    Use the constancy of the design’s propensity: hp : ∀ i, D.propensity i = p replaces D.propensity i by p under the binder, which rw could not do. This is the only assumption about the design in the whole lemma.

    Goal
    ⊢ 1 / ↑n * ∑ x, (p * (P.y1 x / q) - (1 - p) * (P.y0 x / (1 - q))) =
        p / q * mean P.y1 - (1 - p) / (1 - q) * mean P.y0
  5. rw [Finset.sum_congr rfl fun i _ => 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]

    Regroup each summand so that the bias factors stand in front of the outcomes. show … from by ring is an inline anonymous lemma: it states the summand identity and proves it by ring. Writing it inline avoids a have whose statement would have to be repeated.

    Goal
    ⊢ 1 / ↑n * ∑ i, (p / q * P.y1 i - (1 - p) / (1 - q) * P.y0 i) =
        p / q * mean P.y1 - (1 - p) / (1 - q) * mean P.y0
  6. Finset.sum_sub_distrib, ← Finset.mul_sum, ← Finset.mul_sum

    Split the sum and pull both constant factors out — ← Finset.mul_sum runs c * ∑ i, f i = ∑ i, c * f i from right to left, which is the direction that extracts a constant. Two sums over the population are left.

    Goal
    ⊢ 1 / ↑n * (p / q * ∑ i, P.y1 i - (1 - p) / (1 - q) * ∑ i, P.y0 i) =
        p / q * mean P.y1 - (1 - p) / (1 - q) * mean P.y0
  7. simp only [mean]; ring

    Unfold mean on the right-hand side — the state below — so that both sides are built from ∑ i, P.y1 i, ∑ i, P.y0 i and the constants; ring then matches them. As in HT_unbiased, no hypothesis n≠0 is needed, because ring never divides.

    Goal
    ⊢ 1 / ↑n * (p / q * ∑ i, P.y1 i - (1 - p) / (1 - q) * ∑ i, P.y0 i) =
        p / q * ((∑ i, P.y1 i) / ↑n) - (1 - p) / (1 - q) * ((∑ i, P.y0 i) / ↑n)

A two-unit instance, with fair coins and q=1/4: the true propensity is 1/2, so the bias factors are 1/21/4=2 and 1/23/4=2/3, and against twoUnits (where y¯1=1 and y¯0=0) the expectation is 2⋅1−23⋅0=2.

A concrete biased estimatorNew tabOpens with 850 lines of library preamble — the code above is at the bottom of the editor.
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

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
  1. noncomputable def halfBernoulli : Design 2 := bernoulliDesign 2 (1 / 2) (by norm_num) (by norm_num)

    The design: independent fair coins on two units. The two by norm_num arguments are the proofs of 0 ≤ 1/2 and 1/2 ≤ 1 that bernoulliDesign requires — a definition can take proofs as arguments, and here they are inlined rather than named because nothing else will ever refer to them.

  2. theorem halfBernoulli_propensity (i : Fin 2) : halfBernoulli.propensity i = 1 / 2 := bernoulli_propensity i

    A term-mode restatement of bernoulli_propensity at this design, marked @[simp] so later norm_num calls can use it. It type-checks because halfBernoulli is definitionally that Bernoulli design, not merely equal to it.

  3. rw [expect_HTq halfBernoulli twoUnits (1 / 4) halfBernoulli_propensity, twoUnits_mean_y1, twoUnits_mean_y0]; norm_num

    The computation. The general formula is applied with p=1/2 (supplied implicitly by the hp argument halfBernoulli_propensity) and q=1/4, then the two population means are replaced by 1 and 0, leaving (1/2) / (1/4) * 1 - (1 - 1/2) / (1 - 1/4) * 0 = 2 for norm_num to check.

  4. rw [expect_HTq_halfBernoulli, twoUnits_tau]; norm_num

    And the bias, recorded as a disequality: the two rewrites reduce the goal to (2 : ℝ) ≠ 1, which norm_num decides. 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 n1=|z| is itself random, and it can be 0 or n. Both degenerate assignments have strictly positive probability:

Pr(z=∅)=(1−p)n>0,Pr(z={1,…,n})=pn>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.

Proposition. For the two-unit population with y1i=1 and y0i=0 under independent fair coins, E[τ^DiM]=3/4≠1=τ.

Proof.

There are four assignments and each has probability 1/4. On ∅ both group averages degenerate: the treated mean is 0/0=0 and the control mean is y¯0=0, so the estimator is 0. On each of the other three the treated mean is 1 and the control mean is 0 — the control mean of an empty control group is again 0/0=0, and here that is also the right answer — so the estimator is 1. Hence the expectation is 14(0+1+1+1)=34, while τ=1.∎

Difference in means is biased under a Bernoulli designNew tabOpens with 865 lines of library preamble — the code above is at the bottom of the editor.
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

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

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

theorem DiM_biased_halfBernoulli :
    halfBernoulli.expect (DiM twoUnits) ≠ twoUnits.tau := by
  rw [expect_DiM_halfBernoulli, twoUnits_tau]
  norm_num
  1. rw [show (univ : Finset (Assignment 2)) = {∅, {0}, {1}, {0, 1}} from by decide]

    sum_assignments_two is the whole trick. Assignment 2 is Finset (Fin 2), which has four elements, and decide proves the enumeration by evaluating both sides — there is no cleverness to it, just a finite check the kernel is willing to run. Wrapping it in show … from supplies the statement, since decide needs a closed proposition to evaluate.

  2. repeat rw [Finset.sum_insert (by decide)]

    Peel the four terms off one at a time. Finset.sum_insert requires that the head is not already in the tail, and that side condition is decidable too, which is why by decide discharges it; repeat applies the rewrite until it no longer matches. Finset.sum_singleton then handles the last element and ring reassociates the additions.

    Goal
    ⊢ f ∅ + (f {0} + (f {1} + ∑ x ∈ {{0, 1}}, f x)) = f ∅ + f {0} + f {1} + f {0, 1}
  3. rw [DiM]; simp only [twoUnits_y1, twoUnits_y0, Finset.sum_const, nsmul_eq_mul, mul_one, mul_zero, zero_div, sub_zero]

    DiM_twoUnits next, and this have is where the population’s flatness is spent. Both outcome vectors are constant, so Finset.sum_const turns ∑ i ∈ z, (1 : ℝ) into #z • 1 and nsmul_eq_mul into (#z : ℝ) * 1; mul_one cleans that up, mul_zero and zero_div annihilate the control half, and sub_zero deletes it. What is left is the claim that the estimator is (#z : ℝ) / (#z : ℝ) — the treated sample mean of a vector of ones.

  4. rw [hz]; split

    The goal is now an if, so split does a case analysis on its condition, naming the two branch hypotheses with next h => …. The point of the whole lemma is that x / x is not unconditionally 1.

    Goal
    ⊢ ↑(#z) / ↑(#z) = if z = ∅ then 0 else 1
  5. next h => simp only [h, Finset.card_empty, Nat.cast_zero, div_zero]

    The empty branch, where h : z = ∅. Rewriting with it makes the cardinality #(∅), which Finset.card_empty sends to 0 and Nat.cast_zero to (0 : ℝ), and then div_zero gives 0 / 0 = 0. This is the assignment that costs the estimator its unbiasedness.

  6. next h => exact div_self (Nat.cast_ne_zero.mpr (Finset.card_ne_zero.mpr (Finset.nonempty_iff_ne_empty.mpr h)))

    The other branch, where h : ¬z = ∅, and the goal is (#z : ℝ) / (#z : ℝ) = 1. div_self needs the denominator non-zero, and the argument is that chain read inside out: non-empty from h by Finset.nonempty_iff_ne_empty, hence #z ≠ 0 by Finset.card_ne_zero, hence (#z : ℝ) ≠ 0 by Nat.cast_ne_zero. Three lemmas to move one fact across two type boundaries is typical.

  7. rw [Design.expect, sum_assignments_two (fun z => halfBernoulli.prob z * DiM twoUnits z)]

    expect_DiM_halfBernoulli puts the pieces together. Unfolding expect exposes the sum over assignments, and the enumeration lemma — applied to the specific summand prob z * DiM twoUnits z — turns it into four named terms.

    Goal
    ⊢ halfBernoulli.prob ∅ * DiM twoUnits ∅ + halfBernoulli.prob {0} * DiM twoUnits {0} +
            halfBernoulli.prob {1} * DiM twoUnits {1} +
          halfBernoulli.prob {0, 1} * DiM twoUnits {0, 1} =
        3 / 4
  8. simp only [DiM_twoUnits, halfBernoulli, bernoulliDesign_prob]; norm_num

    Evaluate. DiM_twoUnits replaces each estimator value by its if, bernoulliDesign_prob replaces each probability by p ^ #z * (1 - p) ^ (n - #z) with p = 1/2, and norm_num computes the four cardinalities, the four conditionals and the arithmetic, arriving at 3 / 4. The final DiM_biased_halfBernoulli is then two rewrites and a norm_num on (3 : ℝ) / 4 ≠ 1.

So E[τ^DiM]=3/4≠1=τ, with the 3/4 checked by the kernel rather than by me. Both degenerate assignments have a 0/0 in them, but only ∅ costs anything: on the all-treated assignment the missing control mean would have been 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:

  1. State a different estimand. E[τ^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.
  2. Condition on the realized count. The Bernoulli pmf is constant on {z:|z|=n1} — every such assignment has probability pn1(1−p)n−n1 — so the Bernoulli design conditioned on |z|=n1 is exactly complete randomization with n1 treated units. Hence E[τ^DiM∣|z|=n1]=τ for every 0<n1<n, and by the law of total expectation E[τ^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.
  3. 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 770 lines of library preamble — the code above is at the bottom of the editor.
#print axioms HT_unbiased

#print axioms DiM_unbiased_completeRandomization

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: the treated half of a single summand has expectation exactly y1 i, given only that π i ≠ 0.

New tabOpens with 896 lines of library preamble — the code above is at the bottom of the editor.
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
Hint

Z z i * y1 i / π i is the same function as Z z i * (y1 i / π i) (mul_div_assoc), and Design.expect_congr is how you rewrite under the fun z => … binder; then Design.expect_mul_const pulls the constant out and Design.expect_Z turns E[Z i] into π i. The goal at that point is D.propensity i * (P.y1 i / D.propensity i) = P.y1 i, which is mul_div_cancel₀ _ hi.

Show solution
New tabOpens with 896 lines of library preamble — the code above is at the bottom of the editor.
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)

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 required, which is the point: the statement is about every assignment.

New tabOpens with 904 lines of library preamble — the code above is at the bottom of the editor.
theorem ex_DiM_empty (P : Population n) : DiM P ∅ = -mean P.y0 := by
  sorry
Hint

unfold DiM mean, then let one simp only do the arithmetic with 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.

Show solution
New tabOpens with 904 lines of library preamble — the code above is at the bottom of the editor.
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)

A weighted average of unbiased estimators is unbiased, provided the weights sum to one. 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.

New tabOpens with 909 lines of library preamble — the code above is at the bottom of the editor.
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
Hint

Exchange E and ∑ with D.expect_sum s (fun k z => w k * T k z). Then establish ∀ k ∈ s, D.expect (fun z => w k * T k z) = w k * θ — one D.expect_const_mul and hT k hk each — apply it inside the sum with Finset.sum_congr rfl, and ← Finset.sum_mul, hw and one_mul finish.

Show solution
New tabOpens with 909 lines of library preamble — the code above is at the bottom of the editor.
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)

Under a Bernoulli design with p<1, the assignment that treats nobody has strictly positive probability, (1−p)n>0. This is the formal content of “Bernoulli difference-in-means is not covered”.

New tabOpens with 919 lines of library preamble — the code above is at the bottom of the editor.
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
Hint

rw [bernoulliDesign_prob] unfolds the pmf; after simp only with 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.

Show solution
New tabOpens with 919 lines of library preamble — the code above is at the bottom of the editor.
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)

The treated half of Horvitz–Thompson, on its own, is unbiased for mean y1. Note which hypothesis you do not need: nothing here mentions π i < 1, because the control weight never appears.

New tabOpens with 925 lines of library preamble — the code above is at the bottom of the editor.
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
Hint

This is HT_unbiased’s proof skeleton with half the algebra removed: pull the 1 / n out with D.expect_const_mul, exchange E and ∑ with D.expect_sum Finset.univ, and apply the warm-up ex_expect_one_weight D P i (hlo i) to every term with Finset.sum_congr rfl; then simp only [mean] and ring.

Show solution
New tabOpens with 925 lines of library preamble — the code above is at the bottom of the editor.
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.

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.