Mathlib.Probability.Certified­Finite­Markov­Expectation.Certified­Normal­CDFEnclosure

theorem reflectInterval_sound {I : RatInterval} {x : ℝ} (hx : I.Contains x) : (reflectInterval I).Contains (1 - x) := by rcases hx with ⟨hlo, hhi⟩ simp only [reflectInterval, RatInterval.sub, RatInterval.add, RatInterval.neg, RatInterval.point, RatInterval.Contains, Rat.cast_add, Rat.cast_neg, Rat.cast_one] exact ⟨sub_le_sub_left hhi 1, sub_le_sub_left hlo 1⟩

theorem reflectInterval_stdNormalCDF_sound {I : RatInterval} {x : ℝ} (hx : I.Contains (stdNormalCDF x)) : (reflectInterval I).Contains (stdNormalCDF (-x)) := by rw [stdNormalCDF_neg] exact reflectInterval_sound hx

theorem expTaylorInterval_sound (z : ℚ) (hz : |z| ≤ 1) (degree : ℕ) : (expTaylorInterval z degree).Contains (exp (z : ℝ)) := by -- Follow the reduced-input part of Causalean's `expRaw_sound`: apply -- `Real.exp_bound`, then identify its partial sum and remainder with -- `Transcendental.expPartial` and `expRemainder` by exact casts. have hzR : |(z : ℝ)| ≤ 1 := by exact_mod_cast hz have hb := Real.exp_bound hzR (n := degree + 1) (Nat.succ_pos degree) have hp : ((Transcendental.expPartial z degree : ℚ) : ℝ) = ∑ k ∈ range (degree + 1), (z : ℝ) ^ k / (k.factorial : ℝ) := by simp [Transcendental.expPartial] have hr : ((Transcendental.expRemainder z degree : ℚ) : ℝ) = |(z : ℝ)| ^ (degree + 1) * ((degree + 2 : ℕ) : ℝ) / (((degree + 1).factorial : ℝ) * (degree + 1 : ℕ)) := by simp [Transcendental.expRemainder] rw [← hp] at hb have hb' : |exp (z : ℝ) - (Transcendental.expPartial z degree : ℝ)| ≤ (Transcendental.expRemainder z degree : ℝ) := by rw [hr] convert hb using 1 <;> norm_num [Nat.cast_succ] <;> ring rw [abs_le] at hb' simp only [expTaylorInterval, Transcendental.expReducedRaw, RatInterval.Contains, Rat.cast_sub, Rat.cast_add] constructor <;> linarith [hb'.1, hb'.2]

theorem expTaylorInterval_width (z : ℚ) (degree : ℕ) : (expTaylorInterval z degree).width = 2 * Transcendental.expRemainder z degree := by -- Unfold the interval endpoints and its width, then use ring arithmetic. simp only [expTaylorInterval, Transcendental.expReducedRaw, RatInterval.width] ring

theorem two_pow_le_factorial_succ (n : ℕ) : 2 ^ n ≤ (n + 1).factorial := by -- Induct on `n`; rewrite both recurrences and use `2 ≤ n+2`. induction n with | zero => norm_num | succ n ih => rw [pow_succ, Nat.factorial_succ] calc 2 ^ n * 2 ≤ (n + 1).factorial * 2 := Nat.mul_le_mul_right 2 ih _ ≤ (n + 1).factorial * (n + 1 + 1) := by gcongr <;> omega _ = (n + 1 + 1) * (n + 1).factorial := Nat.mul_comm _ _

