Stat.Bootstrap.Asymptotic­Linear

theorem bootstrapQuantile_eq_affine (est : (n : ℕ) -> (Fin n -> X) -> ℝ) (hest : forall n, Measurable (est n)) {n : ℕ} (hn : n != 0) {beta : ℝ} (hbeta0 : 0 < beta) (hbeta1 : beta < 1) (x : Fin n -> X) : bootstrapQuantile est n beta x = sqrt (n : ℝ) * (bootstrapLowerQuantile (est n) beta x - est n x) := by -- `sqrt n` is positive. Rewrite both lower quantiles as `sInf` sets and transport the CDF -- threshold set through the increasing affine map `y ↦ sqrt n * (y - est n x)`. have hn0 : n ≠ 0 := bne_iff_ne.mp hn have hsqrt : 0 < sqrt (n : ℝ) := Real.sqrt_pos.2 (by exact_mod_cast Nat.pos_of_ne_zero hn0) have hcenter := measurable_centeredEstimatorBootstrapStatistic est hest n x have hbeta1' : beta ≤ 1 := hbeta1.le have hcdf (t : ℝ) : bootstrapCDF (centeredEstimatorBootstrapStatistic est n x) x (sqrt (n : ℝ) * (t - est n x)) = bootstrapCDF (est n) x t := by rw [bootstrapCDF_eq_average_indicators _ hcenter x _, bootstrapCDF_eq_average_indicators _ (hest n) x t] congr 1 apply Finset.sum_congr rfl intro j hj have hiff : centeredEstimatorBootstrapStatistic est n x (x ∘ j) ≤ sqrt (n : ℝ) * (t - est n x) ↔ est n (x ∘ j) ≤ t := by unfold centeredEstimatorBootstrapStatistic rw [mul_le_mul_iff_of_pos_left hsqrt, sub_le_sub_iff_right] simp only [hiff] apply le_antisymm · change bootstrapLowerQuantile (centeredEstimatorBootstrapStatistic est n x) beta x ≤ _ rw [bootstrapLowerQuantile_le_iff _ hcenter x hn0 beta _ hbeta0 hbeta1', hcdf] exact (bootstrapLowerQuantile_le_iff (est n) (hest n) x hn0 beta _ hbeta0 hbeta1').mp le_rfl · let qstar := bootstrapQuantile est n beta x have hreach : beta ≤ bootstrapCDF (centeredEstimatorBootstrapStatistic est n x) x qstar := (bootstrapLowerQuantile_le_iff _ hcenter x hn0 beta qstar hbeta0 hbeta1').mp le_rfl have hcdf' : bootstrapCDF (centeredEstimatorBootstrapStatistic est n x) x qstar = bootstrapCDF (est n) x (qstar / sqrt (n : ℝ) + est n x) := by have hsqrt_ne : sqrt (n : ℝ) ≠ 0 := ne_of_gt hsqrt rw [← hcdf (qstar / sqrt (n : ℝ) + est n x)] congr 2 field_simp [hsqrt_ne] <;> ring have hq : bootstrapLowerQuantile (est n) beta x ≤ qstar / sqrt (n : ℝ) + est n x := (bootstrapLowerQuantile_le_iff (est n) (hest n) x hn0 beta _ hbeta0 hbeta1').mpr (by rwa [← hcdf']) change sqrt (n : ℝ) * (bootstrapLowerQuantile (est n) beta x - est n x) ≤ qstar calc sqrt (n : ℝ) * (bootstrapLowerQuantile (est n) beta x - est n x) ≤ sqrt (n : ℝ) * ((qstar / sqrt (n : ℝ) + est n x) - est n x) := mul_le_mul_of_nonneg_left (sub_le_sub_right hq _) hsqrt.le _ = qstar := by field_simp [ne_of_gt hsqrt] <;> ring

theorem measurable_bootstrapQuantile (est : (n : ℕ) -> (Fin n -> X) -> ℝ) (hest : forall n, Measurable (est n)) (n : ℕ) {beta : ℝ} (hbeta0 : 0 < beta) (hbeta1 : beta < 1) : Measurable (bootstrapQuantile est n beta) := by -- Split on `n = 0` (the statistic is constant) and otherwise use -- `bootstrapQuantile_eq_affine` plus `measurable_bootstrapLowerQuantile`. by_cases hn : n = 0 · subst n let x0 : Fin 0 → X := fun i ↦ Fin.elim0 i have heq : bootstrapQuantile est 0 beta = fun _ ↦ bootstrapQuantile est 0 beta x0 := by funext x congr 1 exact Subsingleton.elim _ _ rw [heq] exact measurable_const · rw [show bootstrapQuantile est n beta = fun x ↦ sqrt (n : ℝ) * (bootstrapLowerQuantile (est n) beta x - est n x) by funext x exact bootstrapQuantile_eq_affine est hest (bne_iff_ne.mpr hn) hbeta0 hbeta1 x] exact measurable_const.mul ((measurable_bootstrapLowerQuantile (est n) (hest n) hn beta).sub (hest n))

theorem percentileCI_coverage (h : BootstrapAsymLinear S est theta0 psi) {alpha : ℝ} (halpha0 : 0 < alpha) (halpha1 : alpha < 1) : Tendsto (fun n => mu.real {omega | theta0 ∈ percentileCI est n alpha (S.sampleVector n omega)}) atTop (nhds (1 - alpha)) := by -- Rewrite target membership as the rescaled sampling statistic lying between the negatives of -- the upper/lower bootstrap quantiles. The random endpoints converge by `quantile_tendsto`; -- apply random-endpoint portmanteau and the reflected Gaussian equal-tail mass lemma. letI : IsProbabilityMeasure mu := S.indep.isProbabilityMeasure have hvReal : 0 < h.asymptoticVariance := h.var_pos_finite.1 have hv : 0 < h.asymptoticVariance.toNNReal := by rwa [Real.toNNReal_pos] let G : ProbabilityMeasure ℝ := ⟨gaussianReal 0 h.asymptoticVariance.toNNReal, inferInstance⟩ have hG : h.asymptoticGaussian = G := by apply ProbabilityMeasure.toMeasure_injective rfl let T : ℕ → Omega → ℝ := fun n omega => sqrt (n : ℝ) * (est n (S.sampleVector n omega) - theta0) have hx (n : ℕ) : Measurable (fun omega => S.sampleVector n omega) := by apply measurable_pi_iff.mpr intro i exact S.meas i have hTMeas (n : ℕ) : AEMeasurable (T n) mu := by exact (measurable_const.mul (((h.meas.1 n).comp (hx n)).sub measurable_const)).aemeasurable have hT : Modes.TendstoInLaw (fun _ : ℕ => mu) T atTop (G : Measure ℝ) := by apply (Tendsto_dist_iff _ _ _ hTMeas).2 have hlaw : (fun n => (⟨mu.map (T n), Measure.isProbabilityMeasure_map (hTMeas n)⟩ : ProbabilityMeasure ℝ)) = samplingEstimatorLaw S est theta0 h.meas.1 := by funext n apply ProbabilityMeasure.toMeasure_injective rfl rw [hlaw, ← hG] exact sampling_tendsto_gaussian h have ha2_pos : 0 < alpha / 2 := by linarith have ha2_lt_one : alpha / 2 < 1 := by linarith have h1a2_pos : 0 < 1 - alpha / 2 := by linarith have h1a2_lt_one : 1 - alpha / 2 < 1 := by linarith let qlo : ℕ → Omega → ℝ := fun n omega => bootstrapQuantile est n (alpha / 2) (S.sampleVector n omega) let qhi : ℕ → Omega → ℝ := fun n omega => bootstrapQuantile est n (1 - alpha / 2) (S.sampleVector n omega) let alo := sqrt h.asymptoticVariance * probit (alpha / 2) let ahi := sqrt h.asymptoticVariance * probit (1 - alpha / 2) have hqloMeas (n : ℕ) : Measurable (qlo n) := (measurable_bootstrapQuantile est h.meas.1 n ha2_pos ha2_lt_one).comp (hx n) have hqhiMeas (n : ℕ) : Measurable (qhi n) := (measurable_bootstrapQuantile est h.meas.1 n h1a2_pos h1a2_lt_one).comp (hx n) have hqlo : Modes.TendstoInProbability (fun _ : ℕ => mu) qlo atTop (fun _ _ => alo) := quantile_tendsto h ha2_pos ha2_lt_one have hqhi : Modes.TendstoInProbability (fun _ : ℕ => mu) qhi atTop (fun _ _ => ahi) := quantile_tendsto h h1a2_pos h1a2_lt_one have hnegqhi : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => -qhi n omega) atTop (fun _ _ => -ahi) := by simpa only [neg_zero] using Modes.TendstoInProbability.comp_continuousAt (g := fun x : ℝ => -x) continuousAt_neg hqhi have hnegqlo : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => -qlo n omega) atTop (fun _ _ => -alo) := by simpa only [neg_zero] using Modes.TendstoInProbability.comp_continuousAt (g := fun x : ℝ => -x) continuousAt_neg hqlo have hprobit : probit (alpha / 2) ≤ probit (1 - alpha / 2) := by apply Causalean.Mathlib.stdNormalCDF_strictMono.le_iff_le.mp rw [stdNormalCDF_probit ha2_pos ha2_lt_one, stdNormalCDF_probit h1a2_pos h1a2_lt_one] linarith have halohi : alo ≤ ahi := by dsimp [alo, ahi] exact mul_le_mul_of_nonneg_left hprobit (Real.sqrt_nonneg _) have hboundary : (G : Measure ℝ) (frontier (Icc (-ahi) (-alo))) = 0 := by haveI : NullSingletonClass (G : Measure ℝ) := by dsimp [G] exact nullSingletonClass_gaussianReal hv.ne' rw [frontier_Icc (neg_le_neg halohi)] apply le_antisymm · calc (G : Measure ℝ) ({-ahi, -alo} : Set ℝ) ≤ (G : Measure ℝ) ({-ahi} : Set ℝ) + (G : Measure ℝ) ({-alo} : Set ℝ) := by rw [show ({-ahi, -alo} : Set ℝ) = {-ahi} ∪ {-alo} by ext x; simp [or_comm]] exact measure_union_le ({-ahi} : Set ℝ) ({-alo} : Set ℝ) _ = 0 := by simp only [NullSingletonClass.measure_singleton, add_zero] · exact zero_le have hcoverage := tendsto_dist_random_Icc_coverage hTMeas hT (fun n omega => -qhi n omega) (fun n omega => -qlo n omega) (-ahi) (-alo) (fun n => (hqhiMeas n).neg) (fun n => (hqloMeas n).neg) hnegqhi hnegqlo (neg_le_neg halohi) hboundary have hmass : (G : Measure ℝ).real (Icc (-ahi) (-alo)) = 1 - alpha := by dsimp [G, alo, ahi] simpa only [Real.coe_toNNReal _ hvReal.le] using gaussianReal_reflectedQuantile_Icc_toReal hv halpha0 halpha1 rw [hmass] at hcoverage apply hcoverage.congr' filter_upwards [eventually_ne_atTop 0] with n hn congr 1 ext omega have hsqrt : 0 < sqrt (n : ℝ) := Real.sqrt_pos.2 (by exact_mod_cast Nat.pos_of_ne_zero hn) change (-qhi n omega ≤ T n omega ∧ T n omega ≤ -qlo n omega) ↔ (est n (S.sampleVector n omega) + qlo n omega / sqrt (n : ℝ) ≤ theta0 ∧ theta0 ≤ est n (S.sampleVector n omega) + qhi n omega / sqrt (n : ℝ)) rw [iff_comm] dsimp [T] constructor · rintro ⟨hlow, hhigh⟩ constructor · have hz : (theta0 - est n (S.sampleVector n omega)) * sqrt (n : ℝ) ≤ qhi n omega := (le_div_iff₀ hsqrt).mp (by linarith) nlinarith · have hz : qlo n omega ≤ (theta0 - est n (S.sampleVector n omega)) * sqrt (n : ℝ) := (div_le_iff₀ hsqrt).mp (by linarith) nlinarith · rintro ⟨hlow, hhigh⟩ constructor · have hz : qlo n omega / sqrt (n : ℝ) ≤ theta0 - est n (S.sampleVector n omega) := (div_le_iff₀ hsqrt).mpr (by nlinarith) linarith · have hz : theta0 - est n (S.sampleVector n omega) ≤ qhi n omega / sqrt (n : ℝ) := (le_div_iff₀ hsqrt).mpr (by nlinarith) linarith

