Panel.Weighted

lemma ip_self_eq_zero_iff (c : WeightedSupport R) (A : R → ℝ) : c.ip A A = 0 ↔ ∀ r ∈ c.observed, A r = 0 := by unfold ip rw [Finset.sum_eq_zero_iff_of_nonneg (c.ip_self_summand_nonneg A)] constructor · -- From `ω_r · A_r · A_r = 0` and `ω_r > 0`, deduce `A_r = 0`. intro h r hr have hwpos : 0 < c.weight r := c.weight_pos r hr have hsum : c.weight r * A r * A r = 0 := h r hr have hsum' : c.weight r * (A r * A r) = 0 := by rw [← mul_assoc]; exact hsum have hne : c.weight r ≠ 0 := ne_of_gt hwpos have hsq : A r * A r = 0 := (mul_eq_zero.mp hsum').resolve_left hne exact mul_self_eq_zero.mp hsq · -- `A_r = 0` ⇒ `ω_r · A_r · A_r = 0`. intro h r hr have hAr : A r = 0 := h r hr rw [hAr]; ring

lemma residualize_in_orthogonal (c : WeightedSupport R) (H : Submodule ℝ (R → ℝ)) (X : R → ℝ) {h : R → ℝ} (hH : h ∈ H) : c.ip (c.tildeX H X) h = 0 := by simp only [tildeX_eq] exact c.proj_orthogonal H X hH