theorem expTaylorInterval_width_le_geometric (z : ℚ) (hz : |z| ≤ 1) (degree : ℕ) (hdegree : 1 ≤ degree) : (expTaylorInterval z degree).width ≤ 3 / (2 : ℚ) ^ degree := by -- Rewrite the width with `expTaylorInterval_width`. Bound `|z|^(n+1)` by -- one, `(n+2)/(n+1)` by `3/2`, and use -- `two_pow_le_factorial_succ`; clear only manifestly positive denominators. have habsPow : |z| ^ (degree + 1) ≤ (1 : ℚ) := pow_le_one₀ (abs_nonneg z) hz have hratioNat : 2 * (degree + 2) ≤ 3 * (degree + 1) := by omega have hratio : (2 : ℚ) * |z| ^ (degree + 1) * (degree + 2 : ℕ) ≤ 3 * (degree + 1 : ℕ) := by calc (2 : ℚ) * |z| ^ (degree + 1) * (degree + 2 : ℕ) ≤ 2 * 1 * (degree + 2 : ℕ) := by gcongr _ ≤ 3 * (degree + 1 : ℕ) := by exact_mod_cast hratioNat have hfactorialQ : (2 : ℚ) ^ degree ≤ ((degree + 1).factorial : ℚ) := by exact_mod_cast two_pow_le_factorial_succ degree have hcross : ((2 : ℚ) * |z| ^ (degree + 1) * (degree + 2 : ℕ)) * 2 ^ degree ≤ 3 * (((degree + 1).factorial : ℚ) * (degree + 1 : ℕ)) := by calc ((2 : ℚ) * |z| ^ (degree + 1) * (degree + 2 : ℕ)) * 2 ^ degree ≤ (3 * (degree + 1 : ℕ)) * 2 ^ degree := by gcongr _ ≤ (3 * (degree + 1 : ℕ)) * ((degree + 1).factorial : ℚ) := by gcongr _ = 3 * (((degree + 1).factorial : ℚ) * (degree + 1 : ℕ)) := by ring rw [expTaylorInterval_width, Transcendental.expRemainder] rw [show (2 : ℚ) * (|z| ^ (degree + 1) * (degree + 2 : ℕ) / (((degree + 1).factorial : ℚ) * (degree + 1 : ℕ))) = ((2 : ℚ) * |z| ^ (degree + 1) * (degree + 2 : ℕ)) / (((degree + 1).factorial : ℚ) * (degree + 1 : ℕ)) by ring] apply (div_le_div_iff₀ (by positivity) (by positivity)).2 exact hcross

theorem normalizationInterval_sound (s : NormalizationSchedule) (hcheck : normalizationScheduleCheck s = true) : (normalizationInterval s).Contains (1 / sqrt (2 * pi)) := by -- Unfold the executable guard twice. The successful `decide` supplies the -- `AwayFromZero` proof (proof irrelevance identifies it with the proof used -- by `normalizationInterval`). Then reuse `piInterval_sound`, -- `sqrtInterval_sound`, and `RatInterval.inv_sound`, as in Causalean's -- `normalDensityScaleInterval_sound`. unfold normalizationScheduleCheck at hcheck split at hcheck next htwoPi => dsimp only at hcheck have hroot := of_decide_eq_true hcheck rw [normalizationInterval, dif_pos htwoPi, dif_pos hroot] rw [one_div] apply RatInterval.inv_sound hroot apply RatInterval.sqrtInterval_sound htwoPi simpa [twoPiInterval] using RatInterval.mul_sound (RatInterval.point_sound 2) (Complex.Transcendental.piInterval_sound s.fuel) next => contradiction

theorem fineNormalCDFCheck_sound {q : ℚ} {method : MethodCertificate q} {reported : RatInterval} (hcheck : fineNormalCDFCheck q method reported = true) : reported.Contains (stdNormalCDF (q : ℝ)) ∧ |q| ≤ supportedEndpointBound ∧ reported.width ≤ targetWidth := by simp only [fineNormalCDFCheck, Bool.and_eq_true] at hcheck have hc := hcheck have hlimits := of_decide_eq_true hc.2 exact ⟨normalCDFCheck_sound hc.1, hlimits⟩

theorem FineNormalCDFCertificate.sound_range_and_width {q : ℚ} (c : FineNormalCDFCertificate q) : c.enclosure.Contains (stdNormalCDF (q : ℝ)) ∧ |q| ≤ supportedEndpointBound ∧ c.enclosure.width ≤ targetWidth := fineNormalCDFCheck_sound c.checked