theorem basicCI_coverage (h : BootstrapAsymLinear S est theta0 psi) {alpha : ℝ} (halpha0 : 0 < alpha) (halpha1 : alpha < 1) : Tendsto (fun n => mu.real {omega | theta0 ∈ basicCI est n alpha (S.sampleVector n omega)}) atTop (nhds (1 - alpha)) := by -- Rewrite target membership as the rescaled sampling statistic lying between the lower and -- upper bootstrap quantiles. Apply random-endpoint portmanteau, the sampling CLT, quantile -- convergence, and `gaussianReal_quantile_Icc_toReal`. letI : IsProbabilityMeasure mu := S.indep.isProbabilityMeasure have hvReal : 0 < h.asymptoticVariance := h.var_pos_finite.1 have hv : 0 < h.asymptoticVariance.toNNReal := by rwa [Real.toNNReal_pos] let G : ProbabilityMeasure ℝ := ⟨gaussianReal 0 h.asymptoticVariance.toNNReal, inferInstance⟩ have hG : h.asymptoticGaussian = G := by apply ProbabilityMeasure.toMeasure_injective rfl let T : ℕ → Omega → ℝ := fun n omega => sqrt (n : ℝ) * (est n (S.sampleVector n omega) - theta0) have hx (n : ℕ) : Measurable (fun omega => S.sampleVector n omega) := by apply measurable_pi_iff.mpr intro i exact S.meas i have hTMeas (n : ℕ) : AEMeasurable (T n) mu := by exact (measurable_const.mul (((h.meas.1 n).comp (hx n)).sub measurable_const)).aemeasurable have hT : Modes.TendstoInLaw (fun _ : ℕ => mu) T atTop (G : Measure ℝ) := by apply (Tendsto_dist_iff _ _ _ hTMeas).2 have hlaw : (fun n => (⟨mu.map (T n), Measure.isProbabilityMeasure_map (hTMeas n)⟩ : ProbabilityMeasure ℝ)) = samplingEstimatorLaw S est theta0 h.meas.1 := by funext n apply ProbabilityMeasure.toMeasure_injective rfl rw [hlaw, ← hG] exact sampling_tendsto_gaussian h have ha2_pos : 0 < alpha / 2 := by linarith have ha2_lt_one : alpha / 2 < 1 := by linarith have h1a2_pos : 0 < 1 - alpha / 2 := by linarith have h1a2_lt_one : 1 - alpha / 2 < 1 := by linarith let qlo : ℕ → Omega → ℝ := fun n omega => bootstrapQuantile est n (alpha / 2) (S.sampleVector n omega) let qhi : ℕ → Omega → ℝ := fun n omega => bootstrapQuantile est n (1 - alpha / 2) (S.sampleVector n omega) let alo := sqrt h.asymptoticVariance * probit (alpha / 2) let ahi := sqrt h.asymptoticVariance * probit (1 - alpha / 2) have hqloMeas (n : ℕ) : Measurable (qlo n) := (measurable_bootstrapQuantile est h.meas.1 n ha2_pos ha2_lt_one).comp (hx n) have hqhiMeas (n : ℕ) : Measurable (qhi n) := (measurable_bootstrapQuantile est h.meas.1 n h1a2_pos h1a2_lt_one).comp (hx n) have hqlo : Modes.TendstoInProbability (fun _ : ℕ => mu) qlo atTop (fun _ _ => alo) := quantile_tendsto h ha2_pos ha2_lt_one have hqhi : Modes.TendstoInProbability (fun _ : ℕ => mu) qhi atTop (fun _ _ => ahi) := quantile_tendsto h h1a2_pos h1a2_lt_one have hprobit : probit (alpha / 2) ≤ probit (1 - alpha / 2) := by apply Causalean.Mathlib.stdNormalCDF_strictMono.le_iff_le.mp rw [stdNormalCDF_probit ha2_pos ha2_lt_one, stdNormalCDF_probit h1a2_pos h1a2_lt_one] linarith have halohi : alo ≤ ahi := by dsimp [alo, ahi] exact mul_le_mul_of_nonneg_left hprobit (Real.sqrt_nonneg _) have hboundary : (G : Measure ℝ) (frontier (Icc alo ahi)) = 0 := by haveI : NullSingletonClass (G : Measure ℝ) := by dsimp [G] exact nullSingletonClass_gaussianReal hv.ne' rw [frontier_Icc halohi] apply le_antisymm · calc (G : Measure ℝ) ({alo, ahi} : Set ℝ) ≤ (G : Measure ℝ) ({alo} : Set ℝ) + (G : Measure ℝ) ({ahi} : Set ℝ) := by rw [show ({alo, ahi} : Set ℝ) = {alo} ∪ {ahi} by ext x; simp [or_comm]] exact measure_union_le ({alo} : Set ℝ) ({ahi} : Set ℝ) _ = 0 := by simp only [NullSingletonClass.measure_singleton, add_zero] · exact zero_le have hcoverage := tendsto_dist_random_Icc_coverage hTMeas hT qlo qhi alo ahi hqloMeas hqhiMeas hqlo hqhi halohi hboundary have hmass : (G : Measure ℝ).real (Icc alo ahi) = 1 - alpha := by dsimp [G, alo, ahi] simpa only [Real.coe_toNNReal _ hvReal.le] using gaussianReal_quantile_Icc_toReal hv halpha0 halpha1 rw [hmass] at hcoverage apply hcoverage.congr' filter_upwards [eventually_ne_atTop 0] with n hn congr 1 ext omega have hsqrt : 0 < sqrt (n : ℝ) := Real.sqrt_pos.2 (by exact_mod_cast Nat.pos_of_ne_zero hn) change (qlo n omega ≤ T n omega ∧ T n omega ≤ qhi n omega) ↔ (est n (S.sampleVector n omega) - qhi n omega / sqrt (n : ℝ) ≤ theta0 ∧ theta0 ≤ est n (S.sampleVector n omega) - qlo n omega / sqrt (n : ℝ)) rw [iff_comm] dsimp [T] constructor · rintro ⟨hlow, hhigh⟩ constructor · have hz : qlo n omega ≤ (est n (S.sampleVector n omega) - theta0) * sqrt (n : ℝ) := (div_le_iff₀ hsqrt).mp (by linarith) nlinarith · have hz : (est n (S.sampleVector n omega) - theta0) * sqrt (n : ℝ) ≤ qhi n omega := (le_div_iff₀ hsqrt).mp (by linarith) nlinarith · rintro ⟨hlow, hhigh⟩ constructor · have hz : est n (S.sampleVector n omega) - theta0 ≤ qhi n omega / sqrt (n : ℝ) := (le_div_iff₀ hsqrt).mpr (by nlinarith) linarith · have hz : qlo n omega / sqrt (n : ℝ) ≤ est n (S.sampleVector n omega) - theta0 := (div_le_iff₀ hsqrt).mpr (by nlinarith) linarith

