Stat.Nonparametric.Moment­Problems.Finite­Moment­Near­Gaussian­Perturbation

theorem hamburgerCarlemanSeries_eq (ν : Measure ℝ) : hamburgerCarlemanSeries ν = ∑' s : ℕ, (ofReal |∫ x, x ^ (2 * (s + 1)) ∂ν|).rpow (-(1 : ℝ) / (2 * (s + 1))) := rfl

theorem hamburgerCarlemanSeries_eq_top_of_evenMoment_le (ν : Measure ℝ) (hmoment : ∀ n : ℕ, 0 < n → |rawMoment ν (2 * n)| ≤ 2 * (2 * n : ℝ) ^ n) : hamburgerCarlemanSeries ν = ⊤ := by -- Put `n = s + 1`. Antitonicity of `rpow` at the negative exponent turns the -- moment upper bound into a lower bound for each summand. Bound the extra -- factor `2^(-1/(2n))` uniformly below and compare with the divergent real -- p-series `n^(-1/2)`, transported to `ENNReal`. let f : ℕ → NNReal := fun s => (1 / 2 : NNReal) * ((2 * (s + 1) : ℕ) : NNReal) ^ (-(1 : ℝ) / 2) have hf : ¬ Summable f := by let p : ℝ := -(1 : ℝ) / 2 let c : NNReal := (1 / 2 : NNReal) * (2 : NNReal) ^ p have hfull : ¬ Summable (fun n : ℕ => (n : NNReal) ^ p) := by rw [NNReal.summable_rpow] norm_num [p] have hshift : ¬ Summable (fun s : ℕ => ((s + 1 : ℕ) : NNReal) ^ p) := by intro hs apply hfull apply (NNReal.summable_nat_add_iff 1).mp simpa [Nat.add_comm] using hs have hc : c ≠ 0 := by positivity intro hfs apply hshift apply (summable_mul_left_iff hc).mp simpa [f, c, p, NNReal.mul_rpow, mul_assoc] using hfs have hftop : (∑' s : ℕ, (f s : ENNReal)) = ⊤ := by apply ENNReal.tsum_coe_eq_top_iff_not_summable_coe.mpr rwa [NNReal.summable_coe] rw [hamburgerCarlemanSeries] apply top_unique rw [← hftop] apply ENNReal.tsum_le_tsum intro s let n := s + 1 have hn : 0 < n := by omega let m := |rawMoment ν (2 * n)| have hmnonneg : 0 ≤ m := abs_nonneg _ have hterm : (f s : ENNReal) ≤ (ofReal m).rpow (-(1 : ℝ) / (2 * n)) := by by_cases hmzero : m = 0 · have he : -(1 : ℝ) / (2 * n) < 0 := div_neg_of_neg_of_pos (by norm_num) (by positivity) simp [hmzero, he] · have hmpos : 0 < m := lt_of_le_of_ne hmnonneg (Ne.symm hmzero) have hreal := real_half_evenPow_rpow_le n hn m hmpos (hmoment n hn) rw [show (ofReal m).rpow (-(1 : ℝ) / (2 * n)) = ofReal (m ^ (-(1 : ℝ) / (2 * n))) by exact ENNReal.ofReal_rpow_of_pos hmpos] rw [ENNReal.coe_nnreal_eq] apply ENNReal.ofReal_le_ofReal simpa [f, n] using hreal simpa [n, m] using hterm

theorem sourceCumulant_eq_of_rawMoment_eq_up_to (μ ν : Measure ℝ) (K : ℕ) (hm : ∀ j, j ≤ K → rawMoment μ j = rawMoment ν j) : ∀ k, k ≤ K → sourceCumulant μ (id : ℝ → ℝ) k = sourceCumulant ν (id : ℝ → ℝ) k := by intro k hk rw [sourceCumulant_eq_cumFromMom, sourceCumulant_eq_cumFromMom] exact cumFromMom_congr_up_to hm hk

theorem gaussianPerturbation_spec {K : ℕ} {h : ℝ → ℝ} {ε rho : ℝ} (hmeas : Measurable h) (hbound : ∀ x, |h x| ≤ 1) (hnonzero : 0 < ∫ x, |h x| ∂gaussianReal 0 1) (horth : ∀ k, k ≤ K → Integrable (fun x => x ^ k * h x) (gaussianReal 0 1) ∧ (∫ x, x ^ k * h x ∂gaussianReal 0 1) = 0) (hεpos : 0 < ε) (hεone : ε < 1) (hερho : ε < rho) : let F := gaussianPerturbation h ε IsProbabilityMeasure F ∧ F ≠ gaussianReal 0 1 ∧ tvDist F (gaussianReal 0 1) < rho ∧ (∀ k, k ≤ K → rawMoment F k = rawMoment (gaussianReal 0 1) k) ∧ (∀ k, Integrable (fun x : ℝ => |x| ^ k) F) ∧ (∀ n : ℕ, 0 < n → |rawMoment F (2 * n)| ≤ 2 * (2 * n : ℝ) ^ n) := by let μ : Measure ℝ := gaussianReal 0 1 let d : ℝ → ℝ := fun x => 1 + ε * h x let f : ℝ → ℝ≥0∞ := fun x => ofReal (d x) have hd_pos (x : ℝ) : 0 < d x := by have hx := (abs_le.mp (hbound x)).1 dsimp [d] nlinarith have hd_nonneg (x : ℝ) : 0 ≤ d x := (hd_pos x).le have hd_le (x : ℝ) : d x ≤ 2 := by have hx := (abs_le.mp (hbound x)).2 dsimp [d] nlinarith have hf_meas : Measurable f := by exact (measurable_const.add (measurable_const.mul hmeas)).ennreal_ofReal have hf_top : ∀ᵐ x ∂μ, f x < ∞ := by exact Filter.Eventually.of_forall fun x => by simp [f] have hf_toReal (x : ℝ) : (f x).toReal = d x := by simp [f, ENNReal.toReal_ofReal (hd_nonneg x)] obtain ⟨hh_int0, hh_mean0⟩ := horth 0 (zero_le K) have hh_int : Integrable h μ := by simpa [μ] using hh_int0 have habsh_int : Integrable (fun x => |h x|) μ := by simpa [Real.norm_eq_abs] using hh_int.norm have hh_mean : ∫ x, h x ∂μ = 0 := by simpa [μ] using hh_mean0 let F : Measure ℝ := gaussianPerturbation h ε have hF_def : F = μ.withDensity f := by rfl have hprob : IsProbabilityMeasure F := by dsimp [F, gaussianPerturbation, μ] change IsProbabilityMeasure (tiltMeasure (gaussianReal 0 1) h ε) apply isProbabilityMeasure_tiltMeasure hmeas hbound · simpa [μ] using hh_mean · rw [abs_of_pos hεpos] linarith letI : IsProbabilityMeasure F := hprob have hpower_int (k : ℕ) : Integrable (fun x : ℝ => x ^ k) μ := by apply (integrable_norm_iff (by fun_prop)).mp simpa [μ, Real.norm_eq_abs, abs_pow] using standardGaussian_abs_pow_integrable k have habs_mul_h_int (k : ℕ) : Integrable (fun x : ℝ => |x| ^ k * h x) μ := by apply (standardGaussian_abs_pow_integrable k).mono' · exact ((measurable_id.abs.pow_const k).mul hmeas).aestronglyMeasurable · filter_upwards [] with x rw [Real.norm_eq_abs, abs_mul, abs_of_nonneg (pow_nonneg (abs_nonneg x) k)] exact mul_le_of_le_one_right (by positivity) (hbound x) have habs_F_int (k : ℕ) : Integrable (fun x : ℝ => |x| ^ k) F := by rw [hF_def, integrable_withDensity_iff hf_meas hf_top] have hi : Integrable (fun x : ℝ => |x| ^ k + ε * (|x| ^ k * h x)) μ := (standardGaussian_abs_pow_integrable k).add ((habs_mul_h_int k).const_mul ε) apply hi.congr filter_upwards [] with x rw [hf_toReal] dsimp [d] ring have hmoment_match (k : ℕ) (hk : k ≤ K) : rawMoment F k = rawMoment μ k := by obtain ⟨hkh_int, hkh_zero⟩ := horth k hk rw [rawMoment_eq_integral, rawMoment_eq_integral, hF_def, integral_withDensity_eq_integral_toReal_smul hf_meas hf_top] have heq : (fun x : ℝ => (f x).toReal • x ^ k) = fun x => x ^ k + ε * (x ^ k * h x) := by funext x rw [hf_toReal] dsimp [d] change (1 + ε * h x) * x ^ k = x ^ k + ε * (x ^ k * h x) ring rw [heq, integral_add (hpower_int k) (hkh_int.const_mul ε), integral_const_mul, hkh_zero] ring have hne : F ≠ μ := by intro heq have hdensEq : μ.withDensity f = μ.withDensity (fun _ : ℝ => (1 : ℝ≥0∞)) := by calc μ.withDensity f = F := hF_def.symm _ = μ := heq _ = μ.withDensity (fun _ : ℝ => (1 : ℝ≥0∞)) := by symm change μ.withDensity (1 : ℝ → ℝ≥0∞) = μ exact withDensity_one have hfdens : f =ᵐ[μ] (1 : ℝ → ℝ≥0∞) := (withDensity_eq_iff_of_sigmaFinite hf_meas.aemeasurable measurable_const.aemeasurable).mp hdensEq have hh_zero : h =ᵐ[μ] 0 := by filter_upwards [hfdens] with x hx have hx' := congrArg toReal hx rw [hf_toReal] at hx' simp only [Pi.one_apply, ENNReal.toReal_one] at hx' dsimp [d] at hx' have : h x = 0 := by nlinarith simpa [this] have : ∫ x, |h x| ∂μ = 0 := by apply integral_eq_zero_of_ae filter_upwards [hh_zero] with x hx simp [hx] exact (by simpa [μ, this] using hnonzero : False) have htv : tvDist F μ < rho := by rw [tvDist] apply lt_of_le_of_lt (ciSup_le fun A => ?_) hερho have hAint : Integrable (A.1.indicator h) μ := hh_int.indicator A.2 have hevent : F.real A.1 - μ.real A.1 = ε * ∫ x in A.1, h x ∂μ := by have hFI : F.real A.1 = ∫ x, A.1.indicator (fun _ => (1 : ℝ)) x ∂F := by rw [integral_indicator A.2, setIntegral_const] simp [measureReal_def] have hμI : μ.real A.1 = ∫ x, A.1.indicator (fun _ => (1 : ℝ)) x ∂μ := by rw [integral_indicator A.2, setIntegral_const] simp [measureReal_def] rw [hFI, hμI, hF_def, integral_withDensity_eq_integral_toReal_smul hf_meas hf_top] have hi : (fun x : ℝ => (f x).toReal • A.1.indicator (fun _ => (1 : ℝ)) x) = fun x => A.1.indicator (fun _ => (1 : ℝ)) x + ε * A.1.indicator h x := by funext x rw [hf_toReal] by_cases hx : x ∈ A.1 <;> simp [d, hx] <;> ring rw [hi, integral_add ((integrable_const 1).indicator A.2) (hAint.const_mul ε), integral_const_mul] rw [← integral_indicator A.2] ring rw [abs_sub_comm, ← neg_sub, hevent, abs_neg, abs_mul, abs_of_pos hεpos] calc ε * |∫ x in A.1, h x ∂μ| ≤ ε * ∫ x in A.1, |h x| ∂μ := by gcongr exact abs_integral_le_integral_abs _ ≤ ε * ∫ x, |h x| ∂μ := by exact mul_le_mul_of_nonneg_left (setIntegral_le_integral habsh_int (Filter.Eventually.of_forall fun x => abs_nonneg (h x))) hεpos.le _ ≤ ε := by have hint_le : ∫ x, |h x| ∂μ ≤ 1 := by calc _ ≤ ∫ _x, (1 : ℝ) ∂μ := integral_mono habsh_int (integrable_const 1) hbound _ = 1 := by simp [μ] nlinarith refine ⟨hprob, ?_, htv, hmoment_match, habs_F_int, ?_⟩ · simpa [F, μ] using hne · intro n _hn rw [rawMoment_eq_integral] calc |∫ x : ℝ, x ^ (2 * n) ∂F| ≤ ∫ x : ℝ, |x ^ (2 * n)| ∂F := abs_integral_le_integral_abs _ = ∫ x : ℝ, |x| ^ (2 * n) ∂F := by congr 1 funext x rw [abs_pow] _ = ∫ x : ℝ, d x * |x| ^ (2 * n) ∂μ := by rw [hF_def, integral_withDensity_eq_integral_toReal_smul hf_meas hf_top] apply integral_congr_ae filter_upwards [] with x rw [hf_toReal] simp [smul_eq_mul, mul_comm] _ ≤ ∫ x : ℝ, 2 * |x| ^ (2 * n) ∂μ := by apply integral_mono · have hi : Integrable (fun x : ℝ => d x * |x| ^ (2 * n)) μ := by apply ((standardGaussian_abs_pow_integrable (2 * n)).const_mul 2).mono' · fun_prop · filter_upwards [] with x rw [Real.norm_eq_abs, abs_mul, abs_of_pos (hd_pos x), abs_of_nonneg (pow_nonneg (abs_nonneg x) (2 * n))] exact mul_le_mul_of_nonneg_right (hd_le x) (by positivity) exact hi · exact (standardGaussian_abs_pow_integrable (2 * n)).const_mul 2 · intro x exact mul_le_mul_of_nonneg_right (hd_le x) (by positivity) _ = 2 * ∫ x : ℝ, x ^ (2 * n) ∂μ := by rw [integral_const_mul] congr 2 funext x rw [← abs_pow, abs_of_nonneg (pow_nonneg (even_two_mul n) x)] _ ≤ 2 * (2 * n : ℝ) ^ n := by gcongr simpa [μ] using standardGaussian_evenMoment_le n

theorem exists_finiteMoment_near_gaussian_perturbation (K : ℕ) (rho : ℝ) (hK : 3 ≤ K) (hrho : 0 < rho) : ∃ F : Measure ℝ, IsProbabilityMeasure F ∧ (∫ x, x ∂F) = 0 ∧ variance id F = 1 ∧ ¬ IsGaussianLaw F ∧ tvDist F (gaussianReal 0 1) < rho ∧ (∀ k, k ≤ K → (∫ x, x ^ k ∂F) = ∫ x, x ^ k ∂gaussianReal 0 1) ∧ (∀ k, Integrable (fun x : ℝ => |x| ^ k) F) ∧ (∑' s : ℕ, (ofReal |∫ x, x ^ (2 * (s + 1)) ∂F|).rpow (-(1 : ℝ) / (2 * (s + 1)))) = ⊤ ∧ (∀ k, k ≤ K → sourceCumulant F (id : ℝ → ℝ) k = sourceCumulant (gaussianReal 0 1) (id : ℝ → ℝ) k) := by obtain ⟨h, hmeas, hbound, hnonzero, horth⟩ := exists_bounded_gaussian_orthogonal_perturbation K let ε : ℝ := min (1 / 2) (rho / 2) have hεpos : 0 < ε := by dsimp [ε] exact lt_min (by norm_num) (half_pos hrho) have hεone : ε < 1 := by have hεle : ε ≤ 1 / 2 := min_le_left _ _ linarith have hερho : ε < rho := by have hεle : ε ≤ rho / 2 := min_le_right _ _ linarith let F : Measure ℝ := gaussianPerturbation h ε have hspec := gaussianPerturbation_spec hmeas hbound hnonzero horth hεpos hεone hερho change IsProbabilityMeasure F ∧ F ≠ gaussianReal 0 1 ∧ tvDist F (gaussianReal 0 1) < rho ∧ (∀ k, k ≤ K → rawMoment F k = rawMoment (gaussianReal 0 1) k) ∧ (∀ k, Integrable (fun x : ℝ => |x| ^ k) F) ∧ (∀ n : ℕ, 0 < n → |rawMoment F (2 * n)| ≤ 2 * (2 * n : ℝ) ^ n) at hspec rcases hspec with ⟨hprob, hne, htv, hmoment, habs, heven⟩ have hmean : (∫ x, x ∂F) = 0 := by have hmoment_one := hmoment 1 (by omega) simpa [rawMoment_eq_integral] using hmoment_one have hvariance : variance id F = 1 := by have hmoment_two := hmoment 2 (by omega) have hgauss_two : rawMoment (gaussianReal 0 1) 2 = 1 := by have hv : variance id (gaussianReal 0 1) = (1 : ℝ) := by simpa using (variance_id_gaussianReal (μ := (0 : ℝ)) (v := (1 : NNReal))) change variance (fun x : ℝ => x) (gaussianReal 0 1) = 1 at hv rw [variance_eq_integral measurable_id'.aemeasurable] at hv simpa [rawMoment_eq_integral, Function.id_def] using hv change variance (fun x : ℝ => x) F = 1 rw [variance_eq_integral measurable_id'.aemeasurable, hmean] simpa [rawMoment_eq_integral] using hmoment_two.trans hgauss_two have hnongaussian : ¬ IsGaussianLaw F := by rintro ⟨m, v, hF⟩ have hm : m = 0 := by calc m = ∫ x, x ∂gaussianReal m v := integral_id_gaussianReal.symm _ = ∫ x, x ∂F := by rw [hF] _ = 0 := hmean have hv_real : (v : ℝ) = 1 := by calc (v : ℝ) = variance id (gaussianReal m v) := variance_id_gaussianReal.symm _ = variance id F := by rw [hF] _ = 1 := hvariance have hv : v = 1 := by exact_mod_cast hv_real apply hne exact hF.trans (gaussianReal_ext_iff.mpr ⟨hm, hv⟩) have hcarleman : hamburgerCarlemanSeries F = ⊤ := hamburgerCarlemanSeries_eq_top_of_evenMoment_le F heven refine ⟨F, hprob, hmean, hvariance, hnongaussian, htv, ?_, habs, ?_, ?_⟩ · intro k hk simpa only [rawMoment_eq_integral] using hmoment k hk · rw [← hamburgerCarlemanSeries_eq] exact hcarleman · exact sourceCumulant_eq_of_rawMoment_eq_up_to F (gaussianReal 0 1) K hmoment

theorem exists_bounded_gaussian_orthogonal_perturbation (K : ℕ) : ∃ h : ℝ → ℝ, Measurable h ∧ (∀ x, |h x| ≤ 1) ∧ 0 < ∫ x, |h x| ∂gaussianReal 0 1 ∧ ∀ k, k ≤ K → Integrable (fun x => x ^ k * h x) (gaussianReal 0 1) ∧ (∫ x, x ^ k * h x ∂gaussianReal 0 1) = 0 := by classical let A : Matrix (Fin (K + 1)) (Fin (K + 2)) ℝ := fun j i => ∫ x in gaussianBlock i, x ^ (j : ℕ) ∂gaussianReal 0 1 let T : (Fin (K + 2) → ℝ) →ₗ[ℝ] (Fin (K + 1) → ℝ) := A.mulVecLin have hT_not_injective : ¬ Injective T := by intro hT have hdim := LinearMap.finrank_le_finrank_of_injective hT simp only [Module.finrank_pi, Fintype.card_fin] at hdim omega obtain ⟨c, d, hTcd, hcd⟩ := Function.not_injective_iff.mp hT_not_injective let v : Fin (K + 2) → ℝ := c - d have hv_ne : v ≠ 0 := sub_ne_zero.mpr hcd have hTv : T v = 0 := by dsimp [v] rw [map_sub, hTcd, sub_self] have hv_norm_pos : 0 < ‖v‖ := norm_pos_iff.mpr hv_ne let a : Fin (K + 2) → ℝ := fun i => v i / ‖v‖ let h : ℝ → ℝ := fun x => ∑ i : Fin (K + 2), (gaussianBlock i).indicator (fun _ => a i) x have ha_le (i : Fin (K + 2)) : |a i| ≤ 1 := by have hvi : |v i| ≤ ‖v‖ := by rw [Pi.norm_def] exact_mod_cast Finset.le_sup (f := fun j => ‖v j‖₊) (Finset.mem_univ i) change |v i / ‖v‖| ≤ 1 rw [abs_div, abs_of_nonneg (norm_nonneg v)] exact (div_le_one hv_norm_pos).mpr hvi have hh_meas : Measurable h := by dsimp [h] apply measurable_sum intro i _ exact measurable_const.indicator (gaussianBlock_measurable i) have hh_bound : ∀ x, |h x| ≤ 1 := by intro x by_cases hx : ∃ i : Fin (K + 2), x ∈ gaussianBlock i · obtain ⟨i, hxi⟩ := hx change |∑ j : Fin (K + 2), (gaussianBlock j).indicator (fun _ => a j) x| ≤ 1 rw [gaussianBlock_indicator_sum_eq a i hxi] exact ha_le i · have hz : h x = 0 := by dsimp [h] apply Finset.sum_eq_zero intro i _ rw [Set.indicator_apply, if_neg] exact fun hxi => hx ⟨i, hxi⟩ simp [hz] have hh_integrable : Integrable (fun x => |h x|) (gaussianReal 0 1) := by apply (integrable_const (μ := gaussianReal 0 1) (1 : ℝ)).mono' · exact hh_meas.abs.aestronglyMeasurable · filter_upwards [] with x simpa only [Real.norm_eq_abs, abs_abs] using hh_bound x have hh_pos : 0 < ∫ x, |h x| ∂gaussianReal 0 1 := by obtain ⟨i, hvi⟩ : ∃ i, v i ≠ 0 := by simpa [Function.ne_iff] using hv_ne let g : ℝ → ℝ := (gaussianBlock i).indicator (fun _ => |a i|) have hg_integrable : Integrable g (gaussianReal 0 1) := by exact (integrable_const (μ := gaussianReal 0 1) |a i|).indicator (gaussianBlock_measurable i) have hg_le : ∀ x, g x ≤ |h x| := by intro x by_cases hx : x ∈ gaussianBlock i · simp only [g, Set.indicator_of_mem hx] change |a i| ≤ |∑ j : Fin (K + 2), (gaussianBlock j).indicator (fun _ => a j) x| rw [gaussianBlock_indicator_sum_eq a i hx] · simp [g, hx] have hg_pos : 0 < ∫ x, g x ∂gaussianReal 0 1 := by rw [integral_indicator (gaussianBlock_measurable i), setIntegral_const] exact mul_pos (ENNReal.toReal_pos (ne_of_gt (gaussianBlock_measure_pos i)) (MeasureTheory.measure_ne_top _ _)) (abs_pos.mpr (div_ne_zero hvi (ne_of_gt hv_norm_pos))) exact hg_pos.trans_le (integral_mono hg_integrable hh_integrable hg_le) refine ⟨h, hh_meas, hh_bound, hh_pos, ?_⟩ intro k hk have hterm (i : Fin (K + 2)) : Integrable (fun x => a i * (gaussianBlock i).indicator (fun x : ℝ => x ^ k) x) (gaussianReal 0 1) := (integrable_gaussianBlock_pow i).const_mul (a i) have hfun : (fun x => x ^ k * h x) = fun x => ∑ i : Fin (K + 2), a i * (gaussianBlock i).indicator (fun x : ℝ => x ^ k) x := by funext x dsimp [h] rw [Finset.mul_sum] apply Finset.sum_congr rfl intro i _ by_cases hxi : x ∈ gaussianBlock i · simp [hxi, mul_comm] · simp [hxi] have hint : Integrable (fun x => x ^ k * h x) (gaussianReal 0 1) := by rw [hfun] simpa only [Finset.sum_const_zero, Finset.sum_attach] using integrable_finsetSum univ (fun i _ => hterm i) refine ⟨hint, ?_⟩ rw [hfun, integral_finsetSum univ (fun i _ => hterm i)] simp_rw [integral_const_mul, integral_indicator (gaussianBlock_measurable _)] have hj : k < K + 1 := by omega have hTvj : (T v) (⟨k, hj⟩ : Fin (K + 1)) = 0 := by rw [hTv] rfl change (∑ i : Fin (K + 2), (v i / ‖v‖) * A ⟨k, hj⟩ i) = 0 calc (∑ i : Fin (K + 2), (v i / ‖v‖) * A ⟨k, hj⟩ i) = (1 / ‖v‖) * (T v) ⟨k, hj⟩ := by simp only [T, Matrix.mulVecLin_apply, mulVec, dotProduct] rw [Finset.mul_sum] apply Finset.sum_congr rfl intro i _ ring _ = 0 := by rw [hTvj, mul_zero]