theorem normalCDFDifference_sound {a b : ℚ} (ca : NormalCDFCertificate a) (cb : NormalCDFCertificate b) : (cb.enclosure.sub ca.enclosure).Contains (stdNormalCDF (b : ℝ) - stdNormalCDF (a : ℝ)) := by exact RatInterval.sub_sound cb.sound ca.sound

theorem fineNormalCDFDifference_sound {a b : ℚ} (ca : FineNormalCDFCertificate a) (cb : FineNormalCDFCertificate b) : (cb.enclosure.sub ca.enclosure).Contains (stdNormalCDF (b : ℝ) - stdNormalCDF (a : ℝ)) := by exact normalCDFDifference_sound ca.toNormalCDFCertificate cb.toNormalCDFCertificate

theorem gaussianIntegrand_hasSum (t : ℝ) : HasSum (fun k => gaussianIntegrandTerm k t) (exp (-(t ^ 2) / 2)) := by rw [Real.exp_eq_exp_ℝ] have hfun : (fun k => gaussianIntegrandTerm k t) = (fun k => ((k.factorial : ℝ)⁻¹ • (-(t ^ 2) / 2) ^ k)) := by funext k simp [gaussianIntegrandTerm, div_pow] ring rw [hfun] exact NormedSpace.exp_series_hasSum_exp' (𝕂 := ℝ) (x := (-(t ^ 2) / 2 : ℝ))

theorem gaussianIntegrandTerm_integral (x : ℝ) (hx : 0 ≤ x) (k : ℕ) : (∫ t in (0 : ℝ)..x, gaussianIntegrandTerm k t) = (-1 : ℝ) ^ k * gaussianIntegralMagnitude x k := by unfold gaussianIntegrandTerm gaussianIntegralMagnitude simp_rw [div_eq_mul_inv] rw [show (fun t : ℝ => (-1) ^ k * t ^ (2 * k) * (2 ^ k * (k.factorial : ℝ))⁻¹) = (fun t : ℝ => ((-1) ^ k * (2 ^ k * (k.factorial : ℝ))⁻¹) * t ^ (2 * k)) by funext t ring] rw [intervalIntegral.integral_const_mul, integral_pow] simp ring

theorem gaussianIntegrandTerm_norm_le (x : ℝ) (hx : 0 ≤ x) (k : ℕ) {t : ℝ} (ht : t ∈ Icc (0 : ℝ) x) : ‖gaussianIntegrandTerm k t‖ ≤ x ^ (2 * k) / ((2 : ℝ) ^ k * (k.factorial : ℝ)) := by -- Expand the norm and use `0 ≤ t ≤ x` to compare the even powers. have habs : |t| ≤ x := by simpa [abs_of_nonneg ht.1] using ht.2 have hp : |t| ^ (2 * k) ≤ x ^ (2 * k) := pow_le_pow_left₀ (abs_nonneg t) habs _ have hden : 0 ≤ (2 : ℝ) ^ k * (k.factorial : ℝ) := by positivity rw [gaussianIntegrandTerm, Real.norm_eq_abs, abs_div, abs_mul, abs_pow, abs_pow, abs_of_nonneg hden, abs_neg, abs_one, one_pow, one_mul] exact div_le_div_of_nonneg_right hp hden

theorem gaussianIntegrandMajorant_summable (x : ℝ) (hx : 0 ≤ x) : Summable (fun k => x ^ (2 * k) / ((2 : ℝ) ^ k * (k.factorial : ℝ))) := by -- Identify this nonnegative series with the absolute exponential series at -- `x²/2`, using `Real.summable_pow_div_factorial`. refine (Real.summable_pow_div_factorial (x ^ 2 / 2)).congr ?_ intro k rw [div_pow, ← pow_mul] ring

theorem gaussianIntegralMagnitude_nonneg (x : ℝ) (hx : 0 ≤ x) (k : ℕ) : 0 ≤ gaussianIntegralMagnitude x k := by -- Unfold; every factor in the numerator and denominator is nonnegative and -- the natural-number factors in the denominator are strictly positive. unfold gaussianIntegralMagnitude positivity