theorem conditionalSlutsky_cdfKolmogorov_inProb {Omega Y : Type*} [MeasurableSpace Omega] [MeasurableSpace Y] {mu : Measure Omega} [IsProbabilityMeasure mu] (nu : ℕ -> Omega -> ProbabilityMeasure Y) (linear remainder : ℕ -> Omega -> Y -> ℝ) (hlinearMeas : forall n omega, Measurable (linear n omega)) (hremainderMeas : forall n omega, Measurable (remainder n omega)) (G : ProbabilityMeasure ℝ) (hGcont : Continuous (cdf (G : Measure ℝ))) (hlinear : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => cdfKolmogorov ((nu n omega).map (hlinearMeas n omega).aemeasurable) G) atTop (fun _ _ => 0)) (hremainder : forall epsilon : ℝ, 0 < epsilon -> Tendsto (fun n => mu.real {omega | epsilon < ((nu n omega : ProbabilityMeasure Y) : Measure Y).real {y | epsilon < abs (remainder n omega y)}}) atTop (nhds 0)) : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => cdfKolmogorov ((nu n omega).map ((hlinearMeas n omega).add (hremainderMeas n omega)).aemeasurable) G) atTop (fun _ _ => 0) := by rw [Tendsto_inProb_iff, tendstoInMeasure_iff_measureReal_dist] rw [Tendsto_inProb_iff, tendstoInMeasure_iff_measureReal_dist] at hlinear intro epsilon hepsilon let eta := epsilon / 4 have heta : 0 < eta := div_pos hepsilon (by norm_num) have hshift : TendstoUniformly (fun d t : ℝ => cdf (G : Measure ℝ) (t + d)) (cdf (G : Measure ℝ)) (𝓝 0) := by refine tendstoUniformly_of_monotone_of_tendsto_at (F := fun d t : ℝ => cdf (G : Measure ℝ) (t + d)) (G := cdf (G : Measure ℝ)) ?_ ?_ ?_ (monotone_cdf (G : Measure ℝ)) hGcont (tendsto_cdf_atBot (G : Measure ℝ)) (tendsto_cdf_atTop (G : Measure ℝ)) ?_ · intro d x y hxy exact monotone_cdf (G : Measure ℝ) (by linarith) · exact fun d t => cdf_nonneg (G : Measure ℝ) (t + d) · exact fun d t => cdf_le_one (G : Measure ℝ) (t + d) · intro t have htadd : Tendsto (fun d : ℝ => t + d) (𝓝 0) (𝓝 (t + 0)) := tendsto_const_nhds.add tendsto_id change Tendsto ((cdf (G : Measure ℝ)) ∘ fun d : ℝ => t + d) (𝓝 0) (𝓝 (cdf (G : Measure ℝ) t)) simpa only [add_zero] using hGcont.continuousAt.tendsto.comp htadd have hshiftEventually : ∀ᶠ d : ℝ in 𝓝 0, ∀ t : ℝ, dist (cdf (G : Measure ℝ) t) (cdf (G : Measure ℝ) (t + d)) < eta := (Metric.tendstoUniformly_iff.1 hshift) eta heta obtain ⟨r, hr, hball⟩ := Metric.mem_nhds_iff.1 hshiftEventually let delta := min (r / 2) (eta / 2) have hdelta : 0 < delta := lt_min (half_pos hr) (half_pos heta) have hdeltalt : delta < eta := (min_le_right _ _).trans_lt (half_lt_self heta) have hdeltaBall : delta ∈ ball (0 : ℝ) r := by rw [Metric.mem_ball, Real.dist_eq, sub_zero, abs_of_pos hdelta] exact (min_le_left _ _).trans_lt (half_lt_self hr) have hoscPlus (t : ℝ) : dist (cdf (G : Measure ℝ) (t + delta)) (cdf (G : Measure ℝ) t) < eta := by simpa [dist_comm] using hball hdeltaBall t have hoscMinus (t : ℝ) : dist (cdf (G : Measure ℝ) t) (cdf (G : Measure ℝ) (t - delta)) < eta := by simpa [sub_add_cancel, dist_comm] using hoscPlus (t - delta) have hcdfError (rho sigma : ProbabilityMeasure ℝ) (t : ℝ) : |cdf (rho : Measure ℝ) t - cdf (sigma : Measure ℝ) t| ≤ cdfKolmogorov rho sigma := by unfold cdfKolmogorov apply le_ciSup (c := t) refine ⟨1, ?_⟩ rintro _ ⟨x, rfl⟩ rw [abs_le] constructor <;> linarith [cdf_nonneg (rho : Measure ℝ) x, cdf_le_one (rho : Measure ℝ) x, cdf_nonneg (sigma : Measure ℝ) x, cdf_le_one (sigma : Measure ℝ) x] have hupperCDF (n : ℕ) (omega : Omega) (t : ℝ) : cdf (((nu n omega).map ((hlinearMeas n omega).add (hremainderMeas n omega)).aemeasurable) : Measure ℝ) t ≤ cdf (((nu n omega).map (hlinearMeas n omega).aemeasurable) : Measure ℝ) (t + delta) + ((nu n omega : ProbabilityMeasure Y) : Measure Y).real {y | delta < abs (remainder n omega y)} := by rw [cdf_eq_real, cdf_eq_real, measureReal_def, measureReal_def, ProbabilityMeasure.toMeasure_map, ProbabilityMeasure.toMeasure_map, Measure.map_apply_of_aemeasurable ((hlinearMeas n omega).add (hremainderMeas n omega)).aemeasurable measurableSet_Iic, Measure.map_apply_of_aemeasurable (hlinearMeas n omega).aemeasurable measurableSet_Iic] refine (measureReal_mono ?_).trans (measureReal_union_le _ _) intro y hy by_cases htail : delta < abs (remainder n omega y) · exact inr htail · left simp only [mem_preimage, mem_Iic] at hy ⊢ rw [not_lt] at htail have hremLower := (abs_le.mp htail).1 change linear n omega y + remainder n omega y ≤ t at hy linarith have hlowerCDF (n : ℕ) (omega : Omega) (t : ℝ) : cdf (((nu n omega).map (hlinearMeas n omega).aemeasurable) : Measure ℝ) (t - delta) ≤ cdf (((nu n omega).map ((hlinearMeas n omega).add (hremainderMeas n omega)).aemeasurable) : Measure ℝ) t + ((nu n omega : ProbabilityMeasure Y) : Measure Y).real {y | delta < abs (remainder n omega y)} := by rw [cdf_eq_real, cdf_eq_real, measureReal_def, measureReal_def, ProbabilityMeasure.toMeasure_map, ProbabilityMeasure.toMeasure_map, Measure.map_apply_of_aemeasurable (hlinearMeas n omega).aemeasurable measurableSet_Iic, Measure.map_apply_of_aemeasurable ((hlinearMeas n omega).add (hremainderMeas n omega)).aemeasurable measurableSet_Iic] refine (measureReal_mono ?_).trans (measureReal_union_le _ _) intro y hy by_cases htail : delta < abs (remainder n omega y) · exact inr htail · left simp only [mem_preimage, mem_Iic] at hy ⊢ rw [not_lt] at htail have hremUpper := (abs_le.mp htail).2 change linear n omega y + remainder n omega y ≤ t linarith have hsumBound (n : ℕ) (omega : Omega) : cdfKolmogorov ((nu n omega).map ((hlinearMeas n omega).add (hremainderMeas n omega)).aemeasurable) G ≤ cdfKolmogorov ((nu n omega).map (hlinearMeas n omega).aemeasurable) G + ((nu n omega : ProbabilityMeasure Y) : Measure Y).real {y | delta < abs (remainder n omega y)} + eta := by unfold cdfKolmogorov refine ciSup_le fun t => ?_ have hu := hupperCDF n omega t have hl := hlowerCDF n omega t have hePlus := hcdfError ((nu n omega).map (hlinearMeas n omega).aemeasurable) G (t + delta) have heMinus := hcdfError ((nu n omega).map (hlinearMeas n omega).aemeasurable) G (t - delta) unfold cdfKolmogorov at hePlus heMinus rw [abs_le] at hePlus heMinus have hoPlus := hoscPlus t have hoMinus := hoscMinus t rw [Real.dist_eq, abs_lt] at hoPlus hoMinus rw [abs_le] constructor <;> linarith have hsubset (n : ℕ) : {omega | epsilon ≤ dist (cdfKolmogorov ((nu n omega).map ((hlinearMeas n omega).add (hremainderMeas n omega)).aemeasurable) G) 0} ⊆ {omega | eta ≤ dist (cdfKolmogorov ((nu n omega).map (hlinearMeas n omega).aemeasurable) G) 0} ∪ {omega | delta < ((nu n omega : ProbabilityMeasure Y) : Measure Y).real {y | delta < abs (remainder n omega y)}} := by intro omega hbad by_contra hnot rw [mem_union, not_or] at hnot simp only [mem_ofPred_eq] at hbad hnot have hlinNonneg := cdfKolmogorov_nonneg ((nu n omega).map (hlinearMeas n omega).aemeasurable) G have hsumNonneg := cdfKolmogorov_nonneg ((nu n omega).map ((hlinearMeas n omega).add (hremainderMeas n omega)).aemeasurable) G have hlinlt : cdfKolmogorov ((nu n omega).map (hlinearMeas n omega).aemeasurable) G < eta := by rw [Real.dist_eq, sub_zero, abs_of_nonneg hlinNonneg] at hnot exact lt_of_not_ge hnot.1 have htaile : ((nu n omega : ProbabilityMeasure Y) : Measure Y).real {y | delta < abs (remainder n omega y)} ≤ delta := le_of_not_gt hnot.2 have hsumb := hsumBound n omega rw [Real.dist_eq, sub_zero, abs_of_nonneg hsumNonneg] at hbad dsimp [eta] at * linarith have hlinearEta := hlinear eta heta have hremainderDelta := hremainder delta hdelta refine squeeze_zero (fun _ => measureReal_nonneg) (fun n => ?_) (by simpa only [zero_add] using hlinearEta.add hremainderDelta) exact (measureReal_mono (hsubset n)).trans (measureReal_union_le _ _)