theorem fwl_identity (c : WeightedSupport R) (H : Submodule ℝ (R → ℝ)) (X : Fin K → R → ℝ) (Y : R → ℝ) (hRank : c.RankCondition H X) : ∀ β : Fin K → ℝ, ∀ α : R → ℝ, α ∈ H → (∀ β' : Fin K → ℝ, ∀ α' : R → ℝ, α' ∈ H → c.ip (Y - (∑ k, β k • X k) - α) (Y - (∑ k, β k • X k) - α) ≤ c.ip (Y - (∑ k, β' k • X k) - α') (Y - (∑ k, β' k • X k) - α')) → β = thetaHat c X H Y := by intro β α hα hmin -- Step 1: Set up the residual `R₀ := Y - ∑ k β k • X k - α`. set R₀ : R → ℝ := Y - (∑ k, β k • X k) - α with hR₀_def -- Step 2: Derive joint orthogonality from joint minimality. -- -- (a) Varying α while fixing β: for any `g ∈ H`, the test point -- `(β, α + t • g) ∈ Fin K → ℝ × H` gives a one-parameter family. -- The residual at the test point is `R₀ - t • g`. -- Minimality gives `⟨R₀, g⟩_ω = 0` for all `g ∈ H`. have h_orth_H : ∀ g ∈ H, c.ip R₀ g = 0 := by intro g hg -- Apply the WLS first-order condition: `R₀ = Y - fitted β X - α`, -- treat as `(Y - fitted β X)` projected on `H` with candidate `α`. -- Specialize hmin to varying α only. have hmin_alpha : ∀ α' ∈ H, c.ip (Y - (∑ k, β k • X k) - α) (Y - (∑ k, β k • X k) - α) ≤ c.ip (Y - (∑ k, β k • X k) - α') (Y - (∑ k, β k • X k) - α') := by intro α' hα' exact hmin β α' hα' -- Use residualize_orth_iff_argmin with X := Y - fitted β X, p := α. have := (c.residualize_orth_iff_argmin H (Y - (∑ k, β k • X k)) α hα).mpr hmin_alpha exact this g hg -- (b) Varying β_j while fixing α: for any `j : Fin K` and `t : ℝ`, -- the test point `(β + t • e_j, α)` gives residual `R₀ - t • X j`. -- The minimality + perturbation argument forces `⟨R₀, X j⟩_ω = 0`. have h_orth_X : ∀ j : Fin K, c.ip R₀ (X j) = 0 := by intro j -- Build the perturbation: replace β by β + t • e_j. -- Specialize hmin to such perturbations. have hkey : ∀ t : ℝ, c.ip R₀ R₀ ≤ c.ip (R₀ - t • X j) (R₀ - t • X j) := by intro t -- β' := β + t • (Pi.single j 1). Then `∑ k β' k • X k = (∑ k β k • X k) + t • X j`. let δ : Fin K → ℝ := Pi.single j t let β' : Fin K → ℝ := β + δ have hsum_eq : (∑ k, β' k • X k) = (∑ k, β k • X k) + t • X j := by have hsplit : ∀ k, β' k • X k = β k • X k + δ k • X k := by intro k; simp only [β', Pi.add_apply, add_smul] rw [Finset.sum_congr rfl (fun k _ => hsplit k)] rw [Finset.sum_add_distrib] congr 1 rw [Finset.sum_eq_single j] · simp [δ, Pi.single_eq_same] · intro k _ hk simp [δ, Pi.single_eq_of_ne hk] · intro h; exact absurd (Finset.mem_univ j) h have hmin' := hmin β' α hα have heq : Y - (∑ k, β' k • X k) - α = R₀ - t • X j := by rw [hsum_eq, hR₀_def] ext s; simp; ring rw [heq] at hmin' exact hmin' -- Now apply the perturbation argument: quadratic-in-t inequality. have hquad : ∀ t : ℝ, 0 ≤ - (2 * t * c.ip R₀ (X j)) + t^2 * c.ip (X j) (X j) := by intro t have := hkey t have hexp := c.ip_sub_smul_expand R₀ 0 (X j) t -- ip_sub_smul_expand: `c.ip (X - p - t • h) ... = ⟨X-p,X-p⟩ - 2t⟨X-p,h⟩ + t²⟨h,h⟩`. -- Specialize at X := R₀, p := 0, h := X j. simp only [sub_zero] at hexp linarith -- Same case analysis as in WLS. by_cases hzero : c.ip (X j) (X j) = 0 · -- X j vanishes on observed, so `⟨R₀, X j⟩_ω = 0` directly. have hvan : ∀ r ∈ c.observed, X j r = 0 := (c.ip_self_eq_zero_iff (X j)).mp hzero exact c.ip_eq_zero_of_zero_on_observed R₀ (X j) hvan · have hpos : 0 < c.ip (X j) (X j) := lt_of_le_of_ne (c.ip_self_nonneg (X j)) (Ne.symm hzero) set a : ℝ := c.ip R₀ (X j) with ha set b : ℝ := c.ip (X j) (X j) with hb have htest := hquad (a / b) have hsimp : - (2 * (a / b) * a) + (a / b)^2 * b = - a^2 / b := by field_simp; ring rw [hsimp] at htest have ha2 : a^2 ≤ 0 := by have : - a^2 ≥ 0 := by have := (div_nonneg_iff.mp htest).resolve_right ?_ · exact this.1 · push_neg; intro _; exact hpos linarith have ha2nn : 0 ≤ a^2 := sq_nonneg a have ha2eq : a^2 = 0 := le_antisymm ha2 ha2nn exact sq_eq_zero_iff.mp ha2eq -- Step 3: Derive the normal equations `Q_XX β = rhsVec`. -- For each k: `⟨tildeX (X k), Y - ∑ j β j X j⟩_ω = 0` -- (since `tildeX (X k) ⊥ H` absorbs α, and `tildeX (X k) ⊥ ?`). have h_normal : ∀ k : Fin K, ∑ j, β j * c.ip (c.tildeX H (X j)) (c.tildeX H (X k)) = c.ip (c.tildeX H (X k)) Y := by intro k -- Use `ip_tildeX_eq_ip_tildeX_residual` to rewrite both sides into the -- tildeX-tildeX form, then apply orthogonality of residuals. -- Strategy: from h_orth_X applied to j-coordinate, we have -- c.ip R₀ (X j) = 0. -- Expand R₀ = Y - (∑ β j • X j) - α. -- So `c.ip Y (X k) - ∑ j β j • c.ip (X j) (X k) - c.ip α (X k) = 0`. -- But we want to work with tildeX. Use the lemma: residualizing one side -- of an inner product equals taking ⟨tildeX, tildeX⟩. -- -- Cleaner route: compute c.ip R₀ (tildeX H (X k)). -- Since tildeX H (X k) ⊥ H, c.ip α (tildeX H (X k)) = 0. -- Since tildeX H (X k) = X k - proj H (X k), and `c.ip R₀ g = 0` for g ∈ H: -- c.ip R₀ (tildeX H (X k)) = c.ip R₀ (X k) - c.ip R₀ (proj H (X k)) -- = 0 - 0 = 0. have hR₀_tilde : c.ip R₀ (c.tildeX H (X k)) = 0 := by -- tildeX H (X k) = X k + (- proj H (X k)) have hsplit : c.tildeX H (X k) = X k + (- c.proj H (X k)) := by simp [tildeX_eq, sub_eq_add_neg] rw [hsplit, c.ip_add_right] have h2 : c.ip R₀ (-c.proj H (X k)) = - c.ip R₀ (c.proj H (X k)) := by have heq : (-c.proj H (X k) : R → ℝ) = (-1 : ℝ) • c.proj H (X k) := by ext s; simp rw [heq, c.ip_smul_right]; ring rw [h_orth_X k, h2, h_orth_H (c.proj H (X k)) (c.proj_mem H (X k))]; ring -- Now expand c.ip R₀ (tildeX H (X k)) using bilinearity. have h_expand : c.ip R₀ (c.tildeX H (X k)) = c.ip Y (c.tildeX H (X k)) - ∑ j, β j * c.ip (X j) (c.tildeX H (X k)) - c.ip α (c.tildeX H (X k)) := by have hR₀_split : R₀ = Y - (∑ j, β j • X j) - α := hR₀_def have hsubeq : Y - (∑ j, β j • X j) - α = Y + (-(∑ j, β j • X j)) + (-α) := by ext s; simp [sub_eq_add_neg] rw [hR₀_split, hsubeq] rw [c.ip_add_left, c.ip_add_left] -- Term 1: c.ip Y (tildeX H (X k)) -- Term 2: c.ip (-(∑ j, β j • X j)) (tildeX H (X k)) -- = -c.ip (∑ j, β j • X j) (tildeX H (X k)) -- = -∑ j, β j * c.ip (X j) (tildeX H (X k)) -- Term 3: c.ip (-α) (tildeX H (X k)) = -c.ip α (tildeX H (X k)) have hneg_sum : c.ip (-(∑ j, β j • X j)) (c.tildeX H (X k)) = -∑ j, β j * c.ip (X j) (c.tildeX H (X k)) := by have heq : (-(∑ j, β j • X j) : R → ℝ) = (-1 : ℝ) • (∑ j, β j • X j) := by ext s; simp rw [heq, c.ip_smul_left, c.ip_sum_left] have hinner : ∀ j, c.ip (β j • X j) (c.tildeX H (X k)) = β j * c.ip (X j) (c.tildeX H (X k)) := by intro j; rw [c.ip_smul_left] rw [Finset.sum_congr rfl (fun j _ => hinner j)] ring have hneg_α : c.ip (-α) (c.tildeX H (X k)) = -c.ip α (c.tildeX H (X k)) := by have : (-α : R → ℝ) = (-1 : ℝ) • α := by ext s; simp rw [this, c.ip_smul_left]; ring rw [hneg_sum, hneg_α]; ring -- c.ip α (tildeX H (X k)) = 0 because tildeX ⊥ H and α ∈ H. have hα_tilde : c.ip α (c.tildeX H (X k)) = 0 := by rw [c.ip_symm] exact c.residualize_in_orthogonal H (X k) hα -- c.ip (X j) (tildeX H (X k)) = c.ip (tildeX H (X j)) (tildeX H (X k)). have hX_tilde : ∀ j, c.ip (X j) (c.tildeX H (X k)) = c.ip (c.tildeX H (X j)) (c.tildeX H (X k)) := by intro j rw [c.ip_symm (X j), c.ip_symm (c.tildeX H (X j))] exact c.ip_tildeX_eq_ip_tildeX_residual H X (X j) k -- c.ip Y (tildeX H (X k)) = c.ip (tildeX H (X k)) Y. have hY_symm : c.ip Y (c.tildeX H (X k)) = c.ip (c.tildeX H (X k)) Y := c.ip_symm _ _ -- Substitute everything into h_expand combined with hR₀_tilde = 0. rw [hα_tilde, sub_zero] at h_expand rw [hY_symm] at h_expand have : ∑ j, β j * c.ip (X j) (c.tildeX H (X k)) = ∑ j, β j * c.ip (c.tildeX H (X j)) (c.tildeX H (X k)) := by refine Finset.sum_congr rfl ?_ intro j _; rw [hX_tilde j] rw [this] at h_expand linarith [hR₀_tilde] -- Step 4: convert h_normal to `Q_XX.mulVec β = rhsVec`. have hmat : (Q_XX c X H).mulVec β = rhsVec c X H Y := by funext k simp only [mulVec, dotProduct, Q_XX_apply, rhsVec_apply] -- dotProduct: ∑ j, Q_XX k j * β j = ∑ j, c.ip (tildeX (X k)) (tildeX (X j)) * β j have hk := h_normal k -- h_normal: ∑ j, β j * c.ip (tildeX (X j)) (tildeX (X k)) = c.ip (tildeX (X k)) Y. -- Need: ∑ j, c.ip (tildeX (X k)) (tildeX (X j)) * β j = c.ip (tildeX (X k)) Y. -- These match by symmetry of c.ip and commutativity of multiplication. rw [← hk] refine Finset.sum_congr rfl ?_ intro j _ rw [c.ip_symm (c.tildeX H (X k)) (c.tildeX H (X j))] ring -- Step 5: β = Q_XX⁻¹ (Q_XX β) = Q_XX⁻¹ rhsVec = thetaHat. unfold thetaHat have : (Q_XX c X H)⁻¹.mulVec ((Q_XX c X H).mulVec β) = (Q_XX c X H)⁻¹.mulVec (rhsVec c X H Y) := by rw [hmat] rw [Matrix.mulVec_mulVec] at this rw [Matrix.nonsing_inv_mul _ hRank, Matrix.one_mulVec] at this exact this

lemma cellMean_mul_cellMass {𝒢 : Type*} [DecidableEq 𝒢] (c : WeightedSupport R) (F : R → ℝ) (G : R → 𝒢) (g : 𝒢) (hpos : 0 < c.cellMass G g) : c.cellMean F G g hpos * c.cellMass G g = ∑ r ∈ c.observed, c.weight r * F r * cellIndicator G g r := by unfold cellMean exact div_mul_cancel₀ _ (ne_of_gt hpos)

lemma weighted_center_cov {K : Type*} [Field K] (p x y : ι → K) (hp : ∑ i, p i = 1) (h2 : (2 : K) ≠ 0) : ∑ i, p i * (x i - ∑ j, p j * x j) * (y i - ∑ j, p j * y j) = (1 / 2) * ∑ i, ∑ j, p i * p j * (x i - x j) * (y i - y j) := by rw [weighted_center_cov_left p x y hp, pairwise_cov_right p x y hp] simp [h2]

theorem ip_eq_integral (μ : Measure R) [IsProbabilityMeasure μ] (A B : R → ℝ) : (ofProbabilityMeasure μ).ip A B = ∫ r, A r * B r ∂μ := by classical -- The integral over a finite measurable space is the atom-weighted sum. rw [MeasureTheory.integral_fintype (Integrable.of_finite)] -- `ip` sums over `observed`; extend to all of `R` since off-support weights -- are zero. unfold ip rw [Finset.sum_subset (Finset.subset_univ _)] · refine Finset.sum_congr rfl ?_ intro r _ simp only [ofProbabilityMeasure_weight, Measure.real, smul_eq_mul] ring · intro r _ hr have hzero : (μ {r}).toReal = 0 := by by_contra hne have hpos : 0 < (μ {r}).toReal := lt_of_le_of_ne ENNReal.toReal_nonneg (Ne.symm hne) exact hr (Finset.mem_filter.mpr ⟨Finset.mem_univ r, hpos⟩) simp [ofProbabilityMeasure_weight, hzero]

theorem scalar_fwl_of_normalEqs (c : WeightedSupport R) (H : Submodule ℝ (R → ℝ)) (X Y : R → ℝ) (β : ℝ) (α : R → ℝ) (hα : α ∈ H) (hden : c.ip (c.tildeX H X) (c.tildeX H X) ≠ 0) (h_normal_X : c.ip (Y - β • X - α) X = 0) (h_normal_H : ∀ h ∈ H, c.ip (Y - β • X - α) h = 0) : β = c.ip (c.tildeX H X) Y / c.ip (c.tildeX H X) (c.tildeX H X) := by set Xt : R → ℝ := c.tildeX H X with hXt set R₀ : R → ℝ := Y - β • X - α with hR₀ -- The residualized regressor is orthogonal to the regression residual. have hXt_R₀ : c.ip Xt R₀ = 0 := by have hsplit : Xt = X + (-1 : ℝ) • c.proj H X := by rw [hXt, tildeX_eq]; ext r; simp [sub_eq_add_neg] rw [hsplit, c.ip_add_left] -- ⟨X, R₀⟩ = ⟨R₀, X⟩ = 0 have h1 : c.ip X R₀ = 0 := by rw [c.ip_symm]; exact h_normal_X -- ⟨(-1)•proj X, R₀⟩ = (-1) * ⟨R₀, proj X⟩ = 0 have h2 : c.ip ((-1 : ℝ) • c.proj H X) R₀ = 0 := by rw [c.ip_smul_left, c.ip_symm] rw [h_normal_H (c.proj H X) (c.proj_mem H X)] ring rw [h1, h2, add_zero] -- Expand ⟨Xt, R₀⟩ = ⟨Xt, Y⟩ − β⟨Xt, X⟩ − ⟨Xt, α⟩, and kill ⟨Xt, α⟩. have hα_zero : c.ip Xt α = 0 := c.residualize_in_orthogonal H X hα have hXt_X : c.ip Xt X = c.ip Xt Xt := ip_tildeX_self c H X have hexpand : c.ip Xt R₀ = c.ip Xt Y - β * c.ip Xt Xt := by have hrw : R₀ = Y - β • X - α := hR₀ rw [hrw] -- Y - β•X - α = Y + (-(β•X)) + (-α) have hsub : Y - β • X - α = Y + (-(β • X)) + (-α) := by ext r; simp [sub_eq_add_neg] rw [hsub, c.ip_add_right, c.ip_add_right] have hneg_smul : c.ip Xt (-(β • X)) = - (β * c.ip Xt X) := by have : (-(β • X) : R → ℝ) = (-β) • X := by ext r; simp rw [this, c.ip_smul_right]; ring have hneg_α : c.ip Xt (-α) = - c.ip Xt α := by have : (-α : R → ℝ) = (-1 : ℝ) • α := by ext r; simp rw [this, c.ip_smul_right]; ring rw [hneg_smul, hneg_α, hα_zero, hXt_X]; ring -- Combine: 0 = ⟨Xt, Y⟩ − β⟨Xt, Xt⟩. rw [hexpand] at hXt_R₀ have hne : c.ip Xt Xt ≠ 0 := by simpa [hXt] using hden field_simp at hXt_R₀ ⊢ linarith [hXt_R₀]

theorem proj_eq_argmin (c : WeightedSupport R) (H : Submodule ℝ (R → ℝ)) (X h : R → ℝ) (hH : h ∈ H) : c.ip (X - c.proj H X) (X - c.proj H X) ≤ c.ip (X - h) (X - h) := by refine c.ip_self_sub_le_of_orth H X (c.proj H X) (c.proj_mem H X) ?_ hH intro g hg exact c.proj_orthogonal H X hg

theorem residualize_orth_iff_argmin (c : WeightedSupport R) (H : Submodule ℝ (R → ℝ)) (X p : R → ℝ) (hp : p ∈ H) : (∀ h ∈ H, c.ip (X - p) h = 0) ↔ (∀ h ∈ H, c.ip (X - p) (X - p) ≤ c.ip (X - h) (X - h)) := by constructor · -- Forward direction: standard expansion. intro horth h hh exact c.ip_self_sub_le_of_orth H X p hp horth hh · -- Reverse direction: perturbation argument. intro hmin h hh -- For any `t : ℝ`, `p + t • h ∈ H`. -- Note: `X - (p + t • h) = X - p - t • h`. have hkey : ∀ t : ℝ, c.ip (X - p) (X - p) ≤ c.ip (X - p - t • h) (X - p - t • h) := by intro t have hpth : p + t • h ∈ H := H.add_mem hp (H.smul_mem t hh) have := hmin (p + t • h) hpth have heq : X - (p + t • h) = X - p - t • h := by ext s; simp; ring rw [heq] at this exact this -- Combined with the quadratic expansion, this gives -- 0 ≤ -2 t · ⟨X-p, h⟩_ω + t² · ⟨h, h⟩_ω for all `t : ℝ`. have hquad : ∀ t : ℝ, 0 ≤ - (2 * t * c.ip (X - p) h) + t^2 * c.ip h h := by intro t have := hkey t have hexp := c.ip_sub_smul_expand X p h t linarith -- Case split on whether `⟨h, h⟩_ω > 0`. by_cases hzero : c.ip h h = 0 · -- `h` vanishes on observed, so `⟨X - p, h⟩_ω = 0` directly. have hvan : ∀ r ∈ c.observed, h r = 0 := (c.ip_self_eq_zero_iff h).mp hzero exact c.ip_eq_zero_of_zero_on_observed (X - p) h hvan · -- `⟨h, h⟩_ω > 0`. Plug `t = ⟨X-p, h⟩_ω / ⟨h, h⟩_ω`. have hpos : 0 < c.ip h h := lt_of_le_of_ne (c.ip_self_nonneg h) (Ne.symm hzero) set a : ℝ := c.ip (X - p) h with ha set b : ℝ := c.ip h h with hb -- The minimum of `t ↦ -2 a t + b t²` over `t : ℝ` is `-a² / b` at -- `t = a / b`. have htest := hquad (a / b) -- `- (2 * (a/b) * a) + (a/b)² * b = -2 a² / b + a² / b = -a² / b`. have hsimp : - (2 * (a / b) * a) + (a / b)^2 * b = - a^2 / b := by field_simp ring rw [hsimp] at htest -- So `-a²/b ≥ 0`. Combined with `b > 0`, this forces `a² ≤ 0`, -- hence `a = 0`. have hineq : - a^2 / b ≥ 0 := htest have ha2 : a^2 ≤ 0 := by have : - a^2 ≥ 0 := by have := (div_nonneg_iff.mp hineq).resolve_right ?_ · exact this.1 · push_neg intro _ exact hpos linarith have ha2nn : 0 ≤ a^2 := sq_nonneg a have ha2eq : a^2 = 0 := le_antisymm ha2 ha2nn have : a = 0 := by have := sq_eq_zero_iff.mp ha2eq exact this exact this