theorem gaussianIntegralMagnitude_antitone_from (x : ℝ) (hx : 0 ≤ x) (n : ℕ) (hmono : x ^ 2 ≤ 2 * (n + 2 : ℕ)) : ∀ k, n + 1 ≤ k → gaussianIntegralMagnitude x (k + 1) ≤ gaussianIntegralMagnitude x k := by intro k hk have hk' : n + 2 ≤ k + 1 := by omega have hbound : x ^ 2 ≤ 2 * ((k + 1 : ℕ) : ℝ) := by apply hmono.trans exact_mod_cast Nat.mul_le_mul_left 2 hk' have h₁ : x ^ 2 / (2 * (k + 1 : ℕ)) ≤ (1 : ℝ) := (div_le_one (by positivity)).2 hbound have h₁nonneg : 0 ≤ x ^ 2 / (2 * (k + 1 : ℕ) : ℝ) := by positivity have h₂ : ((2 * k + 1 : ℕ) : ℝ) / (2 * k + 3 : ℕ) ≤ (1 : ℝ) := by apply (div_le_one (by positivity)).2 norm_num have h₂nonneg : 0 ≤ ((2 * k + 1 : ℕ) : ℝ) / (2 * k + 3 : ℕ) := by positivity have hid : gaussianIntegralMagnitude x (k + 1) = gaussianIntegralMagnitude x k * (x ^ 2 / (2 * (k + 1 : ℕ))) * (((2 * k + 1 : ℕ) : ℝ) / (2 * k + 3 : ℕ)) := by unfold gaussianIntegralMagnitude rw [Nat.factorial_succ] push_cast field_simp ring rw [hid] have hm := gaussianIntegralMagnitude_nonneg x hx k calc gaussianIntegralMagnitude x k * (x ^ 2 / (2 * ↑(k + 1))) * (↑(2 * k + 1) / ↑(2 * k + 3)) ≤ gaussianIntegralMagnitude x k * 1 * (↑(2 * k + 1) / ↑(2 * k + 3)) := mul_le_mul_of_nonneg_right (mul_le_mul_of_nonneg_left h₁ hm) h₂nonneg _ ≤ gaussianIntegralMagnitude x k * 1 := by simpa only [mul_one] using mul_le_mul_of_nonneg_left h₂ (mul_nonneg hm zero_le_one) _ = gaussianIntegralMagnitude x k := by ring

theorem gaussianIntegral_hasSum (x : ℝ) (hx : 0 ≤ x) : HasSum (fun k => (-1 : ℝ) ^ k * gaussianIntegralMagnitude x k) (∫ t in (0 : ℝ)..x, exp (-(t ^ 2) / 2)) := by -- Rewrite `exp (-(t^2)/2)` by `Real.hasSum_exp`, justify termwise -- integration on `[0,x]` using summable sup norms, then integrate monomials. -- `intervalIntegral.hasSum_intervalIntegral_of_summable_norm` is available. let f : ℕ → C(ℝ, ℝ) := fun k => ⟨gaussianIntegrandTerm k, by unfold gaussianIntegrandTerm fun_prop⟩ have hnorm : Summable (fun k => ‖(f k).restrict (⟨Set.uIcc (0 : ℝ) x, isCompact_uIcc⟩ : TopologicalSpace.Compacts ℝ)‖) := by apply Summable.of_nonneg_of_le (fun _ => norm_nonneg _) (fun k => ?_) (gaussianIntegrandMajorant_summable x hx) apply (ContinuousMap.norm_le _ (by positivity)).2 intro t change ‖gaussianIntegrandTerm k t‖ ≤ _ simpa only [ContinuousMap.restrict_apply] using gaussianIntegrandTerm_norm_le x hx k (by simpa [Set.uIcc_of_le hx] using t.property) have hi := intervalIntegral.hasSum_intervalIntegral_of_summable_norm hnorm have hterms : (fun k => ∫ t in (0 : ℝ)..x, f k t) = (fun k => (-1 : ℝ) ^ k * gaussianIntegralMagnitude x k) := by funext k change (∫ t in (0 : ℝ)..x, gaussianIntegrandTerm k t) = _ exact gaussianIntegrandTerm_integral x hx k rw [hterms] at hi rw [show (∫ t in (0 : ℝ)..x, exp (-(t ^ 2) / 2)) = ∫ t in (0 : ℝ)..x, ∑' k, f k t by apply intervalIntegral.integral_congr intro t _ change exp (-(t ^ 2) / 2) = ∑' k, gaussianIntegrandTerm k t exact (gaussianIntegrand_hasSum t).tsum_eq.symm] exact hi