theorem quantile_tendsto_inProb_of_cdfKolmogorov {Omega : Type*} [MeasurableSpace Omega] {mu : Measure Omega} [IsProbabilityMeasure mu] (nu : ℕ -> Omega -> ProbabilityMeasure ℝ) (G : ProbabilityMeasure ℝ) {beta : ℝ} (hbeta0 : 0 < beta) (hbeta1 : beta < 1) (hGcont : ContinuousAt (cdf (G : Measure ℝ)) (quantile (G : Measure ℝ) beta)) (hGstrict : StrictlyIncreasingAt (cdf (G : Measure ℝ)) (quantile (G : Measure ℝ) beta)) (hdist : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => cdfKolmogorov (nu n omega) G) atTop (fun _ _ => 0)) : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => quantile (nu n omega : Measure ℝ) beta) atTop (fun _ _ => quantile (G : Measure ℝ) beta) := by rw [Tendsto_inProb_iff, tendstoInMeasure_iff_measureReal_dist] rw [Tendsto_inProb_iff, tendstoInMeasure_iff_measureReal_dist] at hdist intro epsilon hepsilon let q := quantile (G : Measure ℝ) beta have hcdfq : cdf (G : Measure ℝ) q = beta := by exact cdf_quantile_eq (G : Measure ℝ) hbeta0 hbeta1 hGcont have hleft : cdf (G : Measure ℝ) (q - epsilon / 2) < beta := by have h := hGstrict.1 (q - epsilon / 2) (by linarith) simpa only [q, hcdfq] using h have hright : beta < cdf (G : Measure ℝ) (q + epsilon / 2) := by have h := hGstrict.2 (q + epsilon / 2) (by linarith) simpa only [q, hcdfq] using h let delta := min (beta - cdf (G : Measure ℝ) (q - epsilon / 2)) (cdf (G : Measure ℝ) (q + epsilon / 2) - beta) / 2 have hdelta : 0 < delta := by dsimp [delta] positivity have hcdfError (rho sigma : ProbabilityMeasure ℝ) (t : ℝ) : |cdf (rho : Measure ℝ) t - cdf (sigma : Measure ℝ) t| ≤ cdfKolmogorov rho sigma := by unfold cdfKolmogorov apply le_ciSup (c := t) refine ⟨1, ?_⟩ rintro _ ⟨x, rfl⟩ rw [abs_le] constructor <;> linarith [cdf_nonneg (rho : Measure ℝ) x, cdf_le_one (rho : Measure ℝ) x, cdf_nonneg (sigma : Measure ℝ) x, cdf_le_one (sigma : Measure ℝ) x] have hsubset (n : ℕ) : {omega | epsilon ≤ dist (quantile (nu n omega : Measure ℝ) beta) q} ⊆ {omega | delta ≤ dist (cdfKolmogorov (nu n omega) G) 0} := by intro omega hbad have hkolNonneg := cdfKolmogorov_nonneg (nu n omega) G by_contra hnot have hkoldist : dist (cdfKolmogorov (nu n omega) G) 0 < delta := lt_of_not_ge hnot have hkol : cdfKolmogorov (nu n omega) G < delta := by simpa [Real.dist_eq, abs_of_nonneg hkolNonneg] using hkoldist have hdeltaLeft : delta ≤ (beta - cdf (G : Measure ℝ) (q - epsilon / 2)) / 2 := by dsimp [delta] exact div_le_div_of_nonneg_right (min_le_left _ _) (by norm_num) have hdeltaRight : delta ≤ (cdf (G : Measure ℝ) (q + epsilon / 2) - beta) / 2 := by dsimp [delta] exact div_le_div_of_nonneg_right (min_le_right _ _) (by norm_num) have herrorLeft := hcdfError (nu n omega) G (q - epsilon / 2) have herrorRight := hcdfError (nu n omega) G (q + epsilon / 2) rw [abs_le] at herrorLeft herrorRight have hnuLeft : cdf (nu n omega : Measure ℝ) (q - epsilon / 2) < beta := by linarith have hnuRight : beta < cdf (nu n omega : Measure ℝ) (q + epsilon / 2) := by linarith have hqLower : q - epsilon / 2 < quantile (nu n omega : Measure ℝ) beta := by by_contra h have := (quantile_le_iff hbeta0 hbeta1).1 (le_of_not_gt h) linarith have hqUpper : quantile (nu n omega : Measure ℝ) beta ≤ q + epsilon / 2 := (quantile_le_iff hbeta0 hbeta1).2 hnuRight.le change epsilon ≤ dist (quantile (nu n omega : Measure ℝ) beta) q at hbad rw [Real.dist_eq] at hbad have hgood : |quantile (nu n omega : Measure ℝ) beta - q| < epsilon := by rw [abs_lt] constructor <;> linarith linarith have hdistDelta := hdist delta hdelta refine squeeze_zero (fun _ ↦ measureReal_nonneg) (fun n ↦ ?_) hdistDelta exact measureReal_mono (hsubset n)

