Causal inferenceChapter 08
Unbiasedness
Horvitz–Thompson is unbiased under any design with strictly interior propensities; difference in means is unbiased under complete randomization.
Contents
This chapter proves two results using the finite-population definitions from the preceding chapters:
- the Horvitz–Thompson estimator is unbiased for
under any design whose propensity scores are strictly inside ; - 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:
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 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
Proposition (Horvitz–Thompson is a statistic). For every assignment,
where
Compare the two expressions summand by summand and split on whether unit
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]unfold HTReplace
HT D P zon the left by its definition, so that both sides of the goal are now1 / n * ∑ i, …with visibly different summands.unfoldisrwwith a definition instead of a lemma.congr 1Both sides are products of the same factor
1 / ↑nwith a sum.congr 1peels one layer of application off an equation between two applications of the same head symbol, leaving one goal per differing argument: here the1 / ↑nfactors match syntactically and the only goal left, labellede_a, is the equality of the two sums. Usingcongrrather thanfield_simpkeepsnout of the rest of the proof — nothing below needs .Goalcase 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))refine Finset.sum_congr rfl fun i _ => ?_Finset.sum_congrstates that two sums are equal when their index sets are equal and their summands agree on that set. Therfldischarges “same index set” (both areuniv), thefun i _ => ?_leaves the summand goal for a general uniti, and the_swallows the unused proof thati ∈ univ. This is the standard way to get from a statement about a sum to a statement about one term.Goali : 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)by_cases hi : i ∈ z <;> simp [Yobs_eq_ite, Z, hi]The two cases of the informal proof.
by_casesaddshi : i ∈ zin the first branch andhi : i ∉ zin the second, and<;>runs the samesimpon both.Yobs_eq_iterewritesYobs P z itoif i ∈ z then P.y1 i else P.y0 i, unfoldingZexposes the sameif i ∈ z then 1 else 0, andhidecides both conditionals. In the treated branch what is left is1 * y1 i / π i - 0 * y0 i / (1 - π i)on each side, whichsimpfinishes withone_mul,sub_self,zero_mul,zero_divandsub_zero; the control branch is the mirror image, ending inzero_sub.
Difference in means is the estimator everyone actually runs: the average observed outcome among the treated units minus the average among the controls,
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.
example (x : ℝ) : x / 0 = 0 := div_zero xHorvitz–Thompson is unbiased
Informally. Fix a finite population and a design. If every unit has
Formally.
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]
ringEach part of the statement has a specific role.
D : Design nis 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 nis universally quantified too. The theorem holds for every configuration of potential outcomes, which is what “design-based” means: the randomness is entirely inD.D.expect (HT D P)is , i.e.∑ z, D.prob z * HT D P z. It is a sum over all assignments.P.tauis the population ATE,mean P.effect. It is a fixed real number, not a random variable.hloandhhiare 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,
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]theorem twoUnits_y1 : twoUnits.y1 = fun _ => (1 : ℝ) := rflA 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 abouttwoUnitsbenorm_numcalls rather thanunfoldgymnastics:simpwill replacetwoUnits.y1by the constant function without being asked.rw [twoUnits_y1]; exact mean_const (by norm_num) 1The mean of the constant vector.
rwturns the goalmean twoUnits.y1 = 1intomean (fun _ => (1 : ℝ)) = 1, andmean_const : 0 < n → mean (fun _ => c) = ccloses it.mean_constneeds becausemeandivides bynand1 / 0 = 0would make the claim false at ; here , soby norm_numsupplies it.twoUnits_mean_y0is the same two lines with0in place of1.rw [Population.tau_eq_mean_sub_mean, twoUnits_mean_y1, twoUnits_mean_y0, sub_zero]Four rewrites, left to right.
tau_eq_mean_sub_meanis chapter 6’s , and it turns the goal intomean twoUnits.y1 - mean twoUnits.y0 = 1. The next two rewrites replace those means by1and0using the lemmas just proved, leaving1 - 0 = 1, andsub_zerofinishes by rewriting1 - 0to1— after which the goal is1 = 1andrwcloses it withrflautomatically.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
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]nonneg z := by by_cases h : z = ∅ <;> simp [h]The first
Designaxiom. The pmf is anif, so split on its condition: withh : z = ∅the value is1and withh : z ≠ ∅it is0, andsimp [h]reduces the conditional and checks0 ≤ 1and0 ≤ 0. A structure instance is a proof obligation per field, and these twobyblocks are ordinary tactic proofs that happen to be written inside a definition.sum_one := by simpThe second axiom,
∑ z, (if z = ∅ then 1 else 0) = 1. Exactly one term of the sum is non-zero, andsimpknows the lemma that collapses such a sum —Finset.sum_ite_eq', which evaluates∑ z ∈ s, if z = a then f z else 0tof awhena ∈ s— leaving1 = 1aftermem_univdischarges the membership.simp [Design.expect, neverTreat]The proof of
expect_neverTreat, that expectation under this design is evaluation at∅. Unfoldingexpectgives∑ z, (if z = ∅ then 1 else 0) * f z;ite_mul,one_mulandzero_mulpush the conditional outward toif z = ∅ then f z else 0, andFinset.sum_ite_eq'again collapses the sum tof ∅. 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.rw [Design.propensity, expect_neverTreat]; simp [Z]The propensity.
propensity iis defined asE[Z i], so unfolding it and applying the previous lemma leavesZ ∅ i = 0; unfoldingZmakes that(if i ∈ (∅ : Assignment n) then 1 else 0) = 0, andsimpcloses it withFinset.notMem_empty. Nothing is treated, so every propensity is zero — andneverTreattherefore satisfies the hypothesis of no theorem in this chapter.
Under this design the treated half of each Horvitz–Thompson summand is
Proposition. Under the design that never treats anybody,
Every unit has
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_numhave 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
haveover alliis what makes the followingsimp onlyable to fire under the sum’s binder.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.
hzreplaces everyZ ∅ iby0andneverTreat_propensityeveryπ iby0. Thenzero_mulandzero_divkill the treated half (0 * y1 i / 0 = 0);sub_zeroturns1 - 0into1in the control weight;div_oneremoves the resulting/ 1;one_mulremoves the1 *;zero_subturns each summand0 - y0 iinto-(y0 i); andFinset.sum_neg_distribpulls the negation out of the sum. Note that this is asimp only, so nothing beyond these nine lemmas is tried.Goal⊢ 1 / ↑n * -∑ x, P.y0 x = -mean P.y0simp only [mean]; ringUnfold
meanon the right, and the goal is an identity of the commutative field with∑ x, P.y0 xas an opaque atom.ringtreats division as multiplication by an inverse, so it does not need — which matters, because there is no such hypothesis here.Goal⊢ 1 / ↑n * -∑ x, P.y0 x = -((∑ x, P.y0 x) / ↑n)rw [expect_neverTreat, HT_neverTreat_empty, twoUnits_mean_y0, twoUnits_tau]; norm_numThe counterexample itself. The three rewrites turn
(neverTreat 2).expect (HT (neverTreat 2) twoUnits) ≠ twoUnits.tauinto-0 ≠ 1:expect_neverTreatevaluates the expectation at∅,HT_neverTreat_emptyreplaces the estimator by-mean y0, and the twotwoUnitslemmas put in the numbers.norm_numthen decides the numeric disequality. Rewriting inside a≠is legitimate becausea ≠ bis notation fora = b → False, andrwrewrites 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 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 hhi is not vacuous either —
Design.propensity_le_one says < 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
By linearity of expectation it is enough to compute the expectation of a single
summand, since
Now
∎
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
the individual treatment effect of that unit.
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.
have hne : D.propensity i ≠ 0 := ne_of_gt hloConvert the hypothesis
into the form the cancellation lemma wants.ne_of_gt : a < b → b ≠ a— the direction matters, and it is the reason this ishne : D.propensity i ≠ 0rather than0 ≠ D.propensity i.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 hhigives1 ≠ D.propensity ifrom , andsub_ne_zero_of_ne : a ≠ b → a - b ≠ 0turns that into1 - π i ≠ 0. These twohaves are the only placehloandhhiare used in this proof: everything before the final cancellation holds for any propensity at all.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
rwlist: re-associateZ z i * y1 i / π iintoZ z i * (y1 i / π i), so that everything except the indicator is a single constant factor.Design.expect_congrsays two statistics that agree at everyzhave the same expectation, and it is how you rewrite under the binder — you cannotrwinside afun z => …directly, so this lemma is the standard workaround. Its argument is the pointwise proof, discharged bymul_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 iD.expect_subLinearity: 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.D.expect_mul_const, D.expect_mul_constTwice, once per half:
E[f z * c] = E[f] * cpulls the constantsy1 i / π iandy0 i / (1 - π i)out of the expectations, leaving onlyE[Z i]andE[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 iDesign.expect_Z, Design.expect_one_sub_ZThe indicator expectation step:
and . The first isrfl— true by definition, becausepropensityis defined asE[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 imul_div_cancel₀ _ hne, mul_div_cancel₀ _ hne'The cancellations
and .mul_div_cancel₀takes the proof that the denominator is non-zero as its second argument, which is whathneandhne'are for. After both, the goal readsy1 i - y0 i = y1 i - y0 iand therwcloses it.
Step 2: sum over units. With the per-unit fact in hand, the theorem is bookkeeping.
have hexpand : … := D.expect_const_mul _ (1 / (n : ℝ))Pull the leading
outside the expectation, which is the ” is not random” half of the informal first sentence.haveintroduces 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 lemmaexpect_const_mul : E[c * f] = c * E[f].D.expect_sum Finset.univ (fun i z => …)Exchange
and — the other half of the first sentence. This is Fubini for finite sums, and it is proved byFinset.sum_comm, i.e. by swapping the order of a double sum. Nothing deeper is happening. The explicitfun 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 ofz.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.tauFinset.sum_congr rfl fun i _ => expect_HT_term D P i (hlo i) (hhi i)Apply step 1 inside the sum, once per unit.
rflsays the index set is unchanged; the function argument supplies, for eachi, the proof that its summand equalsy1 i - y0 i, with the two pointwise hypotheses obtained by instantiatinghloandhhiat that unit.Goal⊢ 1 / ↑n * ∑ i, (P.y1 i - P.y0 i) = P.tauP.tau_eqUnfold 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)) / ↑nringClose it.
ringproves any identity that holds in a commutative ring, treating division as multiplication by an inverse, so it does not need to know that : and are the same term after normalization whatever 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:
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)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),Pis the population, and the twofuns are the pointwise hypotheses∀ i, 0 < π iand∀ i, π i < 1thatHT_unbiaseddemands — supplied as functions of the unit, because that is what a∀ ihypothesis is.fun i => by rw [bernoulli_propensity]; exact hloThe first hypothesis, one unit at a time. The goal inside the
byis0 < (bernoulliDesign n p hp0 hp1).propensity i;bernoulli_propensity : π i = prewrites it to0 < p, which is exactly the assumptionhlo. Abyblock inside a term is an ordinary tactic proof of whatever that position requires — here, of a proposition about one unit.fun i => by rw [bernoulli_propensity]; exact hhiThe second, identically:
π i < 1becomesp < 1, which ishhi. 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
Formally.
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 hltR0 < n₁ and n₁ < n are needed for the obvious reason. At neverTreat above, written differently — and DiM
returns -mean y0. At DiM
returns mean y1. Neither is 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-
Under complete randomization every unit has the same propensity
Fix an assignment
the treated sample mean, and likewise the control half is
Complete randomization assigns probability zero to every assignment with
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
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
refine Finset.sum_congr rfl fun z _ => ?_Both expectations are sums over all assignments, so compare them one assignment at a time.
rflis “same index set”, and the remaining goal isD.prob z * f z = D.prob z * g zfor a fixedz. Note what this step has bought: the quantifier is gone andzis now a constant the hypothesishcan be applied to.by_cases hz : D.prob z = 0The case split the lemma exists for. The hypothesis
hsays nothing when the probability is zero, so that case has to be handled by arithmetic instead.rw [hz, zero_mul, zero_mul]The zero branch. Rewriting with
hzmakes both sides0 * …, andzero_multurns each into0. No hypothesis aboutforgis needed, which is the entire point of the lemma: off the support the two statistics may disagree wildly and the expectation cannot tell.rw [h z hz]The other branch, whose state on entry is the one shown below. Now
hz : ¬D.prob z = 0is exactly whathwants, soh z hz : f z = g zrewrites the left-hand side into the right and the goal closes by reflexivity.Goalz : 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
The pmf is
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])by_contra hcThe 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. Afterby_contrathe goal isFalseand bothhzandhcare available.Goalhc : ¬#z = n₁ ⊢ Falseexact hz (by simp [completeRandomization_prob, hc])hz : D.prob z ≠ 0is, unfolded, a function fromD.prob z = 0toFalse, so applying it to a proof that the probability is zero closes the goal. That proof is the innerby simp:completeRandomization_probrewritesD.prob ztoif #z = n₁ then (n.choose n₁ : ℝ)⁻¹ else 0, andhcsends 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
The statistical content is one sentence: when every propensity is the same
constant, inverse-probability weighting is just averaging. With
the treated sample mean, since the
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.
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 to , whereomega— a decision procedure for linear arithmetic overℕandℤ— getsn ≠ 0out ofhpos : 0 < n₁andhle : n₁ ≤ n.hn1R : (n₁ : ℝ) ≠ 0is the same line.have hn0R : (n : ℝ) - (n₁ : ℝ) ≠ 0 := by …The third, and the one that has to be done in
: casthltto(n₁ : ℝ) < (n : ℝ)withexact_mod_castand letlinarithconclude(n : ℝ) - (n₁ : ℝ) ≠ 0. Doing it as((n - n₁ : ℕ) : ℝ) ≠ 0instead would be a different, weaker statement, because natural subtraction truncates.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_complis chapter 6’s#zᶜ = n - #z, which lives inℕwhere subtraction is truncated (3 - 5 = 0);hzreplaces#zbyn₁; andNat.cast_sub hleis what licenses moving the now-ℕsubtraction across the cast, which it does only becausen₁ ≤ n. This is the single most common source of wrong autoformalized statements about finite populations, and it is why the hypothesishleis threaded everywhere.intro iNow the algebra on one summand, stated as
htermfor alliso that it can later be applied under the sum’s binder.intro ifixes the unit; the goal is the -th summand of Horvitz–Thompson with already substituted, against the -th summand of the difference-in-means form.Goali : 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₁)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 / bturns dividing by into multiplying by and dividing by , andone_sub_div, whose hypothesis ishnR, rewrites1 - n₁ / nto(n - n₁) / n. Both moves are exactly the two sentences of the informal proof.Goali : 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₁)field_simpClear the denominators and finish the summand.
field_simpneeds to know that each denominator is non-zero, and it findshnR,hn1Randhn0Rin the context — which is the only reason those threehaves are at the top of the proof.unfold DiM HTBack 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ᶜand1 / n * ∑ i, (…).simp only [completeRandomization_propensity]Substitute
. It has to be asimp onlyrather than arw, because the propensity occurs under the sum’s binder, whererwcannot reach;simp onlyrewrites everywhere, including insidefuns 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))Finset.mul_sumFirst entry of the closing
rwlist:c * ∑ i, f i = ∑ i, c * f imoves the1 / ninside the sum, so thathterm— which is a statement about one summand, with the1 / nincluded — becomes applicable.Finset.sum_congr rfl fun i _ => hterm iApply the summand identity to every unit. The right-hand side is now the sum of the difference-in-means summands, with the constants
and 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₁))Finset.sum_sub_distrib, ← Finset.sum_div, ← Finset.sum_divSplit the sum of differences into a difference of sums, then pull each constant denominator back out of its sum — the
←runsFinset.sum_div : (∑ i, f i) / c = ∑ i, f i / cbackwards, 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 , minus another over .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₁)sum_mem_eq_sum_Z_mul, sum_compl_eq_sum_one_sub_Z_mul, hz, hcomplThe last four rewrites convert the left side.
sum_mem_eq_sum_Z_mulis∑ i ∈ z, v i = ∑ i, Z z i * v iandsum_compl_eq_sum_one_sub_Z_mulis 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. Thenhzandhcomplreplace the two group sizes byn₁and(n : ℝ) - (n₁ : ℝ), and the two sides are syntactically equal.
With that, DiM_unbiased_completeRandomization has nothing left to do.
have hnR : (0 : ℝ) < (n : ℝ) := by exact_mod_cast Nat.lt_of_lt_of_le hpos hleThe three numeric facts the last two lines will need, all obtained by casting the natural-number hypotheses into
: (from ), thenhn1R : 0 < (n₁ : ℝ)andhltR : (n₁ : ℝ) < (n : ℝ)byexact_mod_castfromhposandhlt.exact_mod_castcloses a goal that differs from a hypothesis only by casts.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_zeroextracts#z = n₁and feeds it toDiM_eq_HT_of_card, which returns the identity at that assignment. This onehaveis the whole “agree on the support” step.rw [(completeRandomization n n₁ hle).expect_congr_of_support hsupport]Transfer the expectation. The goal changes from
expect (DiM P) = P.tautoexpect (HT … P) = P.tau, and no sum over assignments has been touched.refine HT_unbiased _ P (fun i => ?_) (fun i => ?_)Quote the theorem.
refineapplies it but leaves the two pointwise hypotheses as goals — that is what the?_holes mean — one for0 < π iand one forπ i < 1, each for a general uniti.rw [completeRandomization_propensity]; exact div_pos hn1R hnRThe first hole. Rewriting the propensity to
n₁ / nreduces it to0 < (n₁ : ℝ) / (n : ℝ), which isdiv_posapplied to the two positivity facts from the first step.rw [completeRandomization_propensity]; exact (div_lt_one hnR).mpr hltRThe second hole.
div_lt_one hnR : a / b < 1 ↔ a < b(for0 < b), so the goaln₁ / n < 1follows fromhltR : (n₁ : ℝ) < (n : ℝ). This is the step that would fail at , and it is the only placehltis used in this theorem.
How wrong can the weights be?
HT_unbiased needs the true propensity scores. Suppose an analyst uses some
other number
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]
ringProposition. Let the design have constant propensity
By linearity, and because
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_termisexpect_HT_term’s proof with the twohaves 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.D.expect_sub, D.expect_mul_const, D.expect_mul_const, Design.expect_Z, Design.expect_one_sub_ZLinearity, then the two constants out, then the two indicator expectations — and stop. The goal is now closed, because the statement of
expect_HTq_termkeepsπ 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 .rw [hexpand, D.expect_sum Finset.univ (fun i z => …), Finset.sum_congr rfl fun i _ => expect_HTq_term D P q i]expect_HTqnow repeats step 2 ofHT_unbiasedverbatim:hexpandpulls the1 / nout withexpect_const_mul,expect_sumexchanges and , andsum_congrapplies 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.y0simp only [hp]Use the constancy of the design’s propensity:
hp : ∀ i, D.propensity i = preplacesD.propensity ibypunder the binder, whichrwcould 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.y0rw [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 ringis an inline anonymous lemma: it states the summand identity and proves it byring. Writing it inline avoids ahavewhose 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.y0Finset.sum_sub_distrib, ← Finset.mul_sum, ← Finset.mul_sumSplit the sum and pull both constant factors out —
← Finset.mul_sumrunsc * ∑ i, f i = ∑ i, c * f ifrom 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.y0simp only [mean]; ringUnfold
meanon the right-hand side — the state below — so that both sides are built from∑ i, P.y1 i,∑ i, P.y0 iand the constants;ringthen matches them. As inHT_unbiased, no hypothesis is needed, becauseringnever 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 twoUnits (where
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_numnoncomputable 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_numarguments are the proofs of0 ≤ 1/2and1/2 ≤ 1thatbernoulliDesignrequires — a definition can take proofs as arguments, and here they are inlined rather than named because nothing else will ever refer to them.theorem halfBernoulli_propensity (i : Fin 2) : halfBernoulli.propensity i = 1 / 2 := bernoulli_propensity iA term-mode restatement of
bernoulli_propensityat this design, marked@[simp]so laternorm_numcalls can use it. It type-checks becausehalfBernoulliis definitionally that Bernoulli design, not merely equal to it.rw [expect_HTq halfBernoulli twoUnits (1 / 4) halfBernoulli_propensity, twoUnits_mean_y1, twoUnits_mean_y0]; norm_numThe computation. The general formula is applied with
(supplied implicitly by thehpargumenthalfBernoulli_propensity) and , then the two population means are replaced by1and0, leaving(1/2) / (1/4) * 1 - (1 - 1/2) / (1 - 1/4) * 0 = 2fornorm_numto check.rw [expect_HTq_halfBernoulli, twoUnits_tau]; norm_numAnd the bias, recorded as a disequality: the two rewrites reduce the goal to
(2 : ℝ) ≠ 1, whichnorm_numdecides. 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
On the first, DiM returns -mean y0; on the second, mean y1. Neither has
any reason to equal
Two units, 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
There are four assignments and each has probability
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_numrw [show (univ : Finset (Assignment 2)) = {∅, {0}, {1}, {0, 1}} from by decide]sum_assignments_twois the whole trick.Assignment 2isFinset (Fin 2), which has four elements, anddecideproves the enumeration by evaluating both sides — there is no cleverness to it, just a finite check the kernel is willing to run. Wrapping it inshow … fromsupplies the statement, sincedecideneeds a closed proposition to evaluate.repeat rw [Finset.sum_insert (by decide)]Peel the four terms off one at a time.
Finset.sum_insertrequires that the head is not already in the tail, and that side condition is decidable too, which is whyby decidedischarges it;repeatapplies the rewrite until it no longer matches.Finset.sum_singletonthen handles the last element andringreassociates the additions.Goal⊢ f ∅ + (f {0} + (f {1} + ∑ x ∈ {{0, 1}}, f x)) = f ∅ + f {0} + f {1} + f {0, 1}rw [DiM]; simp only [twoUnits_y1, twoUnits_y0, Finset.sum_const, nsmul_eq_mul, mul_one, mul_zero, zero_div, sub_zero]DiM_twoUnitsnext, and thishaveis where the population’s flatness is spent. Both outcome vectors are constant, soFinset.sum_constturns∑ i ∈ z, (1 : ℝ)into#z • 1andnsmul_eq_mulinto(#z : ℝ) * 1;mul_onecleans that up,mul_zeroandzero_divannihilate the control half, andsub_zerodeletes it. What is left is the claim that the estimator is(#z : ℝ) / (#z : ℝ)— the treated sample mean of a vector of ones.rw [hz]; splitThe goal is now an
if, sosplitdoes a case analysis on its condition, naming the two branch hypotheses withnext h => …. The point of the whole lemma is thatx / xis not unconditionally1.Goal⊢ ↑(#z) / ↑(#z) = if z = ∅ then 0 else 1next h => simp only [h, Finset.card_empty, Nat.cast_zero, div_zero]The empty branch, where
h : z = ∅. Rewriting with it makes the cardinality#(∅), whichFinset.card_emptysends to0andNat.cast_zeroto(0 : ℝ), and thendiv_zerogives0 / 0 = 0. This is the assignment that costs the estimator its unbiasedness.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_selfneeds the denominator non-zero, and the argument is that chain read inside out: non-empty fromhbyFinset.nonempty_iff_ne_empty, hence#z ≠ 0byFinset.card_ne_zero, hence(#z : ℝ) ≠ 0byNat.cast_ne_zero. Three lemmas to move one fact across two type boundaries is typical.rw [Design.expect, sum_assignments_two (fun z => halfBernoulli.prob z * DiM twoUnits z)]expect_DiM_halfBernoulliputs the pieces together. Unfoldingexpectexposes the sum over assignments, and the enumeration lemma — applied to the specific summandprob 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 / 4simp only [DiM_twoUnits, halfBernoulli, bernoulliDesign_prob]; norm_numEvaluate.
DiM_twoUnitsreplaces each estimator value by itsif,bernoulliDesign_probreplaces each probability byp ^ #z * (1 - p) ^ (n - #z)withp = 1/2, andnorm_numcomputes the four cardinalities, the four conditionals and the arithmetic, arriving at3 / 4. The finalDiM_biased_halfBernoulliis then two rewrites and anorm_numon(3 : ℝ) / 4 ≠ 1.
So div_zero guesses right by luck. That is what a bias of exactly
Three honest repairs, in increasing order of usefulness:
- State a different estimand.
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 for every . - Condition on the realized count. The Bernoulli pmf is constant on
— every such assignment has probability — so the Bernoulli design conditioned on is exactly complete randomization with treated units. Hence for every , and by the law of total expectation 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
HTinstead.HT_unbiased_bernoulliholds unconditionally for , because Horvitz–Thompson divides by the design’s , which is not random, rather than by the realized , 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:
#print axioms HT_unbiased
#print axioms DiM_unbiased_completeRandomizationThree names come back, and they are the same three that come back for almost every theorem in Mathlib:
propext— propositional extensionality: ifp ↔ qthenp = 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 whyby_casesandClassical.emare 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,Finsetis built onMultisetwhich is a quotient ofList. 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.
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
sorryHint
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
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 y0, so DiM P ∅ = -mean y0. No hypothesis is
required, which is the point: the statement is about every assignment.
theorem ex_DiM_empty (P : Population n) : DiM P ∅ = -mean P.y0 := by
sorryHint
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
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.
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
sorryHint
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
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
theorem ex_bernoulli_prob_empty_pos {p : ℝ} (hp0 : 0 ≤ p) (hp1 : p ≤ 1) (hlt : p < 1) :
0 < (bernoulliDesign n p hp0 hp1).prob ∅ := by
sorryHint
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
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) nExerciseThe 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.
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
sorryHint
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
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]
ringReading 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.