theorem gaussianIntegralRemainder_geometric_step (x : ℚ) (hx0 : 0 ≤ x) (hx8 : x ≤ centralCutoff) (n : ℕ) (hn : 64 ≤ n) : gaussianIntegralRemainder x (n + 1) ≤ gaussianIntegralRemainder x n / 2 := by -- Divide consecutive explicit remainders. With `x ≤ 8` and `n ≥ 64`, the -- ratio is bounded by `32/(n+2) ≤ 1/2`; all denominators are positive. have hx8' : x ≤ 8 := by simpa [centralCutoff] using hx8 have hx2 : x ^ 2 ≤ 64 := by nlinarith have hn2 : (66 : ℚ) ≤ (n + 2 : ℕ) := by exact_mod_cast (show 66 ≤ n + 2 by omega) have hratio : x ^ 2 / (2 * (n + 2 : ℕ) : ℚ) ≤ 1 / 2 := by apply (div_le_iff₀ (by positivity)).2 nlinarith have hfactor : ((2 * n + 3 : ℕ) : ℚ) / (2 * n + 5 : ℕ) ≤ 1 := by apply (div_le_one (by positivity)).2 norm_num have hfactor0 : 0 ≤ ((2 * n + 3 : ℕ) : ℚ) / (2 * n + 5 : ℕ) := by positivity have hid : gaussianIntegralRemainder x (n + 1) = gaussianIntegralRemainder x n * (x ^ 2 / (2 * (n + 2 : ℕ))) * (((2 * n + 3 : ℕ) : ℚ) / (2 * n + 5 : ℕ)) := by unfold gaussianIntegralRemainder rw [abs_of_nonneg hx0, Nat.factorial_succ] push_cast field_simp ring rw [hid] have hr0 : 0 ≤ gaussianIntegralRemainder x n := by unfold gaussianIntegralRemainder positivity calc gaussianIntegralRemainder x n * (x ^ 2 / (2 * ↑(n + 2))) * (↑(2 * n + 3) / ↑(2 * n + 5)) ≤ gaussianIntegralRemainder x n * (1 / 2) * (↑(2 * n + 3) / ↑(2 * n + 5)) := mul_le_mul_of_nonneg_right (mul_le_mul_of_nonneg_left hratio hr0) hfactor0 _ ≤ gaussianIntegralRemainder x n * (1 / 2) * 1 := mul_le_mul_of_nonneg_left hfactor (mul_nonneg hr0 (by norm_num)) _ = gaussianIntegralRemainder x n / 2 := by ring

theorem gaussianIntegralRemainder_geometric (x : ℚ) (hx0 : 0 ≤ x) (hx8 : x ≤ centralCutoff) (j : ℕ) : gaussianIntegralRemainder x (64 + j) ≤ gaussianIntegralRemainder x 64 / (2 : ℚ) ^ j := by -- Induct on `j`, invoke the one-step bound at `64+j`, and clear the positive -- power-of-two denominator. induction j with | zero => simp | succ j ih => have hstep := gaussianIntegralRemainder_geometric_step x hx0 hx8 (64 + j) (by omega) rw [show 64 + (j + 1) = (64 + j) + 1 by omega] calc gaussianIntegralRemainder x ((64 + j) + 1) ≤ gaussianIntegralRemainder x (64 + j) / 2 := hstep _ ≤ (gaussianIntegralRemainder x 64 / (2 : ℚ) ^ j) / 2 := div_le_div_of_nonneg_right ih (by norm_num) _ = gaussianIntegralRemainder x 64 / (2 : ℚ) ^ (j + 1) := by rw [pow_succ] ring

