1import Mathlib23set_option linter.style.header false4set_option linter.style.longLine false56/-!7WP28 — machine-checked algebraic core of the proofs (Lean 4 + Mathlib).89What is formalized (deterministic statements only; the probabilistic O_p arguments are not):10 S1 Vandermonde moment constraints for the patterns (1,2), (1,2,4), (1,2,4,8), (1,2,3) and the constants beta_r.11 S3 Exact bias identities of J_2 and J_3 in the exponential model (Theorem 2(ii)).12 S4 Newton map and error identity in the exponential model (Theorem 2(iv)).13 S5 q-fold Newton error identity (Theorem 2(iv)) and the boundary constant algebra (Theorem 2(v)).14 S7 Newton basin: |u_j| = |u_0|^(2^j); convergence iff |u_0| < 1 (Theorem 2(iv)).15 S8 Exponent algebra of the frontier: with alpha(r,q) = 1 - 1/(2^(q+1) r) the scaled bias term is exactly 1 (Theorem 1 proof).16 S2/S6 Gamma inverse moments, E R_n, E R_n^2, Var(theta_hat_n), and the exact variance of J_2 from the Beta/Gamma17 decomposition (Theorem 2(i),(iii)); the moment values are derived from Real.Gamma_add_one.18-/192021namespace WP282223/-! ### S1 Vandermonde weights -/2425theorem vandermonde_12 : ((2:ℚ) + (-1) = 1) ∧ ((2:ℚ) * 1 + (-1) * 2 = 0) := by norm_num2627theorem vandermonde_124 :28 ((8/3:ℚ) + (-2) + 1/3 = 1) ∧ ((8/3:ℚ) * 1 + (-2) * 2 + (1/3) * 4 = 0) ∧29 ((8/3:ℚ) * 1 + (-2) * 4 + (1/3) * 16 = 0) := by norm_num3031theorem vandermonde_1248 :32 ((64/21:ℚ) - 8/3 + 2/3 - 1/21 = 1) ∧ ((64/21:ℚ) * 1 - (8/3) * 2 + (2/3) * 4 - (1/21) * 8 = 0) ∧33 ((64/21:ℚ) * 1 - (8/3) * 4 + (2/3) * 16 - (1/21) * 64 = 0) ∧34 ((64/21:ℚ) * 1 - (8/3) * 8 + (2/3) * 64 - (1/21) * 512 = 0) := by norm_num3536theorem vandermonde_123 :37 ((3:ℚ) - 3 + 1 = 1) ∧ ((3:ℚ) * 1 - 3 * 2 + 1 * 3 = 0) ∧ ((3:ℚ) * 1 - 3 * 4 + 1 * 9 = 0) := by norm_num3839/-- Leading bias coefficients beta_r = sum_j w_j s_j^r. -/40theorem beta_2 : (2:ℚ) * 1 + (-1) * 4 = -2 := by norm_num41theorem beta_3 : (8/3:ℚ) * 1 + (-2) * 8 + (1/3) * 64 = 8 := by norm_num42theorem beta_4 : (64/21:ℚ) * 1 - (8/3) * 16 + (2/3) * 256 - (1/21) * 4096 = -64 := by norm_num43theorem beta_3_pattern123 : (3:ℚ) * 1 - 3 * 8 + 1 * 27 = 6 := by norm_num4445/-- Remainder amplification constants A_r = sum_j |w_j| s_j. -/46theorem amp_2 : |(2:ℚ)| * 1 + |(-1:ℚ)| * 2 = 4 := by norm_num47theorem amp_3 : |(8/3:ℚ)| * 1 + |(-2:ℚ)| * 2 + |(1/3:ℚ)| * 4 = 8 := by norm_num4849/-! ### S3 Exact biases in the exponential model (Theorem 2(ii)) -/5051/-- E J_2 - θ0 = θ0 Σ_j w_j s_j/(m - s_j) = -2θ0/((m-1)(m-2)). -/52theorem exact_bias_J2 (θ m : ℝ) (h1 : m - 1 ≠ 0) (h2 : m - 2 ≠ 0) :53 θ * (2 * 1 / (m - 1) + (-1) * 2 / (m - 2)) = -2 * θ / ((m - 1) * (m - 2)) := by54 field_simp <;> ring5556/-- E J_3 - θ0 = 8θ0/((m-1)(m-2)(m-4)) for the pattern (1,2,4). -/57theorem exact_bias_J3 (θ m : ℝ) (h1 : m - 1 ≠ 0) (h2 : m - 2 ≠ 0) (h4 : m - 4 ≠ 0) :58 θ * ((8/3) * 1 / (m - 1) + (-2) * 2 / (m - 2) + (1/3) * 4 / (m - 4))59 = 8 * θ / ((m - 1) * (m - 2) * (m - 4)) := by60 field_simp <;> ring6162/-! ### S4 Newton identities in the exponential model (Theorem 2(iv)) -/6364/-- Newton map: x - g(x)/H(x) with g = 1/x - xbar, H = -1/x^2. -/65theorem newton_map (x xbar : ℝ) (hx : x ≠ 0) :66 x - (1 / x - xbar) / (-(1 / x ^ 2)) = 2 * x - x ^ 2 * xbar := by67 field_simp <;> ring6869/-- Error identity: N(x) - θ̂ = -(x - θ̂)^2/θ̂ with θ̂ = 1/xbar. -/70theorem newton_error (x xbar : ℝ) (h : xbar ≠ 0) :71 (2 * x - x ^ 2 * xbar) - 1 / xbar = -((x - 1 / xbar) ^ 2) / (1 / xbar) := by72 field_simp <;> ring7374/-! ### S5 q-fold Newton error and boundary constant -/7576/-- Error recursion e_{j+1} = -e_j^2 / t with t = θ̂_N. -/77noncomputable def newtonErr (t : ℝ) : ℕ → ℝ → ℝ78 | 0, e => e79 | n + 1, e => -(newtonErr t n e) ^ 2 / t8081/-- t^(2^q) · e_q = -t · e_0^(2^q) for q ≥ 1, i.e. e_q = -e_0^(2^q)/t^(2^q - 1). -/82theorem newtonErr_closed (t e : ℝ) (ht : t ≠ 0) :83 ∀ q : ℕ, t ^ (2 ^ (q + 1)) * newtonErr t (q + 1) e = -t * e ^ (2 ^ (q + 1)) := by84 intro q85 induction q with86 | zero =>87 simp [newtonErr]88 field_simp <;> ring89 | succ n ih =>90 have h2 : (2:ℕ) ^ (n + 2) = 2 ^ (n + 1) * 2 := by ring91 rw [newtonErr, h2, pow_mul, pow_mul]92 have : (t ^ 2 ^ (n + 1)) ^ 2 * (-(newtonErr t (n + 1) e) ^ 2 / t)93 = -((t ^ 2 ^ (n + 1) * newtonErr t (n + 1) e) ^ 2) / t := by94 field_simp <;> ring95 rw [this, ih]96 field_simp <;> ring9798/-- Boundary constant (Theorem 2(v)), deterministic part: if the pilot error is99 e0 = β θ c^r N^(-1/2^(q+1)) then √N · (−e0^(2^q)/θ^(2^q−1)) = −θ (β c^r)^(2^q). -/100theorem boundary_constant (N θ β c : ℝ) (r q : ℕ) (hN : 0 < N) (hθ : θ ≠ 0) :101 let e0 := β * θ * c ^ r * N ^ (-(1 : ℝ) / 2 ^ (q + 1))102 N ^ ((1:ℝ) / 2) * (-(e0 ^ (2 ^ q)) / θ ^ (2 ^ q - 1)) = -θ * (β * c ^ r) ^ (2 ^ q) := by103 intro e0104 have hq : (2:ℕ) ^ q ≥ 1 := Nat.one_le_two_pow105 have hpow : (N ^ (-(1 : ℝ) / 2 ^ (q + 1))) ^ (2 ^ q) = N ^ (-(1:ℝ) / 2) := by106 rw [← Real.rpow_natCast, ← Real.rpow_mul hN.le]107 congr 1108 push_cast109 field_simp <;> ring110 have hθp : θ ^ (2 ^ q - 1) ≠ 0 := pow_ne_zero _ hθ111 have hroot : N ^ ((1:ℝ) / 2) * N ^ (-(1:ℝ) / 2) = 1 := by112 rw [← Real.rpow_add hN]; norm_num113 have hθs : θ ^ (2 ^ q) = θ * θ ^ (2 ^ q - 1) := by114 rw [← pow_succ']; congr 1; omega115 calc N ^ ((1:ℝ) / 2) * (-(e0 ^ (2 ^ q)) / θ ^ (2 ^ q - 1))116 = N ^ ((1:ℝ) / 2) * (-((β * θ * c ^ r) ^ (2 ^ q) * (N ^ (-(1 : ℝ) / 2 ^ (q + 1))) ^ (2 ^ q)) / θ ^ (2 ^ q - 1)) := by117 simp only [e0, mul_pow]118 _ = -(β * θ * c ^ r) ^ (2 ^ q) / θ ^ (2 ^ q - 1) * (N ^ ((1:ℝ) / 2) * N ^ (-(1:ℝ) / 2)) := by119 rw [hpow]; ring120 _ = -(β * θ * c ^ r) ^ (2 ^ q) / θ ^ (2 ^ q - 1) := by rw [hroot, mul_one]121 _ = -θ * (β * c ^ r) ^ (2 ^ q) := by122 rw [show (β * θ * c ^ r) ^ (2 ^ q) = θ ^ (2 ^ q) * (β * c ^ r) ^ (2 ^ q) by ring, hθs]123 field_simp <;> ring124125/-! ### S7 Newton basin -/126127/-- Normalized error recursion u_{j+1} = -u_j^2. -/128def basinIter : ℕ → ℝ → ℝ129 | 0, u => u130 | n + 1, u => -(basinIter n u) ^ 2131132theorem basinIter_abs (u : ℝ) : ∀ j : ℕ, |basinIter j u| = |u| ^ (2 ^ j) := by133 intro j134 induction j with135 | zero => simp [basinIter]136 | succ n ih =>137 rw [basinIter, abs_neg, abs_pow, ih, ← pow_mul, pow_succ]138139/-- Convergence of the Newton map in the exponential model when the start lies in (0, 2θ̂):140 |u_j| = |u|^(2^j) ≤ |u|^j → 0 for |u| < 1. -/141theorem basin_converges (u : ℝ) (hu : |u| < 1) :142 Filter.Tendsto (fun j => basinIter j u) Filter.atTop (nhds 0) := by143 rw [tendsto_zero_iff_abs_tendsto_zero]144 refine squeeze_zero (fun j => abs_nonneg _) (fun j => ?_)145 (tendsto_pow_atTop_nhds_zero_of_lt_one (abs_nonneg u) hu)146 show |basinIter j u| ≤ |u| ^ j147 rw [basinIter_abs]148 exact pow_le_pow_of_le_one (abs_nonneg u) hu.le (Nat.lt_two_pow_self).le149150/-- Divergence when the start lies outside [0, 2θ̂]: |u_j| = |u|^(2^j) ≥ |u|^j → ∞ for |u| > 1. -/151theorem basin_diverges (u : ℝ) (hu : 1 < |u|) :152 Filter.Tendsto (fun j => |basinIter j u|) Filter.atTop Filter.atTop := by153 refine Filter.tendsto_atTop_mono (fun j => ?_) (tendsto_pow_atTop_atTop_of_one_lt hu)154 rw [basinIter_abs]155 exact pow_le_pow_right₀ hu.le (Nat.lt_two_pow_self).le156157/-! ### S8 Exponent algebra of the frontier (proof of Theorem 1) -/158159/-- With K = N^α, α = 1 − 1/(2^(q+1) r) and m = N/K = N^(1−α), the scaled bias term √N·m^(−r2^q) equals 1:160 the frontier exponent is exactly where the deterministic bias term stops vanishing. -/161theorem frontier_exponent (N : ℝ) (r q : ℕ) (hN : 0 < N) (hr : 0 < r) :162 let α : ℝ := 1 - 1 / (2 ^ (q + 1) * r)163 N ^ ((1:ℝ) / 2) * (N ^ ((1:ℝ) - α)) ^ (-((r : ℝ) * 2 ^ q)) = 1 := by164 intro α165 have hr' : (r : ℝ) ≠ 0 := by exact_mod_cast hr.ne'166 rw [← Real.rpow_mul hN.le, ← Real.rpow_add hN]167 have : (1:ℝ) / 2 + ((1:ℝ) - α) * (-((r : ℝ) * 2 ^ q)) = 0 := by168 simp only [α]169 field_simp <;> ring170 rw [this, Real.rpow_zero]171172/-- Below the frontier, K = N^α' with α' < α, the scaled bias term is N^(−(α−α') r 2^q), a negative power of N. -/173theorem frontier_below (N : ℝ) (r q : ℕ) (α' : ℝ) (hN : 0 < N) (hr : 0 < r) :174 let α : ℝ := 1 - 1 / (2 ^ (q + 1) * r)175 N ^ ((1:ℝ) / 2) * (N ^ ((1:ℝ) - α')) ^ (-((r : ℝ) * 2 ^ q)) = N ^ (-((α - α') * (r * 2 ^ q))) := by176 intro α177 have hr' : (r : ℝ) ≠ 0 := by exact_mod_cast hr.ne'178 rw [← Real.rpow_mul hN.le, ← Real.rpow_add hN]179 congr 1180 simp only [α]181 field_simp <;> ring182183/-! ### S2 / S6 Gamma moments and the exact variance of J_2 -/184185/-- E[Y^{-1}] = n Γ(n−1)/Γ(n) = n/(n−1) for n > 1 (nY ~ Gamma(n,1)). -/186theorem gamma_inv_moment1 (n : ℝ) (hn : 1 < n) :187 n * Real.Gamma (n - 1) / Real.Gamma n = n / (n - 1) := by188 have h : Real.Gamma n = (n - 1) * Real.Gamma (n - 1) := by189 have := Real.Gamma_add_one (s := n - 1) (by linarith)190 rw [show n - 1 + 1 = n by ring] at this; exact this191 have hpos : Real.Gamma (n - 1) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'192 have h1 : n - 1 ≠ 0 := by linarith193 rw [h]; field_simp194195/-- E[Y^{-2}] = n² Γ(n−2)/Γ(n) = n²/((n−1)(n−2)) for n > 2. -/196theorem gamma_inv_moment2 (n : ℝ) (hn : 2 < n) :197 n ^ 2 * Real.Gamma (n - 2) / Real.Gamma n = n ^ 2 / ((n - 1) * (n - 2)) := by198 have h1 : Real.Gamma n = (n - 1) * Real.Gamma (n - 1) := by199 have := Real.Gamma_add_one (s := n - 1) (by linarith)200 rw [show n - 1 + 1 = n by ring] at this; exact this201 have h2 : Real.Gamma (n - 1) = (n - 2) * Real.Gamma (n - 2) := by202 have := Real.Gamma_add_one (s := n - 2) (by linarith)203 rw [show n - 2 + 1 = n - 1 by ring] at this; exact this204 have hpos : Real.Gamma (n - 2) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'205 have e1 : n - 1 ≠ 0 := by linarith206 have e2 : n - 2 ≠ 0 := by linarith207 rw [h1, h2]; field_simp208209/-- Theorem 2(i): with E Y = 1, E Y² = 1 + 1/n, E Y⁻¹ = n/(n−1), E Y⁻² = n²/((n−1)(n−2)):210 E R_n = θ0/(n−1), E R_n² = θ0²(3n+2)/(n(n−1)(n−2)), Var θ̂_n = θ0² n²/((n−1)²(n−2)). -/211theorem exponential_remainder_moments (θ n : ℝ) (hn : 2 < n) :212 θ * (n / (n - 1) - 2 + 1) = θ / (n - 1) ∧213 θ ^ 2 * (n ^ 2 / ((n - 1) * (n - 2)) - 4 * (n / (n - 1)) + 6 - 4 * 1 + (1 + 1 / n))214 = θ ^ 2 * (3 * n + 2) / (n * (n - 1) * (n - 2)) ∧215 θ ^ 2 * (n ^ 2 / ((n - 1) * (n - 2)) - (n / (n - 1)) ^ 2) = θ ^ 2 * n ^ 2 / ((n - 1) ^ 2 * (n - 2)) := by216 have h1 : n - 1 ≠ 0 := by linarith217 have h2 : n - 2 ≠ 0 := by linarith218 have h0 : n ≠ 0 := by linarith219 refine ⟨?_, ?_, ?_⟩ <;> (field_simp; ring)220221/-- Beta(h,h) inverse moments as Gamma ratios reduce to rational functions of h (h > 2):222 E W⁻¹ = (2h−1)/(h−1), E W⁻² = (2h−1)(2h−2)/((h−1)(h−2)), E W⁻¹(1−W)⁻¹ = (2h−1)(2h−2)/(h−1)². -/223theorem beta_inverse_moments (h : ℝ) (hh : 2 < h) :224 Real.Gamma (h - 1) * Real.Gamma h * Real.Gamma (2 * h) / (Real.Gamma h ^ 2 * Real.Gamma (2 * h - 1))225 = (2 * h - 1) / (h - 1) ∧226 Real.Gamma (h - 2) * Real.Gamma h * Real.Gamma (2 * h) / (Real.Gamma h ^ 2 * Real.Gamma (2 * h - 2))227 = (2 * h - 1) * (2 * h - 2) / ((h - 1) * (h - 2)) ∧228 Real.Gamma (h - 1) * Real.Gamma (h - 1) * Real.Gamma (2 * h) / (Real.Gamma h ^ 2 * Real.Gamma (2 * h - 2))229 = (2 * h - 1) * (2 * h - 2) / (h - 1) ^ 2 := by230 have gh : Real.Gamma h = (h - 1) * Real.Gamma (h - 1) := by231 have := Real.Gamma_add_one (s := h - 1) (by linarith)232 rw [show h - 1 + 1 = h by ring] at this; exact this233 have gh1 : Real.Gamma (h - 1) = (h - 2) * Real.Gamma (h - 2) := by234 have := Real.Gamma_add_one (s := h - 2) (by linarith)235 rw [show h - 2 + 1 = h - 1 by ring] at this; exact this236 have g2h : Real.Gamma (2 * h) = (2 * h - 1) * Real.Gamma (2 * h - 1) := by237 have := Real.Gamma_add_one (s := 2 * h - 1) (by linarith)238 rw [show 2 * h - 1 + 1 = 2 * h by ring] at this; exact this239 have g2h1 : Real.Gamma (2 * h - 1) = (2 * h - 2) * Real.Gamma (2 * h - 2) := by240 have := Real.Gamma_add_one (s := 2 * h - 2) (by linarith)241 rw [show 2 * h - 2 + 1 = 2 * h - 1 by ring] at this; exact this242 have p1 : Real.Gamma (h - 2) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'243 have p2 : Real.Gamma (2 * h - 2) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'244 have p3 : Real.Gamma (h - 1) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'245 have p4 : Real.Gamma (2 * h - 1) ≠ 0 := (Real.Gamma_pos_of_pos (by linarith)).ne'246 have p1' : Real.Gamma (-1 + h) ≠ 0 := by rw [show (-1:ℝ) + h = h - 1 by ring]; exact p3247 have p2' : Real.Gamma (-1 + h * 2) ≠ 0 := by rw [show (-1:ℝ) + h * 2 = 2 * h - 1 by ring]; exact p4248 have p3' : Real.Gamma (-2 + h) ≠ 0 := by rw [show (-2:ℝ) + h = h - 2 by ring]; exact p1249 have p4' : Real.Gamma (-2 + h * 2) ≠ 0 := by rw [show (-2:ℝ) + h * 2 = 2 * h - 2 by ring]; exact p2250 have f1 : (-1:ℝ) + h ≠ 0 := by linarith251 have f2 : (-2:ℝ) + h ≠ 0 := by linarith252 have f3 : (-1:ℝ) + h * 2 ≠ 0 := by linarith253 have f4 : (-2:ℝ) + h * 2 ≠ 0 := by linarith254 have q1 : h * 2 - 1 ≠ 0 := by linarith255 have q2 : h * 2 - 2 ≠ 0 := by linarith256 have q3 : Real.Gamma (h * 2 - 1) ≠ 0 := by rw [show h * 2 - 1 = 2 * h - 1 by ring]; exact p4257 have q4 : Real.Gamma (h * 2 - 2) ≠ 0 := by rw [show h * 2 - 2 = 2 * h - 2 by ring]; exact p2258 have q5 : Real.Gamma (2 * (h - 1)) ≠ 0 := by rw [show 2 * (h - 1) = 2 * h - 2 by ring]; exact p2259 have e1 : h - 1 ≠ 0 := by linarith260 have e2 : h - 2 ≠ 0 := by linarith261 have e3 : 2 * h - 1 ≠ 0 := by linarith262 have e4 : 2 * h - 2 ≠ 0 := by linarith263 refine ⟨?_, ?_, ?_⟩264 · rw [gh, g2h]; field_simp265 · rw [gh, gh1, g2h, g2h1]; field_simp266 · rw [gh, gh1, g2h, g2h1]; field_simp267268/-- Theorem 2(iii): exact variance of J_2 (θ0 = 1, h = m/2). With E S⁻¹ = 1/(m−1), E S⁻² = 1/((m−1)(m−2)) and the269 Beta moments above, Var J_2 = E S⁻² E G² − (E S⁻¹ E G)² equals m²(m²−6m+17)/((m−4)(m−2)²(m−1)²),270 where G(W) = 2m − (h/2)(W⁻¹ + (1−W)⁻¹). -/271theorem var_J2_exact (m : ℝ) (hm : 4 < m) :272 (1 / ((m - 1) * (m - 2))) *273 (4 * m ^ 2 - 2 * (2 * m) * (m / 2 / 2) * (2 * ((2 * (m / 2) - 1) / (m / 2 - 1)))274 + ((m / 2) ^ 2 / 4) * (2 * ((2 * (m / 2) - 1) * (2 * (m / 2) - 2) / ((m / 2 - 1) * (m / 2 - 2)))275 + 2 * ((2 * (m / 2) - 1) * (2 * (m / 2) - 2) / (m / 2 - 1) ^ 2)))276 - ((1 / (m - 1)) * (2 * m - (m / 2 / 2) * (2 * ((2 * (m / 2) - 1) / (m / 2 - 1))))) ^ 2277 = m ^ 2 * (m ^ 2 - 6 * m + 17) / ((m - 4) * (m - 2) ^ 2 * (m - 1) ^ 2) := by278 have e1 : m - 1 ≠ 0 := by linarith279 have e2 : m - 2 ≠ 0 := by linarith280 have e4 : m - 4 ≠ 0 := by linarith281 have e5 : m / 2 - 1 ≠ 0 := by intro h; linarith282 have e6 : m / 2 - 2 ≠ 0 := by intro h; linarith283 have f1 : (-1:ℝ) + m ≠ 0 := by linarith284 have f2 : (-2:ℝ) + m ≠ 0 := by linarith285 have f4 : (-4:ℝ) + m ≠ 0 := by linarith286 have g4 : m - 2 ^ 2 ≠ 0 := by norm_num; linarith287 field_simp288 ring289290/-- Mean of J_2 from the same decomposition: E S⁻¹ · E G = 1 − 2/((m−1)(m−2)). -/291theorem mean_J2_exact (m : ℝ) (hm : 4 < m) :292 (1 / (m - 1)) * (2 * m - (m / 2 / 2) * (2 * ((2 * (m / 2) - 1) / (m / 2 - 1))))293 = 1 - 2 / ((m - 1) * (m - 2)) := by294 have e1 : m - 1 ≠ 0 := by linarith295 have e2 : m - 2 ≠ 0 := by linarith296 have e5 : m / 2 - 1 ≠ 0 := by intro h; linarith297 have f1 : (-1:ℝ) + m ≠ 0 := by linarith298 have f2 : (-2:ℝ) + m ≠ 0 := by linarith299 field_simp300 ring301302end WP28303