ReciprocalGaussian.lean

1 import Mathlib.Analysis.SpecialFunctions.Gaussian.GaussianIntegral
2 
3 /-! Author: Komorebi P. A reciprocal Gaussian substitution for the positive Poisson mode. -/
4 namespace KomorebiPaperVI.ReciprocalGaussian
5 open Set MeasureTheory
6 noncomputable def phi (a x : ℝ) : ℝ := x-a/x
7 
8 theorem phi_derivative (a x : ℝ) (hx : x≠0) :
9     HasDerivAt (phi a) (1+a/x^2) x := by
10   convert! (hasDerivAt_id x).sub ((hasDerivAt_const x a).div (hasDerivAt_id x) hx) using 1
11   dsimp
12   ring
13 
14 theorem phi_injective (a : ℝ) (ha : 0 < a) : InjOn (phi a) (Ioi 0) := by
15   apply StrictMonoOn.injOn
16   intro x hx y hy hxy
17   unfold phi
18   have h : a/y < a/x := div_lt_div_of_pos_left ha hx hxy
19   linarith
20 
21 theorem phi_surjective (a : ℝ) (ha : 0 < a) : phi a '' Ioi 0=univ := by
22   apply eq_univ_of_forall
23   intro y
24   let r := Real.sqrt (y^2+4*a)
25   have hr : r^2=y^2+4*a := Real.sq_sqrt (by positivity)
26   have hr0 : 0 ≤ r := Real.sqrt_nonneg _
27   have hpos : 0 < (y+r)/2 := by nlinarith
28   refine ⟨(y+r)/2,hpos,?_⟩
29   unfold phi
30   apply (sub_eq_iff_eq_add).mpr
31   have hne : y+r≠0 := by linarith
32   field_simp [hne]
33   nlinarith
34 
35 theorem gaussian_jacobian_integrable (a : ℝ) (ha : 0 < a) :
36     IntegrableOn (fun x => (1+a/x^2)*Real.exp (-(phi a x)^2)) (Ioi 0) := by
37   have h := (integrableOn_image_iff_integrableOn_abs_deriv_smul measurableSet_Ioi
38     (fun x hx => (phi_derivative a x hx.ne').hasDerivWithinAt) (phi_injective a ha)
39     (fun u : ℝ => Real.exp (-u^2))).mp
40     (by rw [phi_surjective a ha]; simpa using integrable_exp_neg_mul_sq (b := 1) (by norm_num))
41   apply h.congr
42   filter_upwards [ae_restrict_mem measurableSet_Ioi] with x hx
43   simp only [smul_eq_mul,abs_of_pos (by positivity : 0 < 1+a/x^2)]
44 
45 theorem gaussian_jacobian_integral (a : ℝ) (ha : 0 < a) :
46     (∫ x in Ioi 0,(1+a/x^2)*Real.exp (-(phi a x)^2))=Real.sqrt Real.pi := by
47   have h := integral_image_eq_integral_abs_deriv_smul measurableSet_Ioi
48     (fun x hx => (phi_derivative a x hx.ne').hasDerivWithinAt) (phi_injective a ha)
49     (fun u : ℝ => Real.exp (-u^2))
50   rw [phi_surjective a ha,setIntegral_univ] at h
51   have hg : (∫ u : ℝ,Real.exp (-u^2))=Real.sqrt Real.pi := by
52     simpa using integral_gaussian (1:ℝ)
53   rw [hg] at h
54   rw [h]
55   apply setIntegral_congr_fun measurableSet_Ioi
56   intro x hx
57   simp only [smul_eq_mul,abs_of_pos (by positivity : 0 < 1+a/x^2)]
58 
59 theorem reciprocal_derivative (a x : ℝ) (hx : x≠0) :
60     HasDerivAt (fun u => a/u) (-a/x^2) x := by
61   convert! (hasDerivAt_const x a).div (hasDerivAt_id x) hx using 1
62   dsimp
63   ring
64 
65 theorem reciprocal_image (a : ℝ) (ha : 0 < a) : (fun x => a/x) '' Ioi 0=Ioi 0 := by
66   ext y
67   constructor
68   · rintro ⟨x,hx,rfl⟩
69     exact div_pos ha hx
70   · intro hy
71     refine ⟨a/y,div_pos ha hy,?_⟩
72     field_simp
73 
74 theorem reciprocal_injective (a : ℝ) (ha : 0 < a) : InjOn (fun x => a/x) (Ioi 0) := by
75   intro x hx y hy h
76   change 0 < x at hx
77   change 0 < y at hy
78   have hxy : a*y=a*x := (div_eq_div_iff hx.ne' hy.ne').mp h
79   nlinarith
80 
81 theorem reciprocal_phi (a x : ℝ) (ha : a≠0) (hx : x≠0) : phi a (a/x)= -phi a x := by
82   unfold phi
83   field_simp
84   ring
85 
86 theorem reciprocal_gaussian_symmetry (a : ℝ) (ha : 0 < a) :
87     (∫ x in Ioi 0,a/x^2*Real.exp (-(phi a x)^2))=
88       ∫ x in Ioi 0,Real.exp (-(phi a x)^2) := by
89   have h := integral_image_eq_integral_abs_deriv_smul measurableSet_Ioi
90     (fun x hx => (reciprocal_derivative a x hx.ne').hasDerivWithinAt)
91     (reciprocal_injective a ha) (fun x => Real.exp (-(phi a x)^2))
92   rw [reciprocal_image a ha] at h
93   rw [h]
94   apply setIntegral_congr_fun measurableSet_Ioi
95   intro x hx
96   dsimp only
97   rw [reciprocal_phi a x ha.ne' hx.ne',neg_sq]
98   simp only [smul_eq_mul,abs_div,abs_neg,abs_of_pos ha,abs_of_nonneg (sq_nonneg x)]
99 
100 theorem phi_gaussian_integrable (a : ℝ) (ha : 0 < a) :
101     IntegrableOn (fun x => Real.exp (-(phi a x)^2)) (Ioi 0) := by
102   have hc : ContinuousOn (fun x => Real.exp (-(phi a x)^2)) (Ioi 0) := by
103     intro x hx
104     apply ContinuousAt.continuousWithinAt
105     unfold phi
106     fun_prop (disch := exact hx.ne')
107   apply (gaussian_jacobian_integrable a ha).mono' (hc.aestronglyMeasurable measurableSet_Ioi)
108   filter_upwards [ae_restrict_mem measurableSet_Ioi] with x hx
109   rw [Real.norm_of_nonneg (Real.exp_pos _).le]
110   have h : 0 ≤ a/x^2 := by positivity
111   nlinarith [Real.exp_pos (-(phi a x)^2)]
112 
113 theorem phi_gaussian_integral (a : ℝ) (ha : 0 < a) :
114     (∫ x in Ioi 0,Real.exp (-(phi a x)^2))=Real.sqrt Real.pi/2 := by
115   have hi := phi_gaussian_integrable a ha
116   have hw : IntegrableOn (fun x => a/x^2*Real.exp (-(phi a x)^2)) (Ioi 0) := by
117     convert! (gaussian_jacobian_integrable a ha).sub hi using 1
118     funext x
119     dsimp only [Pi.sub_apply]
120     ring
121   have hj := gaussian_jacobian_integral a ha
122   have he : (fun x => (1+a/x^2)*Real.exp (-(phi a x)^2))=
123       fun x => Real.exp (-(phi a x)^2)+a/x^2*Real.exp (-(phi a x)^2) := by funext x; ring
124   rw [he,integral_add hi hw,reciprocal_gaussian_symmetry a ha] at hj
125   linarith
126 
127 theorem reciprocal_square_gaussian (a : ℝ) (ha : 0 ≤ a) :
128     (∫ x in Ioi 0,Real.exp (-x^2-a^2/x^2))=Real.exp (-2*a)*Real.sqrt Real.pi/2 := by
129   rcases ha.eq_or_lt with he | ha
130   · subst a
131     simpa using integral_gaussian_Ioi (1:ℝ)
132   · have he : EqOn (fun x => Real.exp (-x^2-a^2/x^2))
133         (fun x => Real.exp (-2*a)*Real.exp (-(phi a x)^2)) (Ioi 0) := by
134       intro x hx
135       dsimp only
136       rw [← Real.exp_add]
137       congr 1
138       unfold phi
139       have hx0 : x ≠ 0 := ne_of_gt hx
140       field_simp [hx0]
141       ring
142     rw [setIntegral_congr_fun measurableSet_Ioi he,integral_const_mul,phi_gaussian_integral a ha]
143     ring
144 
145 end KomorebiPaperVI.ReciprocalGaussian

Back to source controls