theorem hasDerivAt_stdNormalPDF (x : ℝ) : HasDerivAt stdNormalPDF (-x * stdNormalPDF x) x := by have hpdf : stdNormalPDF = fun y : ℝ => (sqrt (2 * pi))⁻¹ * exp (-y ^ 2 / 2) := by funext y unfold stdNormalPDF gaussianPDFReal congr 2 · norm_num · norm_num rw [hpdf] have hpow : HasDerivAt (fun y : ℝ => -y ^ 2 / 2) (-x) x := by have h := ((hasDerivAt_pow 2 x).div_const 2).fun_neg simpa [neg_div, pow_one] using h.congr_deriv (by ring) have hcomp : HasDerivAt (fun y : ℝ => exp (-y ^ 2 / 2)) (exp (-x ^ 2 / 2) * (-x)) x := (Real.hasDerivAt_exp _).comp x hpow simpa only [mul_assoc, mul_comm, mul_left_comm, neg_mul, mul_neg] using hcomp.const_mul (sqrt (2 * pi))⁻¹

theorem millsRatio_upper (x : ℝ) (hx : 0 < x) : 1 - stdNormalCDF x ≤ stdNormalPDF x / x := by have hcompare : (∫ t in Ioi x, stdNormalPDF t) ≤ ∫ t in Ioi x, (1 / x) * (t * stdNormalPDF t) := by have hpdf : Integrable (fun t : ℝ => stdNormalPDF t) := by change Integrable (gaussianPDFReal 0 1) exact ProbabilityTheory.integrable_gaussianPDFReal 0 1 have hweighted : Integrable (fun t : ℝ => (1 / x) * (t * stdNormalPDF t)) := (integrable_id_mul_stdNormalPDF).const_mul (1 / x) refine MeasureTheory.setIntegral_mono_on_ae hpdf.integrableOn hweighted.integrableOn measurableSet_Ioi ?_ filter_upwards with t intro ht have hxt : x ≤ t := le_of_lt ht have hratio : 1 ≤ t / x := (le_div_iff₀ hx).2 (by simpa using hxt) calc stdNormalPDF t = 1 * stdNormalPDF t := by ring _ ≤ (t / x) * stdNormalPDF t := mul_le_mul_of_nonneg_right hratio (stdNormalPDF_pos t).le _ = (1 / x) * (t * stdNormalPDF t) := by ring calc 1 - stdNormalCDF x = ∫ t in Ioi x, stdNormalPDF t := one_sub_stdNormalCDF_eq_integral_Ioi x _ ≤ ∫ t in Ioi x, (1 / x) * (t * stdNormalPDF t) := hcompare _ = (1 / x) * ∫ t in Ioi x, t * stdNormalPDF t := by rw [MeasureTheory.integral_const_mul] _ = stdNormalPDF x / x := by rw [integral_Ioi_id_mul_stdNormalPDF] ring