theorem tendsto_dist_random_Icc_coverage {Omega : Type*} [MeasurableSpace Omega] {mu : Measure Omega} [IsProbabilityMeasure mu] {Sn : ℕ -> Omega -> ℝ} {Q : Measure ℝ} [IsProbabilityMeasure Q] (hSnMeas : forall n, AEMeasurable (Sn n) mu) (hSn : Modes.TendstoInLaw (fun _ : ℕ => mu) Sn atTop Q) (lower upper : ℕ -> Omega -> ℝ) (a b : ℝ) (hlowerMeas : forall n, Measurable (lower n)) (hupperMeas : forall n, Measurable (upper n)) (hlower : Modes.TendstoInProbability (fun _ : ℕ => mu) lower atTop (fun _ _ => a)) (hupper : Modes.TendstoInProbability (fun _ : ℕ => mu) upper atTop (fun _ _ => b)) (hab : a <= b) (hboundary : Q (frontier (Icc a b)) = 0) : Tendsto (fun n => mu.real {omega | Sn n omega ∈ Icc (lower n omega) (upper n omega)}) atTop (nhds (Q.real (Icc a b))) := by let F : ℝ → (ℝ × ℝ) × ℝ := fun z => ((z, a), b) let E : Set ((ℝ × ℝ) × ℝ) := {p | p.1.2 ≤ p.1.1} ∩ {p | p.1.1 ≤ p.2} have hSn' : TendstoInDistribution Sn atTop id (fun _ => mu) Q := by refine { forall_aemeasurable := hSnMeas aemeasurable_limit := measurable_id.aemeasurable tendsto := ?_ } have ht := (Tendsto_dist_iff Sn Q mu hSnMeas).1 hSn simpa only [Measure.map_id] using ht have hSnLower : TendstoInDistribution (fun n omega => (Sn n omega, lower n omega)) atTop (fun z : ℝ => (z, a)) (fun _ => mu) Q := hSn'.prodMk_of_tendstoInMeasure_const Sn lower id hlower (fun n => (hlowerMeas n).aemeasurable) have hJoint : TendstoInDistribution (fun n omega => ((Sn n omega, lower n omega), upper n omega)) atTop F (fun _ => mu) Q := hSnLower.prodMk_of_tendstoInMeasure_const (fun n omega => (Sn n omega, lower n omega)) upper (fun z : ℝ => (z, a)) hupper (fun n => (hupperMeas n).aemeasurable) have hEmeas : MeasurableSet E := by dsimp [E] exact (measurableSet_le (by fun_prop) (by fun_prop)).inter (measurableSet_le (by fun_prop) (by fun_prop)) have hfront_subset : frontier E ⊆ {p : (ℝ × ℝ) × ℝ | p.1.2 = p.1.1} ∪ {p : (ℝ × ℝ) × ℝ | p.1.1 = p.2} := by calc frontier E ⊆ frontier {p : (ℝ × ℝ) × ℝ | p.1.2 ≤ p.1.1} ∪ frontier {p : (ℝ × ℝ) × ℝ | p.1.1 ≤ p.2} := by dsimp [E] exact (frontier_inter_subset _ _).trans (union_subset_union inter_subset_left inter_subset_right) _ ⊆ {p : (ℝ × ℝ) × ℝ | p.1.2 = p.1.1} ∪ {p : (ℝ × ℝ) × ℝ | p.1.1 = p.2} := union_subset_union (frontier_le_subset_eq (by fun_prop) (by fun_prop)) (frontier_le_subset_eq (by fun_prop) (by fun_prop)) have hlimitBoundary : (Q.map F) (frontier E) = 0 := by rw [Measure.map_apply_of_aemeasurable (by fun_prop) isClosed_frontier.measurableSet] refine measure_mono_null ?_ hboundary intro z hz have hz' := hfront_subset hz rw [frontier_Icc hab] simpa [F, eq_comm] using hz' have hport := ProbabilityMeasure.tendsto_measure_of_null_frontier_of_tendsto hJoint.tendsto (E := E) (by change ((Q.map F) (frontier E)).toNNReal = 0 simp [hlimitBoundary]) have hreal := NNReal.continuous_coe.tendsto _ |>.comp hport have hmap_n (n : ℕ) : (mu.map (fun omega => ((Sn n omega, lower n omega), upper n omega))) E = mu {omega | Sn n omega ∈ Icc (lower n omega) (upper n omega)} := by rw [Measure.map_apply_of_aemeasurable (hJoint.forall_aemeasurable n) hEmeas] rfl have hmap_limit : (Q.map F) E = Q (Icc a b) := by rw [Measure.map_apply_of_aemeasurable (by fun_prop) hEmeas] congr 1 simpa [Function.comp_def, ProbabilityMeasure.coe_mk, ENNReal.coe_toNNReal_eq_toReal, Measure.real, hmap_n, hmap_limit] using hreal

