feat: if X is Gaussian, then X+y and c*X are Gaussian (#7674)

Co-authored-by: Alexander Bentkamp



Co-authored-by: RemyDegenne <remydegenne@gmail.com>
Co-authored-by: Mario Carneiro <di.gama@gmail.com>

Mathlib/Algebra/NeZero.lean

@@ -39,6 +39,9 @@ theorem neZero_iff {n : R} : NeZero n ↔ n ≠ 0 :=
   ⟨fun h ↦ h.out, NeZero.mk⟩
 #align ne_zero_iff neZero_iff
 
+@[simp] lemma neZero_zero_iff_false {α : Type*} [Zero α] : NeZero (0 : α) ↔ False :=
+  ⟨fun h ↦ h.ne rfl, fun h ↦ h.elim⟩
+
 theorem not_neZero {n : R} : ¬NeZero n ↔ n = 0 := by simp [neZero_iff]
 #align not_ne_zero not_neZero
 

Mathlib/MeasureTheory/Function/Jacobian.lean

@@ -1263,4 +1263,72 @@ theorem integral_target_eq_integral_abs_det_fderiv_smul {f : LocalHomeomorph E E
   exact (hf' x hx).hasFDerivWithinAt
 #align measure_theory.integral_target_eq_integral_abs_det_fderiv_smul MeasureTheory.integral_target_eq_integral_abs_det_fderiv_smul
 
+section withDensity
+
+lemma _root_.MeasurableEmbedding.withDensity_ofReal_comap_apply_eq_integral_abs_det_fderiv_mul
+    (hs : MeasurableSet s) (hf : MeasurableEmbedding f)
+    {g : E → ℝ} (hg : ∀ᵐ x ∂μ, x ∈ f '' s → 0 ≤ g x) (hg_int : IntegrableOn g (f '' s) μ)
+    (hf' : ∀ x ∈ s, HasFDerivWithinAt f (f' x) s x) :
+    (μ.withDensity (fun x ↦ ENNReal.ofReal (g x))).comap f s
+      = ENNReal.ofReal (∫ x in s, |(f' x).det| * g (f x) ∂μ) := by
+  rw [Measure.comap_apply f hf.injective (fun t ht ↦ hf.measurableSet_image' ht) _ hs,
+    withDensity_apply _ (hf.measurableSet_image' hs),
+    ← ofReal_integral_eq_lintegral_ofReal hg_int
+      ((ae_restrict_iff' (hf.measurableSet_image' hs)).mpr hg),
+    integral_image_eq_integral_abs_det_fderiv_smul μ hs hf' (hf.injective.injOn _)]
+  simp_rw [smul_eq_mul]
+
+lemma _root_.MeasurableEquiv.withDensity_ofReal_map_symm_apply_eq_integral_abs_det_fderiv_mul
+    (hs : MeasurableSet s) (f : E ≃ᵐ E)
+    {g : E → ℝ} (hg : ∀ᵐ x ∂μ, x ∈ f '' s → 0 ≤ g x) (hg_int : IntegrableOn g (f '' s) μ)
+    (hf' : ∀ x ∈ s, HasFDerivWithinAt f (f' x) s x) :
+    (μ.withDensity (fun x ↦ ENNReal.ofReal (g x))).map f.symm s
+      = ENNReal.ofReal (∫ x in s, |(f' x).det| * g (f x) ∂μ) := by
+  rw [MeasurableEquiv.map_symm,
+    MeasurableEmbedding.withDensity_ofReal_comap_apply_eq_integral_abs_det_fderiv_mul μ hs
+      f.measurableEmbedding hg hg_int hf']
+
+lemma _root_.MeasurableEmbedding.withDensity_ofReal_comap_apply_eq_integral_abs_deriv_mul
+    {f : ℝ → ℝ} (hf : MeasurableEmbedding f) {s : Set ℝ} (hs : MeasurableSet s)
+    {g : ℝ → ℝ} (hg : ∀ᵐ x, x ∈ f '' s → 0 ≤ g x) (hg_int : IntegrableOn g (f '' s))
+    {f' : ℝ → ℝ} (hf' : ∀ x ∈ s, HasDerivWithinAt f (f' x) s x) :
+    (volume.withDensity (fun x ↦ ENNReal.ofReal (g x))).comap f s
+      = ENNReal.ofReal (∫ x in s, |f' x| * g (f x)) := by
+  rw [hf.withDensity_ofReal_comap_apply_eq_integral_abs_det_fderiv_mul volume hs
+    hg hg_int hf']
+  simp only [det_one_smulRight]
+
+lemma _root_.MeasurableEquiv.withDensity_ofReal_map_symm_apply_eq_integral_abs_deriv_mul
+    (f : ℝ ≃ᵐ ℝ) {s : Set ℝ} (hs : MeasurableSet s)
+    {g : ℝ → ℝ} (hg : ∀ᵐ x, x ∈ f '' s → 0 ≤ g x) (hg_int : IntegrableOn g (f '' s))
+    {f' : ℝ → ℝ} (hf' : ∀ x ∈ s, HasDerivWithinAt f (f' x) s x) :
+    (volume.withDensity (fun x ↦ ENNReal.ofReal (g x))).map f.symm s
+      = ENNReal.ofReal (∫ x in s, |f' x| * g (f x)) := by
+  rw [MeasurableEquiv.withDensity_ofReal_map_symm_apply_eq_integral_abs_det_fderiv_mul volume hs
+      f hg hg_int hf']
+  simp only [det_one_smulRight]
+
+lemma _root_.MeasurableEmbedding.withDensity_ofReal_comap_apply_eq_integral_abs_deriv_mul'
+    {f : ℝ → ℝ} (hf : MeasurableEmbedding f) {s : Set ℝ} (hs : MeasurableSet s)
+    {f' : ℝ → ℝ} (hf' : ∀ x, HasDerivAt f (f' x) x)
+    {g : ℝ → ℝ} (hg : 0 ≤ᵐ[volume] g) (hg_int : Integrable g) :
+    (volume.withDensity (fun x ↦ ENNReal.ofReal (g x))).comap f s
+      = ENNReal.ofReal (∫ x in s, |f' x| * g (f x)) :=
+  hf.withDensity_ofReal_comap_apply_eq_integral_abs_deriv_mul hs
+    (by filter_upwards [hg] with x hx using fun _ ↦ hx) hg_int.integrableOn
+    (fun x _ => (hf' x).hasDerivWithinAt)
+
+lemma _root_.MeasurableEquiv.withDensity_ofReal_map_symm_apply_eq_integral_abs_deriv_mul'
+    (f : ℝ ≃ᵐ ℝ) {s : Set ℝ} (hs : MeasurableSet s)
+    {f' : ℝ → ℝ} (hf' : ∀ x, HasDerivAt f (f' x) x)
+    {g : ℝ → ℝ} (hg : 0 ≤ᵐ[volume] g) (hg_int : Integrable g) :
+    (volume.withDensity (fun x ↦ ENNReal.ofReal (g x))).map f.symm s
+      = ENNReal.ofReal (∫ x in s, |f' x| * g (f x)) := by
+  rw [MeasurableEquiv.withDensity_ofReal_map_symm_apply_eq_integral_abs_det_fderiv_mul volume hs
+      f (by filter_upwards [hg] with x hx using fun _ ↦ hx) hg_int.integrableOn
+      (fun x _ => (hf' x).hasDerivWithinAt)]
+  simp only [det_one_smulRight]
+
+end withDensity
+
 end MeasureTheory

Mathlib/MeasureTheory/Measure/AEMeasurable.lean

@@ -44,6 +44,12 @@ theorem aemeasurable_id'' (μ : Measure α) {m : MeasurableSpace α} (hm : m ≤
   @Measurable.aemeasurable α α m0 m id μ (measurable_id'' hm)
 #align probability_theory.ae_measurable_id'' aemeasurable_id''
 
+lemma aemeasurable_of_map_neZero {mβ : MeasurableSpace β} {μ : Measure α}
+    {f : α → β} (h : NeZero (μ.map f)) :
+    AEMeasurable f μ := by
+  by_contra h'
+  simp [h'] at h
+
 namespace AEMeasurable
 
 lemma mono_ac (hf : AEMeasurable f ν) (hμν : μ ≪ ν) : AEMeasurable f μ :=

Mathlib/MeasureTheory/Measure/Dirac.lean

@@ -63,6 +63,16 @@ theorem map_dirac {f : α → β} (hf : Measurable f) (a : α) : (dirac a).map f
   ext fun s hs => by simp [hs, map_apply hf hs, hf hs, indicator_apply]
 #align measure_theory.measure.map_dirac MeasureTheory.Measure.map_dirac
 
+lemma map_const (μ : Measure α) (c : β) : μ.map (fun _ ↦ c) = (μ Set.univ) • dirac c := by
+  ext s hs
+  simp only [aemeasurable_const, measurable_const, smul_toOuterMeasure, OuterMeasure.coe_smul,
+    Pi.smul_apply, dirac_apply' _ hs, smul_eq_mul]
+  classical
+  rw [Measure.map_apply measurable_const hs, Set.preimage_const]
+  by_cases hsc : c ∈ s
+  · rw [(Set.indicator_eq_one_iff_mem _).mpr hsc, mul_one, if_pos hsc]
+  · rw [if_neg hsc, (Set.indicator_eq_zero_iff_not_mem _).mpr hsc, measure_empty, mul_zero]
+
 @[simp]
 theorem restrict_singleton (μ : Measure α) (a : α) : μ.restrict {a} = μ {a} • dirac a := by
   ext1 s hs

Mathlib/MeasureTheory/Measure/MeasureSpace.lean

@@ -1184,6 +1184,7 @@ theorem map_zero (f : α → β) : (0 : Measure α).map f = 0 := by
   by_cases hf : AEMeasurable f (0 : Measure α) <;> simp [map, hf]
 #align measure_theory.measure.map_zero MeasureTheory.Measure.map_zero
 
+@[simp]
 theorem map_of_not_aemeasurable {f : α → β} {μ : Measure α} (hf : ¬AEMeasurable f μ) :
     μ.map f = 0 := by simp [map, hf]
 #align measure_theory.measure.map_of_not_ae_measurable MeasureTheory.Measure.map_of_not_aemeasurable

Mathlib/Probability/Distributions/Gaussian.lean

@@ -1,7 +1,7 @@
 /-
 Copyright (c) 2023 Rémy Degenne. All rights reserved.
 Released under Apache 2.0 license as described in the file LICENSE.
-Authors: Lorenzo Luccioli, Rémy Degenne
+Authors: Lorenzo Luccioli, Rémy Degenne, Alexander Bentkamp
 -/
 import Mathlib.Analysis.SpecialFunctions.Gaussian
 import Mathlib.Probability.Notation
@@ -21,6 +21,13 @@ We define a Gaussian measure over the reals.
   If `v = 0`, this is `dirac μ`, otherwise it is defined as the measure with density
   `gaussianPdf μ v` with respect to the Lebesgue measure.
 
+## Main results
+
+* `gaussianReal_add_const`: if `X` is a random variable with Gaussian distribution with mean `μ` and
+  variance `v`, then `X + y` is Gaussian with mean `μ + y` and variance `v`.
+* `gaussianReal_const_mul`: if `X` is a random variable with Gaussian distribution with mean `μ` and
+  variance `v`, then `c * X` is Gaussian with mean `c * μ` and variance `c^2 * v`.
+
 -/
 
 open scoped ENNReal NNReal Real
@@ -114,6 +121,40 @@ lemma integral_gaussianPdfReal_eq_one (μ : ℝ) {v : ℝ≥0} (hv : v ≠ 0) :
     (ae_of_all _ (gaussianPdfReal_nonneg _ _)), ← ENNReal.ofReal_one] at h
   rwa [← ENNReal.ofReal_eq_ofReal_iff (integral_nonneg (gaussianPdfReal_nonneg _ _)) zero_le_one]
 
+lemma gaussianPdfReal_sub {μ : ℝ} {v : ℝ≥0} (x y : ℝ) :
+    gaussianPdfReal μ v (x - y) = gaussianPdfReal (μ + y) v x := by
+  simp only [gaussianPdfReal]
+  rw [sub_add_eq_sub_sub_swap]
+
+lemma gaussianPdfReal_add {μ : ℝ} {v : ℝ≥0} (x y : ℝ) :
+    gaussianPdfReal μ v (x + y) = gaussianPdfReal (μ - y) v x := by
+  rw [sub_eq_add_neg, ← gaussianPdfReal_sub, sub_eq_add_neg, neg_neg]
+
+lemma gaussianPdfReal_inv_mul {μ : ℝ} {v : ℝ≥0} {c : ℝ} (hc : c ≠ 0) (x : ℝ) :
+    gaussianPdfReal μ v (c⁻¹ * x) = |c| * gaussianPdfReal (c * μ) (⟨c^2, sq_nonneg _⟩ * v) x := by
+  simp only [gaussianPdfReal._eq_1, gt_iff_lt, zero_lt_two, zero_le_mul_left, NNReal.zero_le_coe,
+    Real.sqrt_mul', one_div, mul_inv_rev, NNReal.coe_mul, NNReal.coe_mk, NNReal.coe_pos]
+  rw [← mul_assoc]
+  refine congr_arg₂ _ ?_ ?_
+  · field_simp
+    rw [Real.sqrt_sq_eq_abs]
+    ring_nf
+    calc (Real.sqrt ↑v)⁻¹ * (Real.sqrt 2)⁻¹ * (Real.sqrt π)⁻¹
+      = (Real.sqrt ↑v)⁻¹ * (Real.sqrt 2)⁻¹ * (Real.sqrt π)⁻¹ * (|c| * |c|⁻¹) := by
+          rw [mul_inv_cancel, mul_one]
+          simp only [ne_eq, abs_eq_zero, hc, not_false_eq_true]
+    _ = (Real.sqrt ↑v)⁻¹ * (Real.sqrt 2)⁻¹ * (Real.sqrt π)⁻¹ * |c| * |c|⁻¹ := by ring
+  · congr 1
+    field_simp
+    congr 1
+    ring
+
+lemma gaussianPdfReal_mul {μ : ℝ} {v : ℝ≥0} {c : ℝ} (hc : c ≠ 0) (x : ℝ) :
+    gaussianPdfReal μ v (c * x)
+      = |c⁻¹| * gaussianPdfReal (c⁻¹ * μ) (⟨(c^2)⁻¹, inv_nonneg.mpr (sq_nonneg _)⟩ * v) x := by
+  conv_lhs => rw [← inv_inv c, gaussianPdfReal_inv_mul (inv_ne_zero hc)]
+  simp
+
 /-- The pdf of a Gaussian distribution on ℝ with mean `μ` and variance `v`. -/
 noncomputable
 def gaussianPdf (μ : ℝ) (v : ℝ≥0) (x : ℝ) : ℝ≥0∞ := ENNReal.ofReal (gaussianPdfReal μ v x)
@@ -192,6 +233,124 @@ lemma rnDeriv_gaussianReal (μ : ℝ) (v : ℝ≥0) :
   · rw [gaussianReal_of_var_ne_zero _ hv]
     exact Measure.rnDeriv_withDensity _ (measurable_gaussianPdf μ v)
 
+section Transformations
+
+variable {μ : ℝ} {v : ℝ≥0}
+
+lemma _root_.MeasurableEmbedding.gaussianReal_comap_apply (hv : v ≠ 0)
+    {f : ℝ → ℝ} (hf : MeasurableEmbedding f)
+    {f' : ℝ → ℝ} (h_deriv : ∀ x, HasDerivAt f (f' x) x) {s : Set ℝ} (hs : MeasurableSet s) :
+    (gaussianReal μ v).comap f s
+      = ENNReal.ofReal (∫ x in s, |f' x| * gaussianPdfReal μ v (f x)) := by
+  rw [gaussianReal_of_var_ne_zero _ hv, gaussianPdf_def]
+  exact hf.withDensity_ofReal_comap_apply_eq_integral_abs_deriv_mul' hs h_deriv
+    (ae_of_all _ (gaussianPdfReal_nonneg _ _)) (integrable_gaussianPdfReal _ _)
+
+lemma _root_.MeasurableEquiv.gaussianReal_map_symm_apply (hv : v ≠ 0) (f : ℝ ≃ᵐ ℝ) {f' : ℝ → ℝ}
+    (h_deriv : ∀ x, HasDerivAt f (f' x) x) {s : Set ℝ} (hs : MeasurableSet s) :
+    (gaussianReal μ v).map f.symm s
+      = ENNReal.ofReal (∫ x in s, |f' x| * gaussianPdfReal μ v (f x)) := by
+  rw [gaussianReal_of_var_ne_zero _ hv, gaussianPdf_def]
+  exact f.withDensity_ofReal_map_symm_apply_eq_integral_abs_deriv_mul' hs h_deriv
+    (ae_of_all _ (gaussianPdfReal_nonneg _ _)) (integrable_gaussianPdfReal _ _)
+
+/-- The map of a Gaussian distribution by addition of a constant is a Gaussian. -/
+lemma gaussianReal_map_add_const (y : ℝ) :
+    (gaussianReal μ v).map (· + y) = gaussianReal (μ + y) v := by
+  by_cases hv : v = 0
+  · simp only [hv, ne_eq, not_true, gaussianReal_zero_var]
+    exact Measure.map_dirac (measurable_id'.add_const _) _
+  let e : ℝ ≃ᵐ ℝ := (Homeomorph.addRight y).symm.toMeasurableEquiv
+  have he' : ∀ x, HasDerivAt e ((fun _ ↦ 1) x) x := fun _ ↦ (hasDerivAt_id _).sub_const y
+  change (gaussianReal μ v).map e.symm = gaussianReal (μ + y) v
+  ext s' hs'
+  rw [MeasurableEquiv.gaussianReal_map_symm_apply hv e he' hs']
+  simp only [abs_neg, abs_one, MeasurableEquiv.coe_mk, Equiv.coe_fn_mk, one_mul, ne_eq]
+  rw [gaussianReal_apply_eq_integral _ hv hs']
+  simp [gaussianPdfReal_sub _ y, Homeomorph.addRight, ← sub_eq_add_neg]
+
+/-- The map of a Gaussian distribution by addition of a constant is a Gaussian. -/
+lemma gaussianReal_map_const_add (y : ℝ) :
+    (gaussianReal μ v).map (y + ·) = gaussianReal (μ + y) v := by
+  simp_rw [add_comm y]
+  exact gaussianReal_map_add_const y
+
+/-- The map of a Gaussian distribution by multiplication by a constant is a Gaussian. -/
+lemma gaussianReal_map_const_mul (c : ℝ) :
+    (gaussianReal μ v).map (c * ·) = gaussianReal (c * μ) (⟨c^2, sq_nonneg _⟩ * v) := by
+  by_cases hv : v = 0
+  · simp only [hv, mul_zero, ne_eq, not_true, gaussianReal_zero_var]
+    exact Measure.map_dirac (measurable_id'.const_mul c) μ
+  by_cases hc : c = 0
+  · simp only [hc, zero_mul, ne_eq, abs_zero, mul_eq_zero]
+    rw [Measure.map_const]
+    simp only [ne_eq, measure_univ, one_smul, mul_eq_zero]
+    convert (gaussianReal_zero_var 0).symm
+    simp only [ne_eq, zero_pow', mul_eq_zero, hv, or_false]
+    rfl
+  let e : ℝ ≃ᵐ ℝ := (Homeomorph.mulLeft₀ c hc).symm.toMeasurableEquiv
+  have he' : ∀ x, HasDerivAt e ((fun _ ↦ c⁻¹) x) x := by
+    suffices ∀ x, HasDerivAt (fun x => c⁻¹ * x) (c⁻¹ * 1) x by rwa [mul_one] at this
+    exact fun _ ↦ HasDerivAt.const_mul _ (hasDerivAt_id _)
+  change (gaussianReal μ v).map e.symm = gaussianReal (c * μ) (⟨c^2, sq_nonneg _⟩ * v)
+  ext s' hs'
+  rw [MeasurableEquiv.gaussianReal_map_symm_apply hv e he' hs']
+  simp only [MeasurableEquiv.coe_mk, Equiv.coe_fn_mk, ne_eq, mul_eq_zero]
+  rw [gaussianReal_apply_eq_integral _ _ hs']
+  swap
+  · simp only [ne_eq, mul_eq_zero, hv, or_false]
+    rw [← NNReal.coe_eq]
+    simp [hc]
+  simp only [Homeomorph.mulLeft₀, Equiv.toFun_as_coe, Equiv.mulLeft₀_apply, Equiv.invFun_as_coe,
+    Equiv.mulLeft₀_symm_apply, Homeomorph.toMeasurableEquiv_coe, Homeomorph.homeomorph_mk_coe_symm,
+    Equiv.coe_fn_symm_mk, gaussianPdfReal_inv_mul hc]
+  congr with x
+  suffices |c⁻¹| * |c| = 1 by rw [← mul_assoc, this, one_mul]
+  rw [abs_inv, inv_mul_cancel]
+  rwa [ne_eq, abs_eq_zero]
+
+/-- The map of a Gaussian distribution by multiplication by a constant is a Gaussian. -/
+lemma gaussianReal_map_mul_const (c : ℝ) :
+    (gaussianReal μ v).map (· * c) = gaussianReal (c * μ) (⟨c^2, sq_nonneg _⟩ * v) := by
+  simp_rw [mul_comm _ c]
+  exact gaussianReal_map_const_mul c
+
+variable {Ω : Type} [MeasureSpace Ω]
+
+/-- If `X` is a real random variable with Gaussian law with mean `μ` and variance `v`, then `X + y`
+has Gaussian law with mean `μ + y` and variance `v`. -/
+lemma gaussianReal_add_const {X : Ω → ℝ} (hX : Measure.map X ℙ = gaussianReal μ v) (y : ℝ) :
+    Measure.map (fun ω ↦ X ω + y) ℙ = gaussianReal (μ + y) v := by
+  have hXm : AEMeasurable X := aemeasurable_of_map_neZero (by rw [hX]; infer_instance)
+  change Measure.map ((fun ω ↦ ω + y) ∘ X) ℙ = gaussianReal (μ + y) v
+  rw [← AEMeasurable.map_map_of_aemeasurable (measurable_id'.add_const _).aemeasurable hXm, hX,
+    gaussianReal_map_add_const y]
+
+/-- If `X` is a real random variable with Gaussian law with mean `μ` and variance `v`, then `y + X`
+has Gaussian law with mean `μ + y` and variance `v`. -/
+lemma gaussianReal_const_add {X : Ω → ℝ} (hX : Measure.map X ℙ = gaussianReal μ v) (y : ℝ) :
+    Measure.map (fun ω ↦ y + X ω) ℙ = gaussianReal (μ + y) v := by
+  simp_rw [add_comm y]
+  exact gaussianReal_add_const hX y
+
+/-- If `X` is a real random variable with Gaussian law with mean `μ` and variance `v`, then `c * X`
+has Gaussian law with mean `c * μ` and variance `c^2 * v`. -/
+lemma gaussianReal_const_mul {X : Ω → ℝ} (hX : Measure.map X ℙ = gaussianReal μ v) (c : ℝ) :
+    Measure.map (fun ω ↦ c * X ω) ℙ = gaussianReal (c * μ) (⟨c^2, sq_nonneg _⟩ * v) := by
+  have hXm : AEMeasurable X := aemeasurable_of_map_neZero (by rw [hX]; infer_instance)
+  change Measure.map ((fun ω ↦ c * ω) ∘ X) ℙ = gaussianReal (c * μ) (⟨c^2, sq_nonneg _⟩ * v)
+  rw [← AEMeasurable.map_map_of_aemeasurable (measurable_id'.const_mul c).aemeasurable hXm, hX]
+  exact gaussianReal_map_const_mul c
+
+/-- If `X` is a real random variable with Gaussian law with mean `μ` and variance `v`, then `X * c`
+has Gaussian law with mean `c * μ` and variance `c^2 * v`. -/
+lemma gaussianReal_mul_const {X : Ω → ℝ} (hX : Measure.map X ℙ = gaussianReal μ v) (c : ℝ) :
+    Measure.map (fun ω ↦ X ω * c) ℙ = gaussianReal (c * μ) (⟨c^2, sq_nonneg _⟩ * v) := by
+  simp_rw [mul_comm _ c]
+  exact gaussianReal_const_mul hX c
+
+end Transformations
+
 end GaussianReal
 
 end ProbabilityTheory