theorem millsRatio_lower (x : ℝ) (hx : 0 < x) : x / (x ^ 2 + 1) * stdNormalPDF x ≤ 1 - stdNormalCDF x := by have hpdf_int : Integrable stdNormalPDF := by change Integrable (gaussianPDFReal 0 1) exact ProbabilityTheory.integrable_gaussianPDFReal 0 1 have hpdf_cont : Continuous stdNormalPDF := by exact continuous_iff_continuousAt.2 fun t => (hasDerivAt_stdNormalPDF t).continuousAt have hrem_int : IntegrableOn (fun t : ℝ => stdNormalPDF t / t ^ 2) (Ioi x) := by have hmajor : IntegrableOn (fun t : ℝ => (1 / x ^ 2) * stdNormalPDF t) (Ioi x) := (hpdf_int.const_mul (1 / x ^ 2)).integrableOn refine hmajor.mono' ((hpdf_cont.continuousOn.div (continuous_id.pow 2).continuousOn (fun t ht => pow_ne_zero 2 (ne_of_gt (hx.trans ht)))).aestronglyMeasurable measurableSet_Ioi) ?_ filter_upwards [MeasureTheory.ae_restrict_mem measurableSet_Ioi] with t ht have htpos : 0 < t := hx.trans ht have hxt : x < t := ht have hsq : x ^ 2 ≤ t ^ 2 := by nlinarith have hinv : 1 / t ^ 2 ≤ 1 / x ^ 2 := (div_le_div_iff₀ (pow_pos htpos 2) (pow_pos hx 2)).2 (by simpa using hsq) rw [Real.norm_eq_abs, abs_of_nonneg (div_nonneg (stdNormalPDF_pos t).le (sq_nonneg t))] calc stdNormalPDF t / t ^ 2 = (1 / t ^ 2) * stdNormalPDF t := by ring _ ≤ (1 / x ^ 2) * stdNormalPDF t := mul_le_mul_of_nonneg_right hinv (stdNormalPDF_pos t).le have hpdf_tendsto : Tendsto stdNormalPDF atTop (nhds 0) := by have hpdf : stdNormalPDF = fun y : ℝ => (sqrt (2 * pi))⁻¹ * exp (-y ^ 2 / 2) := by funext y simp [stdNormalPDF, gaussianPDFReal] rw [hpdf] have hsq : Tendsto (fun t : ℝ => t ^ 2) atTop atTop := Filter.tendsto_pow_atTop (α := ℝ) (n := 2) (by norm_num) have hneg : Tendsto (fun t : ℝ => -t ^ 2 / 2) atTop atBot := by apply Filter.Tendsto.atBot_div_const (by norm_num) exact Filter.tendsto_neg_atBot_iff.mpr hsq simpa using (tendsto_const_nhds.mul (Real.tendsto_exp_atBot.comp hneg)) have hquot_tendsto : Tendsto (fun t : ℝ => -stdNormalPDF t / t) atTop (nhds 0) := by simpa [div_eq_mul_inv] using (hpdf_tendsto.mul tendsto_inv_atTop_zero).neg have hderiv : ∀ t ∈ Ici x, HasDerivAt (fun y : ℝ => -stdNormalPDF y / y) (stdNormalPDF t + stdNormalPDF t / t ^ 2) t := by intro t ht have ht0 : t ≠ 0 := ne_of_gt (hx.trans_le ht) change HasDerivAt ((fun y : ℝ => -stdNormalPDF y) / id) (stdNormalPDF t + stdNormalPDF t / t ^ 2) t have h := (hasDerivAt_stdNormalPDF t).fun_neg.div (hasDerivAt_id t) ht0 exact h.congr_deriv (by simp only [id_eq] field_simp [ht0] ring) have htail_int : IntegrableOn stdNormalPDF (Ioi x) := hpdf_int.integrableOn have hsum_eq : (∫ t in Ioi x, stdNormalPDF t + stdNormalPDF t / t ^ 2) = stdNormalPDF x / x := by have h := MeasureTheory.integral_Ioi_of_hasDerivAt_of_tendsto' (f := fun t : ℝ => -stdNormalPDF t / t) (f' := fun t : ℝ => stdNormalPDF t + stdNormalPDF t / t ^ 2) (a := x) (m := 0) hderiv (htail_int.add hrem_int) hquot_tendsto simpa [hx.ne', neg_div] using h have hdecomp : stdNormalPDF x / x = (∫ t in Ioi x, stdNormalPDF t) + ∫ t in Ioi x, stdNormalPDF t / t ^ 2 := by rw [← MeasureTheory.integral_add htail_int hrem_int] exact hsum_eq.symm have hrem_le : (∫ t in Ioi x, stdNormalPDF t / t ^ 2) ≤ (1 / x ^ 2) * ∫ t in Ioi x, stdNormalPDF t := by calc (∫ t in Ioi x, stdNormalPDF t / t ^ 2) ≤ ∫ t in Ioi x, (1 / x ^ 2) * stdNormalPDF t := by refine MeasureTheory.setIntegral_mono_of_nonneg ?_ ?_ ((hpdf_int.const_mul (1 / x ^ 2)).integrableOn) · intro t ht exact div_nonneg (stdNormalPDF_pos t).le (sq_nonneg t) · intro t ht have htpos : 0 < t := hx.trans ht have hxt : x < t := ht have hsq : x ^ 2 ≤ t ^ 2 := by nlinarith have hinv : 1 / t ^ 2 ≤ 1 / x ^ 2 := (div_le_div_iff₀ (pow_pos htpos 2) (pow_pos hx 2)).2 (by simpa using hsq) calc stdNormalPDF t / t ^ 2 = (1 / t ^ 2) * stdNormalPDF t := by ring _ ≤ (1 / x ^ 2) * stdNormalPDF t := mul_le_mul_of_nonneg_right hinv (stdNormalPDF_pos t).le _ = (1 / x ^ 2) * ∫ t in Ioi x, stdNormalPDF t := by rw [MeasureTheory.integral_const_mul] have hmain : stdNormalPDF x / x ≤ (∫ t in Ioi x, stdNormalPDF t) + (1 / x ^ 2) * ∫ t in Ioi x, stdNormalPDF t := by rw [hdecomp] exact add_le_add_right hrem_le _ have hscaled := mul_le_mul_of_nonneg_right hmain (sq_nonneg x) field_simp [hx.ne'] at hscaled rw [one_sub_stdNormalCDF_eq_integral_Ioi] rw [div_mul_eq_mul_div, mul_comm x] apply (div_le_iff₀ (by positivity : 0 < x ^ 2 + 1)).2 exact hscaled

theorem millsRatio_bounds (x : ℝ) (hx : 0 < x) : x / (x ^ 2 + 1) * stdNormalPDF x ≤ 1 - stdNormalCDF x ∧ 1 - stdNormalCDF x ≤ stdNormalPDF x / x := by exact ⟨millsRatio_lower x hx, millsRatio_upper x hx⟩

theorem tailDensityInterval_sound (q : ℚ) (c : TailCertificate q) (hnormalization : normalizationCheck c.normalization = true) (hexponential : expCheck (-|q| ^ 2 / 2) c.exponential = true) : (tailDensityInterval q c).Contains (stdNormalPDF (|q| : ℚ)) := by unfold tailDensityInterval have h := RatInterval.mul_sound (c.normalization.sound hnormalization) (c.exponential.sound hexponential) simpa [stdNormalPDF, gaussianPDFReal] using h

theorem millsTailInterval_sound (x : ℚ) (hx : 0 < x) (density : RatInterval) (hdensity : density.Contains (stdNormalPDF (x : ℝ))) : (millsTailInterval x density).Contains (1 - stdNormalCDF (x : ℝ)) := by have hxR : (0 : ℝ) < (x : ℝ) := by exact_mod_cast hx have hlower := RatInterval.mul_sound (RatInterval.point_sound (x / (x ^ 2 + 1))) hdensity have hupper := RatInterval.mul_sound (RatInterval.point_sound (1 / x)) hdensity have hmills := millsRatio_bounds (x : ℝ) hxR have hlowerLo := hlower.1 have hupperHi := hupper.2 constructor · simp only [millsTailInterval, RatInterval.hull, Rat.cast_min] refine (min_le_left _ _).trans (hlowerLo.trans ?_) simpa only [Rat.cast_div, Rat.cast_add, Rat.cast_pow, Rat.cast_one] using hmills.1 · simp only [millsTailInterval, RatInterval.hull, Rat.cast_max] refine hmills.2.trans ?_ calc stdNormalPDF (x : ℝ) / (x : ℝ) = ((1 / x : ℚ) : ℝ) * stdNormalPDF (x : ℝ) := by simp only [Rat.cast_div, Rat.cast_one] ring _ ≤ (((RatInterval.point (1 / x)).mul density).hi : ℚ) := hupperHi _ ≤ max (((RatInterval.point (x / (x ^ 2 + 1))).mul density).hi : ℝ) (((RatInterval.point (1 / x)).mul density).hi : ℝ) := le_max_right _ _