theorem bootstrap_tendsto_gaussian (h : BootstrapAsymLinear S est theta0 psi) : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => cdfKolmogorov (bootstrapEstimatorLaw est h.meas.1 n (S.sampleVector n omega)) h.asymptoticGaussian) atTop (fun _ _ => 0) := by letI : IsProbabilityMeasure mu := S.indep.isProbabilityMeasure haveI : IsProbabilityMeasure P := by rw [← S.law] exact Measure.isProbabilityMeasure_map (S.meas 0).aemeasurable have hvar : variance psi P = h.asymptoticVariance := by rw [ProbabilityTheory.variance_eq_integral h.meas.2.aemeasurable, h.mean_zero] simp [asymptoticVariance] have hpop_pos : 0 < populationVariance psi P := by rw [populationVariance, Real.toNNReal_pos, hvar] simpa [asymptoticVariance] using h.var_pos_finite.1 let Gmean : ProbabilityMeasure ℝ := ⟨gaussianReal 0 (populationVariance psi P), inferInstance⟩ have hG : Gmean = h.asymptoticGaussian := by apply ProbabilityMeasure.toMeasure_injective simp [Gmean, asymptoticGaussian, gaussianMeasure, asymptoticVariance, populationVariance, hvar] have hGcont : Continuous (cdf (h.asymptoticGaussian : Measure ℝ)) := by rw [← hG] exact continuous_cdf_gaussianReal_zero hpop_pos have hlinearAE : ∀ᵐ omega ∂mu, Tendsto (fun n => cdfKolmogorov (bootstrapMeanLaw S psi h.meas.2 n omega) h.asymptoticGaussian) atTop (nhds 0) := by filter_upwards [bootstrapMeanLaw_tendsto_gaussian_ae S psi h.meas.2 h.var_pos_finite.2] with omega homega rw [← hG] exact tendsto_cdfKolmogorov_of_tendsto homega tendsto_const_nhds (continuous_cdf_gaussianReal_zero hpop_pos) have hlinearLawMeas (n : ℕ) : Measurable (fun omega => cdfKolmogorov (bootstrapMeanLaw S psi h.meas.2 n omega) h.asymptoticGaussian) := by by_cases hn : n = 0 · subst n simp only [bootstrapMeanLaw, ↓reduceDIte] exact measurable_const · rw [show (fun omega => cdfKolmogorov (bootstrapMeanLaw S psi h.meas.2 n omega) h.asymptoticGaussian) = fun omega => ⨆ q : ℚ, |bootstrapCDF (centeredBootstrapSum psi (S.sampleVector n omega)) (S.sampleVector n omega) q - cdf (h.asymptoticGaussian : Measure ℝ) q| by funext omega rw [cdfKolmogorov_eq_iSup_rat] congr 1 funext q let _ : IsProbabilityMeasure (bootstrapLaw (centeredBootstrapSum psi (S.sampleVector n omega)) (S.sampleVector n omega)) := bootstrapLaw_isProbabilityMeasure _ _ hn (measurable_centeredBootstrapSum psi h.meas.2 _) rw [bootstrapMeanLaw_toMeasure hn, cdf_eq_real, Measure.real_def] rfl] apply Measurable.iSup intro q apply Measurable.abs apply Measurable.sub_const rw [show (fun omega => bootstrapCDF (centeredBootstrapSum psi (S.sampleVector n omega)) (S.sampleVector n omega) (q : ℝ)) = fun omega => ((n : ℝ) ^ n)⁻¹ * ∑ j : Fin n → Fin n, if centeredBootstrapSum psi (S.sampleVector n omega) (S.sampleVector n omega ∘ j) ≤ q then 1 else 0 by funext omega exact bootstrapCDF_eq_average_indicators _ (measurable_centeredBootstrapSum psi h.meas.2 _) _ _] apply Measurable.const_mul apply Finset.measurable_sum intro j hj have hmean : Measurable (fun omega => (n : ℝ)⁻¹ * ∑ k : Fin n, psi (S.Z k omega)) := measurable_const.mul (Finset.measurable_sum _ (fun k _ => h.meas.2.comp (S.meas k))) have hstat : Measurable (fun omega => centeredBootstrapSum psi (S.sampleVector n omega) (S.sampleVector n omega ∘ j)) := by unfold centeredBootstrapSum finMean exact measurable_const.mul (Finset.measurable_sum _ (fun i _ => (h.meas.2.comp (S.meas (j i))).sub hmean)) exact Measurable.ite (measurableSet_le hstat measurable_const) measurable_const measurable_const have hlinearProb : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => cdfKolmogorov (bootstrapMeanLaw S psi h.meas.2 n omega) h.asymptoticGaussian) atTop (fun _ _ => 0) := by exact tendstoInMeasure_of_tendsto_ae (fun n => (hlinearLawMeas n).aestronglyMeasurable) hlinearAE let extend (n : ℕ) (omega : Omega) (xstar : Fin n → X) : ℕ → X := fun i => if hi : i < n then xstar ⟨i, hi⟩ else S.Z i omega have hextMeas (n : ℕ) (omega : Omega) : Measurable (extend n omega) := by apply measurable_pi_iff.mpr intro i unfold extend split_ifs · exact measurable_pi_apply _ · exact measurable_const let nu (n : ℕ) (omega : Omega) : ProbabilityMeasure (ℕ → X) := if hn : n = 0 then (⟨Measure.dirac (fun i => S.Z i omega), Measure.dirac.isProbabilityMeasure⟩ : ProbabilityMeasure (ℕ → X)) else let _ : IsProbabilityMeasure (bootstrapResample (S.sampleVector n omega)) := bootstrapResample_isProbabilityMeasure _ hn (⟨(bootstrapResample (S.sampleVector n omega)).map (extend n omega), Measure.isProbabilityMeasure_map (hextMeas n omega).aemeasurable⟩ : ProbabilityMeasure (ℕ → X)) let linear (n : ℕ) (omega : Omega) (y : ℕ → X) : ℝ := centeredBootstrapSum psi (S.sampleVector n omega) (fun i => y i) let remainder (n : ℕ) (omega : Omega) (y : ℕ → X) : ℝ := centeredEstimatorBootstrapStatistic est n (S.sampleVector n omega) (fun i => y i) - linear n omega y have hlinearMeas (n : ℕ) (omega : Omega) : Measurable (linear n omega) := by unfold linear exact (measurable_centeredBootstrapSum psi h.meas.2 _).comp (measurable_pi_iff.mpr fun i => measurable_pi_apply i.val) have hremainderMeas (n : ℕ) (omega : Omega) : Measurable (remainder n omega) := by unfold remainder centeredEstimatorBootstrapStatistic exact (measurable_const.mul (((h.meas.1 n).comp (measurable_pi_iff.mpr fun i => measurable_pi_apply i.val)).sub measurable_const)).sub (hlinearMeas n omega) have hnueq {n : ℕ} (hn : n ≠ 0) (omega : Omega) : (nu n omega : Measure (ℕ → X)) = (bootstrapResample (S.sampleVector n omega)).map (extend n omega) := by simp [nu, hn] have hlinlaw {n : ℕ} (hn : n ≠ 0) (omega : Omega) : (nu n omega).map (hlinearMeas n omega).aemeasurable = bootstrapMeanLaw S psi h.meas.2 n omega := by apply ProbabilityMeasure.toMeasure_injective change Measure.map (linear n omega) (nu n omega : Measure (ℕ → X)) = _ rw [hnueq hn omega, Measure.map_map (hlinearMeas n omega) (hextMeas n omega), bootstrapMeanLaw_toMeasure hn] unfold bootstrapLaw congr 1 funext xstar simp [linear, extend] have hsumlaw {n : ℕ} (hn : n ≠ 0) (omega : Omega) : (nu n omega).map ((hlinearMeas n omega).add (hremainderMeas n omega)).aemeasurable = bootstrapEstimatorLaw est h.meas.1 n (S.sampleVector n omega) := by apply ProbabilityMeasure.toMeasure_injective change Measure.map (fun y => linear n omega y + remainder n omega y) (nu n omega : Measure (ℕ → X)) = _ rw [hnueq hn omega, Measure.map_map (g := fun y => linear n omega y + remainder n omega y) ((hlinearMeas n omega).add (hremainderMeas n omega)) (hextMeas n omega)] simp [bootstrapEstimatorLaw, hn, bootstrapLaw] congr 1 funext xstar simp [remainder, linear, extend] have hlinearMapped : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => cdfKolmogorov ((nu n omega).map (hlinearMeas n omega).aemeasurable) h.asymptoticGaussian) atTop (fun _ _ => 0) := by apply TendstoInMeasure.congr' (f := fun n omega => cdfKolmogorov (bootstrapMeanLaw S psi h.meas.2 n omega) h.asymptoticGaussian) (g := fun _ => 0) _ (Eventually.of_forall fun _ => rfl) hlinearProb filter_upwards [eventually_ne_atTop 0] with n hn exact Filter.Eventually.of_forall fun omega => congrArg (fun law => cdfKolmogorov law h.asymptoticGaussian) (hlinlaw hn omega).symm have htail (epsilon : ℝ) {n : ℕ} (hn : n ≠ 0) (omega : Omega) : (nu n omega : Measure (ℕ → X)).real {y | epsilon < abs (remainder n omega y)} = (bootstrapResample (S.sampleVector n omega)).real {xstar | epsilon < abs (centeredEstimatorBootstrapStatistic est n (S.sampleVector n omega) xstar - centeredBootstrapSum psi (S.sampleVector n omega) xstar)} := by rw [hnueq hn omega, Measure.real_def, Measure.real_def, Measure.map_apply_of_aemeasurable (hextMeas n omega).aemeasurable (measurableSet_lt measurable_const (hremainderMeas n omega).abs)] congr 2 ext xstar simp only [mem_preimage, mem_ofPred_eq] simp [remainder, linear, extend] have hremainder (epsilon : ℝ) (hepsilon : 0 < epsilon) : Tendsto (fun n => mu.real {omega | epsilon < (nu n omega : Measure (ℕ → X)).real {y | epsilon < abs (remainder n omega y)}}) atTop (nhds 0) := by apply (h.boot_linear epsilon hepsilon).congr' filter_upwards [eventually_ne_atTop 0] with n hn congr 1 ext omega simp only [mem_ofPred_eq] rw [htail epsilon hn omega] have hslutsky := conditionalSlutsky_cdfKolmogorov_inProb nu linear remainder hlinearMeas hremainderMeas h.asymptoticGaussian hGcont hlinearMapped hremainder apply TendstoInMeasure.congr' (f := fun n omega => cdfKolmogorov ((nu n omega).map ((hlinearMeas n omega).add (hremainderMeas n omega)).aemeasurable) h.asymptoticGaussian) (g := fun _ => 0) _ (Eventually.of_forall fun _ => rfl) hslutsky filter_upwards [eventually_ne_atTop 0] with n hn exact Filter.Eventually.of_forall fun omega => congrArg (fun law => cdfKolmogorov law h.asymptoticGaussian) (hsumlaw hn omega)

theorem consistent (h : BootstrapAsymLinear S est theta0 psi) : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => cdfKolmogorov (bootstrapEstimatorLaw est h.meas.1 n (S.sampleVector n omega)) (samplingEstimatorLaw S est theta0 h.meas.1 n)) atTop (fun _ _ => 0) := by -- Pólya applied to `sampling_tendsto_gaussian` gives deterministic Kolmogorov convergence of -- the sampling law to the Gaussian law. Combine it with `bootstrap_tendsto_gaussian` and -- `cdfKolmogorov_triangle` (using symmetry of the distance). letI : IsProbabilityMeasure mu := S.indep.isProbabilityMeasure have hv : 0 < h.asymptoticVariance.toNNReal := by rw [Real.toNNReal_pos] exact h.var_pos_finite.1 have hG : h.asymptoticGaussian = (⟨gaussianReal 0 h.asymptoticVariance.toNNReal, inferInstance⟩ : ProbabilityMeasure ℝ) := by apply ProbabilityMeasure.toMeasure_injective rfl have hGcont : Continuous (cdf (h.asymptoticGaussian : Measure ℝ)) := by rw [hG] exact continuous_cdf_gaussianReal_zero hv have hsampling : Tendsto (fun n => cdfKolmogorov (samplingEstimatorLaw S est theta0 h.meas.1 n) h.asymptoticGaussian) atTop (nhds 0) := tendsto_cdfKolmogorov_of_tendsto (sampling_tendsto_gaussian h) tendsto_const_nhds hGcont have hsamplingProb : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n (_ : Omega) => cdfKolmogorov h.asymptoticGaussian (samplingEstimatorLaw S est theta0 h.meas.1 n)) atTop (fun _ _ => 0) := by apply tendstoInMeasure_of_tendsto_ae · intro n exact aestronglyMeasurable_const · filter_upwards [] with omega apply hsampling.congr' filter_upwards [] with n unfold cdfKolmogorov congr 1 funext t rw [abs_sub_comm] have hsum : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => cdfKolmogorov (bootstrapEstimatorLaw est h.meas.1 n (S.sampleVector n omega)) h.asymptoticGaussian + cdfKolmogorov h.asymptoticGaussian (samplingEstimatorLaw S est theta0 h.meas.1 n)) atTop (fun _ _ => 0) := by have hneg := Modes.TendstoInProbability.comp_continuousAt (g := fun x : ℝ => -x) continuousAt_neg hsamplingProb have hsub := Modes.TendstoInProbability.sub (bootstrap_tendsto_gaussian h) hneg simpa only [Pi.zero_apply, neg_zero, sub_neg_eq_add, sub_self] using hsub rw [Tendsto_inProb_iff] at hsum ⊢ rw [tendstoInMeasure_iff_norm] at hsum ⊢ intro epsilon hepsilon refine tendsto_of_tendsto_of_tendsto_of_le_of_le tendsto_const_nhds (hsum epsilon hepsilon) (fun _ => zero_le) (fun n => ?_) apply measure_mono intro omega homega simp only [mem_setOf_eq, sub_zero, Real.norm_eq_abs] at homega ⊢ have htarget_nonneg := cdfKolmogorov_nonneg (bootstrapEstimatorLaw est h.meas.1 n (S.sampleVector n omega)) (samplingEstimatorLaw S est theta0 h.meas.1 n) have hsum_nonneg : 0 ≤ cdfKolmogorov (bootstrapEstimatorLaw est h.meas.1 n (S.sampleVector n omega)) h.asymptoticGaussian + cdfKolmogorov h.asymptoticGaussian (samplingEstimatorLaw S est theta0 h.meas.1 n) := add_nonneg (cdfKolmogorov_nonneg _ _) (cdfKolmogorov_nonneg _ _) rw [abs_of_nonneg htarget_nonneg] at homega rw [abs_of_nonneg hsum_nonneg] exact homega.trans (cdfKolmogorov_triangle _ _ _)

theorem quantile_tendsto (h : BootstrapAsymLinear S est theta0 psi) {beta : ℝ} (hbeta0 : 0 < beta) (hbeta1 : beta < 1) : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => bootstrapQuantile est n beta (S.sampleVector n omega)) atTop (fun _ _ => sqrt h.asymptoticVariance * probit beta) := by -- Apply `quantile_tendsto_inProb_of_cdfKolmogorov` to `bootstrap_tendsto_gaussian`, use the -- measurable-quantile theorem and the eventual identity with `bootstrapEstimatorLaw`, then -- evaluate the Gaussian quantile by `quantile_gaussianReal_zero`. letI : IsProbabilityMeasure mu := S.indep.isProbabilityMeasure have hvReal : 0 < h.asymptoticVariance := h.var_pos_finite.1 have hv : 0 < h.asymptoticVariance.toNNReal := by rwa [Real.toNNReal_pos] let G : ProbabilityMeasure ℝ := ⟨gaussianReal 0 h.asymptoticVariance.toNNReal, inferInstance⟩ have hG : h.asymptoticGaussian = G := by apply ProbabilityMeasure.toMeasure_injective rfl have hquantLaw : Modes.TendstoInProbability (fun _ : ℕ => mu) (fun n omega => quantile (bootstrapEstimatorLaw est h.meas.1 n (S.sampleVector n omega) : Measure ℝ) beta) atTop (fun _ _ => quantile (G : Measure ℝ) beta) := by apply quantile_tendsto_inProb_of_cdfKolmogorov (nu := fun n omega => bootstrapEstimatorLaw est h.meas.1 n (S.sampleVector n omega)) G hbeta0 hbeta1 · exact (continuous_cdf_gaussianReal_zero hv).continuousAt · have hstrict := strictMono_cdf_gaussianReal_zero hv exact ⟨fun x hx => hstrict hx, fun x hx => hstrict hx⟩ · simpa only [← hG] using bootstrap_tendsto_gaussian h apply TendstoInMeasure.congr' (f := fun n omega => quantile (bootstrapEstimatorLaw est h.meas.1 n (S.sampleVector n omega) : Measure ℝ) beta) (g := fun _ => quantile (G : Measure ℝ) beta) _ _ hquantLaw · filter_upwards [eventually_ne_atTop 0] with n hn exact Filter.Eventually.of_forall fun omega => (bootstrapQuantile_eq_quantile_bootstrapEstimatorLaw est h.meas.1 (bne_iff_ne.mpr hn) beta (S.sampleVector n omega)).symm · exact Filter.Eventually.of_forall fun omega => by dsimp [G] rw [quantile_gaussianReal_zero hv hbeta0 hbeta1] rw [Real.coe_toNNReal _ hvReal.le]

theorem BootstrapAsymLinear.sampleMean (S : IIDSample Omega ℝ mu P) (hX2 : Integrable (fun x : ℝ => x ^ 2) P) (hvar : 0 < variance (fun x : ℝ => x) P) : BootstrapAsymLinear S sampleMeanStatistic (∫ x, x ∂P) (fun x => x - ∫ y, y ∂P) := by letI : IsProbabilityMeasure mu := S.indep.isProbabilityMeasure letI : IsProbabilityMeasure P := by rw [← S.law] exact Measure.isProbabilityMeasure_map (S.meas 0).aemeasurable have hIdL2 : MemLp (fun x : ℝ => x) 2 P := (memLp_two_iff_integrable_sq measurable_id.aestronglyMeasurable).2 hX2 have hIdInt : Integrable (fun x : ℝ => x) P := hIdL2.integrable (by norm_num) have hCenterL2 : MemLp (fun x : ℝ => x - ∫ y, y ∂P) 2 P := by convert hIdL2.sub (memLp_const (∫ y, y ∂P)) using 1 ext x rfl refine { meas := ?_ mean_zero := ?_ var_pos_finite := ?_ linear := ?_ boot_linear := ?_ } · constructor · intro n unfold sampleMeanStatistic finMean fun_prop · fun_prop · rw [integral_sub hIdInt (integrable_const _)] simp · constructor · rw [← ProbabilityTheory.variance_eq_integral (X := fun x : ℝ => x) measurable_id.aemeasurable] exact hvar · exact hCenterL2.integrable_sq · refine TendstoInMeasure.congr' (f := fun _ (_ : Omega) => (0 : ℝ)) (f' := fun (n : ℕ) (omega : Omega) => sqrt (n : ℝ) * (sampleMeanStatistic n (S.sampleVector n omega) - ∫ x, x ∂P) - IsAsymLinear.normalizedSum S (fun x => x - ∫ y, y ∂P) (fun m => range m) n omega) ?_ Filter.EventuallyEq.rfl ?_ · exact Filter.eventually_atTop.2 ⟨1, fun n hn => Filter.Eventually.of_forall fun omega => by have hn0 : n ≠ 0 := Nat.ne_of_gt hn simp only unfold sampleMeanStatistic finMean IIDSample.sampleVector IsAsymLinear.normalizedSum simp only [Finset.card_range, Finset.sum_sub_distrib, Finset.sum_const, nsmul_eq_mul, smul_eq_mul] have hsum : (∑ i : Fin n, S.Z i omega) = ∑ i ∈ range n, S.Z i omega := Fin.sum_univ_eq_sum_range (fun i => S.Z i omega) n rw [hsum] have hnR : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr hn0 have hsqrt : sqrt (n : ℝ) ≠ 0 := Real.sqrt_ne_zero'.mpr (Nat.cast_pos.mpr hn) have hsqrt_sq : (sqrt (n : ℝ)) ^ 2 = (n : ℝ) := Real.sq_sqrt (Nat.cast_nonneg n) field_simp rw [hsqrt_sq] ring⟩ · exact tendstoInMeasure_of_tendsto_ae (fun _ => measurable_const.aestronglyMeasurable) (Filter.Eventually.of_forall fun _ => tendsto_const_nhds) · intro epsilon hepsilon have hev : (fun n => mu.real {omega | epsilon < (bootstrapResample (S.sampleVector n omega)).real {xstar | epsilon < abs (centeredEstimatorBootstrapStatistic sampleMeanStatistic n (S.sampleVector n omega) xstar - centeredBootstrapSum (fun x => x - ∫ y, y ∂P) (S.sampleVector n omega) xstar)}}) =ᶠ[atTop] fun _ : ℕ => (0 : ℝ) := by exact Filter.eventually_atTop.2 ⟨1, fun n hn => by have hn0 : n ≠ 0 := Nat.ne_of_gt hn have hrem (omega : Omega) (xstar : Fin n → ℝ) : centeredEstimatorBootstrapStatistic sampleMeanStatistic n (S.sampleVector n omega) xstar - centeredBootstrapSum (fun x => x - ∫ y, y ∂P) (S.sampleVector n omega) xstar = 0 := by rw [centeredBootstrapSum_eq_sqrt_mul_finAverage_sub hn0 (fun x : ℝ => x - ∫ y, y ∂P) (S.sampleVector n omega) xstar] unfold centeredEstimatorBootstrapStatistic sampleMeanStatistic finMean simp only [Finset.sum_sub_distrib, Finset.sum_const, Finset.card_univ, Fintype.card_fin, nsmul_eq_mul] have hnR : (n : ℝ) ≠ 0 := Nat.cast_ne_zero.mpr hn0 field_simp ring simp_rw [hrem] simp [not_lt_of_ge hepsilon.le]⟩ exact Filter.Tendsto.congr' hev.symm tendsto_const_nhds