diff --git a/Statlib.lean b/Statlib.lean index c38d427..8ec626c 100644 --- a/Statlib.lean +++ b/Statlib.lean @@ -1,2 +1,7 @@ import Statlib.Inference import Statlib.QMD +import Statlib.Tweedie.CompoundPoisson +import Statlib.Tweedie.CompoundPoissonCore +import Statlib.Tweedie.GammaConvolution +import Statlib.Tweedie.Tweedie +import Statlib.Tweedie.TweedieAux diff --git a/Statlib/Tweedie/CompoundPoisson.lean b/Statlib/Tweedie/CompoundPoisson.lean new file mode 100644 index 0000000..280f6f3 --- /dev/null +++ b/Statlib/Tweedie/CompoundPoisson.lean @@ -0,0 +1,399 @@ +/- +Copyright (c) 2026 Bjørn Kjos-Hanssen. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Bjørn Kjos-Hanssen, Jireh Loreaux +-/ +import Mathlib +import Statlib.Tweedie.Tweedie +import Statlib.Tweedie.GammaConvolution +import Statlib.Tweedie.CompoundPoissonCore +import Statlib.Tweedie.TweedieAux +/-! The compound-Poisson construction lives in `CompoundPoissonCore`. -/ +open scoped Real +open scoped Pointwise +open scoped NNReal ENNReal +set_option maxRecDepth 4000 +set_option synthInstance.maxSize 128 +set_option relaxedAutoImplicit false +set_option autoImplicit false +set_option grind.warning false + +namespace CompoundPoisson +open MeasureTheory ProbabilityTheory +variable (μ : Measure ℝ) [IsProbabilityMeasure μ] (lam : ℝ≥0) +/-! ### The mass of `{0}` under the Poisson count -/ +lemma poissonMeasure_singleton_zero : + poissonMeasure lam {0} = ENNReal.ofReal (Real.exp (-lam)) := by + unfold poissonMeasure + rw [Measure.sum_apply] + · simp [Pi.single, Function.update] + · exact measurableSet_singleton 0 + +lemma poissonMeasure_singleton_zero_pos : 0 < poissonMeasure lam {0} := by + rw [poissonMeasure_singleton_zero] + simp [ENNReal.ofReal_pos, Real.exp_pos] + +/-! ### The atomless key lemmas -/ + +/- +For a nonempty finite index set, the product of an atomless probability measure assigns zero +mass to the hyperplane `{∑ i, g i = 0}`. This is the statement that an `n`-fold convolution of an +atomless measure (`n ≥ 1`) is again atomless at `0`. +-/ +set_option maxHeartbeats 8000000 in +-- times out +open Classical in +lemma pi_sum_eq_zero {ι : Type*} [Fintype ι] [Nonempty ι] [NoAtoms μ] : + Measure.pi (fun _ : ι => μ) {g : ι → ℝ | ∑ i, g i = 0} = 0 := by + -- Let $j$ be an arbitrary element of $\iota$. + obtain ⟨j, hj⟩ : ∃ j : ι, True := by + exact ⟨Classical.arbitrary ι, trivial⟩ + set p : ι → Prop := fun i => i ≠ j + set t := {q : (({i // p i}) → ℝ) × (({i // ¬ p i}) → ℝ) | (∑ i, q.1 i) + (∑ i, q.2 i) = 0} + -- By the measure-preserving property, the measure of the preimage of $t$ under $e$ is + -- equal to the measure of $t$ under the product measure. + have h_preimage : (Measure.pi (fun _ : ι => μ)) {g : ι → ℝ | ∑ i, g i = 0} = + (Measure.prod (Measure.pi (fun _ : {i // p i} => μ)) + (Measure.pi (fun _ : {i // ¬p i} => μ))) t := by + have h_preimage : (Measure.pi (fun _ : ι => μ)) {g : ι → ℝ | ∑ i, g i = 0} = + (Measure.pi (fun _ : ι => μ)) (Set.preimage (MeasurableEquiv.piEquivPiSubtypeProd + (fun _ : ι => ℝ) p) t) := by + congr with g + simp only [MeasurableEquiv.piEquivPiSubtypeProd_apply, Set.mem_setOf_eq, t] + rw [show ( Finset.univ : Finset ι ) = Finset.image ( fun i : { i // p i } => i.val ) + Finset.univ ∪ { j } from ?_, Finset.sum_union] <;> norm_num [Finset.sum_image, p] + · rw [show ( Finset.univ : Finset { i // ¬p i } ) = { ⟨j, by aesop⟩ } from + Finset.eq_singleton_iff_unique_mem.mpr ⟨Finset.mem_univ _, fun x hx => by aesop⟩ + ]; simp [p] + · ext i; by_cases hi : i = j <;> aesop + rw [h_preimage, MeasureTheory.MeasurePreserving.measure_preimage] + · convert MeasureTheory.measurePreserving_piEquivPiSubtypeProd ( fun _ : ι => μ ) p using 1 + · refine MeasurableSet.nullMeasurableSet ?_ + exact measurableSet_eq_fun ( by measurability ) ( by measurability ) + -- For fixed `a`, `Prod.mk a ⁻¹' t = {b : {i//¬p i} → ℝ | (∑ i, a i) + (∑ i, b i) = 0}`. + -- Since `{i // ¬ p i} = {i // i = j}` is a singleton type with unique element + -- `j₀ := ⟨j, by simp⟩`, `∑ i, b i = b j₀`. + -- Hence the slice is `{b | b j₀ = -(∑ i, a i)}`. + have h_slice : ∀ a : ({i // p i} → ℝ), (Measure.pi (fun _ : {i // ¬p i} => μ)) + {b : {i // ¬p i} → ℝ | (∑ i, a i) + (∑ i, b i) = 0} = 0 := by + intro a + have h_singleton : {b : {i // ¬p i} → ℝ | (∑ i, a i) + (∑ i, b i) = 0} + = {fun _ => -(∑ i, a i)} := by + ext b + simp only [ne_eq, Set.mem_setOf_eq, Set.mem_singleton_iff, p] + rw [show ( Finset.univ : Finset { i // ¬p i } ) = { ⟨j, by aesop⟩ } from + Finset.eq_singleton_iff_unique_mem.mpr ⟨Finset.mem_univ _, + fun x hx => Subtype.ext <| by aesop⟩] + simp only [Finset.sum_singleton, add_eq_zero_iff_eq_neg] + exact ⟨fun h => funext fun x => by rw [ + show x = ⟨j, by aesop⟩ from Subtype.ext <| by aesop]; linarith, fun h => by rw [h]; simp⟩ + aesop + rw [h_preimage, MeasureTheory.Measure.prod_apply] + · aesop + · exact measurableSet_eq_fun ( by measurability ) ( by measurability ) +/- +The infinite product measure assigns zero mass to the event that the first `m+1` jumps sum to +zero, whenever `μ` is atomless. +-/ +lemma infinitePi_sum_range_succ_eq_zero [NoAtoms μ] (m : ℕ) : + Measure.infinitePi (fun _ : ℕ => μ) + {f : ℕ → ℝ | ∑ i ∈ Finset.range (m + 1), f i = 0} = 0 := by + -- Let's denote the finite set $\{0, 1, ..., m\}$ by $I$. + set I : Finset ℕ := Finset.range (m + 1) with hI_def + -- By `MeasureTheory.Measure.infinitePi_map_restrict (fun _ : ℕ => μ)`, we have + -- `Measure.map I.restrict (Measure.infinitePi (fun _ : ℕ => μ)) = Measure.pi (fun i : ↥I => μ)`. + have h_map : Measure.map (fun f : ℕ → ℝ => fun i : I => f i) (Measure.infinitePi (fun _ : ℕ => μ)) + = Measure.pi (fun i : I => μ) := by + convert MeasureTheory.Measure.infinitePi_map_restrict ( fun _ : ℕ => μ ) using 1 + convert congr_arg ( fun μ => μ { g : I → ℝ | ∑ i, g i = 0 } ) h_map using 1 + · rw [MeasureTheory.Measure.map_apply] + · congr! 2 + ext; simp [Finset.sum_attach] + · exact measurable_pi_lambda _ fun _ => measurable_pi_apply _ + · exact measurableSet_eq_fun ( by measurability ) ( by measurability ) + · convert Eq.symm ( pi_sum_eq_zero μ ) using 1 + exact ⟨0, Finset.mem_range.mpr ( Nat.succ_pos _ )⟩ +/-! ### The sanity check -/ +/- +The probability that the compound Poisson variable is `0` equals `ℙ(N = 0)`, provided the jump +distribution `μ` is atomless. +-/ +theorem compoundPoisson_singleton_zero [NoAtoms μ] : + compoundPoisson μ lam {(0 : ℝ)} = poissonMeasure lam {0} := by + have h_integral : ∫⁻ n, (Measure.infinitePi (fun _ : ℕ => μ)) + {f : ℕ → ℝ | ∑ i ∈ Finset.range n, f i = 0} ∂(poissonMeasure lam) + = (poissonMeasure lam) {0} := by + rw [MeasureTheory.lintegral_congr_ae, MeasureTheory.lintegral_indicator_one] + · exact MeasurableSingletonClass.measurableSet_singleton _ + · filter_upwards [] with n + cases n <;> simp_all [infinitePi_sum_range_succ_eq_zero] + convert h_integral using 1 + rw [compoundPoisson, baseMeasure, MeasureTheory.Measure.map_apply] + · convert MeasureTheory.Measure.prod_apply _ using 1 + · infer_instance + · exact measurable_jumpSum ( MeasurableSingletonClass.measurableSet_singleton _ ) + · exact measurable_jumpSum + · exact MeasurableSingletonClass.measurableSet_singleton _ +/-- The probability that the compound Poisson variable is `0` is positive (assuming `μ` atomless), +and in fact equals `e^{-lam}`. -/ +theorem compoundPoisson_singleton_zero_pos [NoAtoms μ] : + 0 < compoundPoisson μ lam {(0 : ℝ)} := by + rw [compoundPoisson_singleton_zero μ lam] + exact poissonMeasure_singleton_zero_pos lam +end CompoundPoisson + +open CompoundPoisson +open GammaConv +-- Now we can construct Tweedie: +/- +The compound-Poisson construction of the Tweedie distribution uses a +`Gamma` jump distribution with **shape** `α = (2 - p)/(p - 1)` and **rate** +`γ = μ^(1-p)/(φ (p-1))` (equivalently scale `φ (p-1) μ^(p-1)`), and Poisson rate +`λ = μ^(2-p)/(φ (2-p))`. + +In Mathlib, `ProbabilityTheory.gammaMeasure a r` takes its **first** argument `a` as the shape +and its **second** argument `r` as the rate (density `r^a/Γ a * x^(a-1) * e^{-r x}`). +-/ + +/-- Realization of the Tweedie distribution as a compound Poisson distribution +for `1 < p < 2`. -/ +noncomputable def tweedieConstruction (μ : ℝ) (hμ : 0 ≤ μ) {φ p : ℝ} + (hφ : 0 ≤ φ) + (hp₂ : p < 2) := + compoundPoisson + (ProbabilityTheory.gammaMeasure + ((2 - p) / (p - 1)) -- shape α + (μ ^ (1 - p) / (φ * (p - 1))) -- rate γ + ) + ⟨μ^(2-p) / (φ * (2-p)), by positivity⟩ -- Poisson rate + +instance (a r : ℝ) : MeasureTheory.NoAtoms + (ProbabilityTheory.gammaMeasure a r) := by + refine { measure_singleton := ?_ } + intro x + simp [ProbabilityTheory.gammaMeasure] + +/-- The gamma distribution with Tweedie parameters is a probability measure. -/ +theorem tweedie_gamma_IsProbabilityMeasure {μ : ℝ} (hμ : 0 < μ) {φ p : ℝ} + (hφ : 0 < φ) (hp₁ : 1 < p) + (hp₂ : p < 2) : MeasureTheory.IsProbabilityMeasure + (ProbabilityTheory.gammaMeasure + ((2 - p) / (p - 1)) (μ ^ (1 - p) / (φ * (p - 1)))) := + ProbabilityTheory.isProbabilityMeasure_gammaMeasure + (div_pos (by linarith) (by linarith)) + (div_pos (Real.rpow_pos_of_pos hμ _) (by nlinarith)) + +/-- Show that the `tweedie_construction` agrees +with the Tweedie distribution at 0. -/ +lemma tweedie_zero_sanity_check (μ : ℝ) (hμ : 0 < μ) {φ p : ℝ} + (hφ : 0 < φ) (hp₁ : 1 < p) + (hp₂ : p < 2) : tweedieConstruction μ (by linarith) + (show 0 ≤ φ by linarith) + hp₂ {0} = tweedieProbZero μ φ p := by + unfold tweedieConstruction + rw [@compoundPoisson_singleton_zero _ + (tweedie_gamma_IsProbabilityMeasure hμ hφ hp₁ hp₂)] + unfold tweedieProbZero + rw [poissonMeasure_singleton_zero] + norm_cast + have : NNReal.toReal ⟨μ ^ (2 - p) / (φ * (2 - p)), + by positivity⟩ + = μ ^ (2 - p) / (φ * (2 - p)) := rfl + rw [this] + have (x : ℝ) (hx : 0 ≤ x) : + ENNReal.ofReal x = ENNReal.ofNNReal (⟨x,hx⟩ : {y // (0:ℝ) ≤ y}) := by + refine (ENNReal.toReal_eq_toReal_iff' ?_ ?_).mp ?_ + · simp + · simp + exact ENNReal.toReal_ofReal hx + rw [this] + · simp + congr + field_simp + · exact Real.exp_nonneg (-(μ ^ (2 - p) / (φ * (2 - p)))) + +/-- The Tweedie measure, which was introduced without explanation, +equals the Tweedie construction, which is compound Poisson, +at least on the set {0}. +-/ +lemma tweedie_zero_sanity_check₂ (μ : ℝ) (hμ : 0 < μ) {φ p : ℝ} + (hφ : 0 < φ) (hp₁ : 1 < p) + (hp₂ : p < 2) : + let hμ₀ := le_of_lt hμ + let hφ₀ := le_of_lt hφ + tweedieConstruction μ hμ₀ hφ₀ hp₂ {0} = + tweedieMeasure μ hφ₀ hp₁ hp₂ {0} := by + intro hμ₀ hφ₀ + rw [tweedie_zero_sanity_check μ hμ hφ hp₁ hp₂] + simp [tweedieMeasure] + +open MeasureTheory ProbabilityTheory in +/-- Closed form of a `twG` summand as a Poisson-weight times a gamma density, valid for `y > 0`. +This is the pointwise correspondence between the `j`-th term of the Tweedie density series and the +`j`-fold-convolution (gamma) term of the compound-Poisson law. -/ +lemma twG_eq_gamma (μ : ℝ) (hμ : 0 < μ) {φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hφ : 0 < φ) + (j : ℕ) {y : ℝ} (hy : 0 < y) : + twG μ φ p j y + = Real.exp (-twZ μ φ p) * (twZ μ φ p) ^ j / (j.factorial) + * gammaPDFReal ((j : ℝ) * ((2 - p) / (p - 1))) (μ ^ (1 - p) / (φ * (p - 1))) y := by + have h1p : (1 - p) ≠ 0 := by linarith + have hpm1 : (p - 1) ≠ 0 := by linarith + rcases Nat.eq_zero_or_pos j with hj0 | hjpos + · subst hj0 + rw [twG_zero, gammaPDFReal, if_pos hy.le] + simp [Real.Gamma_zero] + · have hjR : (0 : ℝ) < (j : ℝ) := by exact_mod_cast hjpos + have hAneg : (2 - p) / (1 - p) < 0 := div_neg_of_pos_of_neg (by linarith) (by linarith) + have hargpos : 0 < -(j : ℝ) * ((2 - p) / (1 - p)) := + mul_pos_of_neg_of_neg (by simpa using hjR) hAneg + have hΓ : Real.Gamma (-(j : ℝ) * ((2 - p) / (1 - p))) ≠ 0 := + ne_of_gt (Real.Gamma_pos_of_pos hargpos) + rw [tw_pt μ hp₁ hφ j hy, gammaPDFReal, if_pos hy.le] + have hAB : (j : ℝ) * ((2 - p) / (p - 1)) = -(j : ℝ) * ((2 - p) / (1 - p)) := by + field_simp; ring + have hdenom : φ ^ ((j : ℝ) * (1 - (2 - p) / (1 - p))) * (2 - p) ^ j ≠ 0 := by positivity + have key := TweedieAux.tw_scalar_identity μ φ p hμ hφ hp₁ hp₂ j + rw [hAB] at key + rw [eq_div_iff hdenom] at key + rw [hAB, ← key] + simp only [twZ] + have hfac : (j.factorial : ℝ) ≠ 0 := by positivity + have hphirpow : φ ^ ((j : ℝ) * (1 - (2 - p) / (1 - p))) ≠ 0 := by positivity + have h2pj : ((2 - p) : ℝ) ^ j ≠ 0 := by positivity + have hphi2p : φ * (2 - p) ≠ 0 := by positivity + field_simp [hΓ, hfac, hphirpow, h2pj, hphi2p] + +open MeasureTheory ProbabilityTheory in +/-- The compound-Poisson law of a Gamma jump distribution splits as an atom at `0` (mass +`e^{-lam}`) plus an absolutely continuous part whose density is the Poisson-weighted sum of the +convolution (gamma) densities. -/ +lemma compound_split (α γ : ℝ) (hα : 0 < α) (hγ : 0 < γ) (lam : ℝ≥0) : + CompoundPoisson.compoundPoisson (gammaMeasure α γ) lam + = ENNReal.ofReal (Real.exp (-(lam : ℝ))) • Measure.dirac 0 + + volume.withDensity (fun z => ∑' n : ℕ, + ENNReal.ofReal (Real.exp (-(lam : ℝ)) * (lam : ℝ) ^ (n + 1) / ((n + 1).factorial)) + * gammaPDF (((n + 1 : ℕ) : ℝ) * α) γ z) := by + haveI : IsProbabilityMeasure (gammaMeasure α γ) := isProbabilityMeasure_gammaMeasure hα hγ + have hgPDFmeas : ∀ a : ℝ, Measurable (fun z => gammaPDF a γ z) := fun a => + (measurable_gammaPDFReal a γ).ennreal_ofReal + set d : ℝ → ℝ≥0∞ := fun z => ∑' n : ℕ, + ENNReal.ofReal (Real.exp (-(lam : ℝ)) * (lam : ℝ) ^ (n + 1) / ((n + 1).factorial)) + * gammaPDF (((n + 1 : ℕ) : ℝ) * α) γ z with hd + have hd_meas : Measurable d := by + apply Measurable.tsum + intro n + exact (hgPDFmeas _).const_mul _ + refine MeasureTheory.Measure.ext_of_lintegral _ (fun f hf => ?_) + rw [compoundPoisson_lintegral _ lam hf] + rw [MeasureTheory.lintegral_add_measure, MeasureTheory.lintegral_smul_measure, + MeasureTheory.lintegral_dirac, + MeasureTheory.lintegral_withDensity_eq_lintegral_mul _ hd_meas hf] + rw [tsum_eq_zero_add' ENNReal.summable] + congr 1 + · rw [GammaConv.convPow_zero, MeasureTheory.lintegral_dirac, smul_eq_mul] + congr 1 + norm_num + · have hmeas_term : ∀ n : ℕ, Measurable (fun z => + ENNReal.ofReal (Real.exp (-(lam : ℝ)) * (lam : ℝ) ^ (n + 1) / ((n + 1).factorial)) + * gammaPDF (((n + 1 : ℕ) : ℝ) * α) γ z * f z) := by + intro n; exact ((hgPDFmeas _).const_mul _).mul hf + have step : ∀ n : ℕ, + ENNReal.ofReal (Real.exp (-(lam : ℝ)) * (lam : ℝ) ^ (n + 1) / ((n + 1).factorial)) + * ∫⁻ z, f z ∂(GammaConv.convPow (gammaMeasure α γ) (n + 1)) + = ∫⁻ z, ENNReal.ofReal (Real.exp (-(lam : ℝ)) * (lam : ℝ) ^ (n + 1) / ((n + 1).factorial)) + * gammaPDF (((n + 1 : ℕ) : ℝ) * α) γ z * f z ∂volume := by + intro n + rw [GammaConv.convPow_gamma α γ hα hγ (n + 1) (Nat.succ_pos n)] + rw [show gammaMeasure (((n + 1 : ℕ) : ℝ) * α) γ + = volume.withDensity (gammaPDF (((n + 1 : ℕ) : ℝ) * α) γ) from rfl] + rw [MeasureTheory.lintegral_withDensity_eq_lintegral_mul _ (hgPDFmeas _) hf] + simp only [Pi.mul_apply] + rw [← MeasureTheory.lintegral_const_mul _ ((hgPDFmeas _).mul hf)] + congr 1; ext z; ring + rw [tsum_congr step] + rw [← MeasureTheory.lintegral_tsum (fun n => (hmeas_term n).aemeasurable)] + congr 1; ext z + simp only [Pi.mul_apply, hd] + rw [← ENNReal.tsum_mul_right] + +open MeasureTheory ProbabilityTheory in +/-- For (Lebesgue-)almost every `y > 0`, the Tweedie term series `j ↦ twG μ φ p j y` is summable. +This follows from the summability of the integral norms (`tw_summable_norm`): the integrable +series has finite total integral, hence is summable pointwise a.e. -/ +lemma twG_ae_summable (μ : ℝ) {φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) : + ∀ᵐ y ∂volume, y ∈ Set.Ioi (0 : ℝ) → Summable (fun j : ℕ => twG μ φ p j y) := by + rw [← MeasureTheory.ae_restrict_iff' measurableSet_Ioi] + exact TweedieAux.ae_summable_of_summable_integral_norm + (fun j => tw_integrable_on hp₁ hp₂ hμ hφ j) + (tw_summable_norm μ φ p hp₁ hp₂ hμ hφ) + +open MeasureTheory ProbabilityTheory in +/-- The Poisson-weighted sum of the gamma densities (shape `(n+1)·α`) equals the Tweedie density, +Lebesgue-almost everywhere. -/ +lemma tweedie_density_series_ae (μ : ℝ) (hμ : 0 < μ) {φ p : ℝ} (hφ : 0 < φ) (hp₁ : 1 < p) + (hp₂ : p < 2) : + (fun z => ∑' n : ℕ, ENNReal.ofReal + (Real.exp (-twZ μ φ p) * (twZ μ φ p) ^ (n + 1) / ((n + 1).factorial)) + * gammaPDF (((n + 1 : ℕ) : ℝ) * ((2 - p) / (p - 1))) (μ ^ (1 - p) / (φ * (p - 1))) z) + =ᵐ[volume] (fun z => tweediePDF' μ hp₁ (le_of_lt hp₂) (le_of_lt hφ) z) := by + have hzne : ∀ᵐ z ∂(volume : Measure ℝ), z ≠ (0 : ℝ) := by + rw [MeasureTheory.ae_iff]; simp + filter_upwards [twG_ae_summable μ hp₁ hp₂ hμ hφ, hzne] with z hzsum hzne0 + have hpdf' : tweediePDF' μ hp₁ (le_of_lt hp₂) (le_of_lt hφ) z + = ENNReal.ofReal (tweediePDF μ φ p z) := by + unfold tweediePDF' + exact (ENNReal.ofReal_eq_coe_nnreal + (tweediePDF_nonneg (le_of_lt hφ) hp₁ (le_of_lt hp₂))).symm + rw [hpdf'] + rcases lt_or_gt_of_ne hzne0 with hzlt | hzgt + · -- z < 0 : both sides vanish + have h1 : tweediePDF μ φ p z = 0 := by + simp only [tweediePDF] + rw [Set.indicator_of_notMem (by simp only [Set.mem_setOf_eq]; linarith)] + rw [h1, ENNReal.ofReal_zero, ENNReal.tsum_eq_zero] + intro n + rw [gammaPDF_of_neg hzlt, mul_zero] + · -- z > 0 : use the pointwise Tweedie series + have hpw : tweediePDF μ φ p z = ∑' j : ℕ, twG μ φ p j z := by + simp only [tweediePDF] + rw [Set.indicator_of_mem (by simp only [Set.mem_setOf_eq]; exact hzgt)] + exact tw_pointwise μ φ p z + rw [hpw, + ENNReal.ofReal_tsum_of_nonneg + (fun j => twG_nonneg μ φ p hp₁ hp₂ hφ j hzgt) (hzsum hzgt)] + conv_rhs => rw [tsum_eq_zero_add' ENNReal.summable] + rw [twG_zero, ENNReal.ofReal_zero, zero_add] + have htwz : 0 ≤ twZ μ φ p := by rw [twZ]; positivity + refine tsum_congr (fun n => ?_) + rw [twG_eq_gamma μ hμ hp₁ hp₂ hφ (n + 1) hzgt, + ENNReal.ofReal_mul + (div_nonneg (mul_nonneg (Real.exp_nonneg _) (pow_nonneg htwz _)) (by positivity)), + gammaPDF] + +theorem tweedie_eq (μ : ℝ) (hμ : 0 < μ) {φ p : ℝ} + (hφ : 0 < φ) (hp₁ : 1 < p) + (hp₂ : p < 2) : + let hμ₀ := le_of_lt hμ + let hφ₀ := le_of_lt hφ + tweedieConstruction μ hμ₀ hφ₀ hp₂ = + tweedieMeasure μ hφ₀ hp₁ hp₂ := by + intro hμ₀ hφ₀ + have hα : (0 : ℝ) < (2 - p) / (p - 1) := div_pos (by linarith) (by linarith) + have hγ : (0 : ℝ) < μ ^ (1 - p) / (φ * (p - 1)) := + div_pos (Real.rpow_pos_of_pos hμ _) (by nlinarith) + unfold tweedieConstruction + rw [compound_split ((2 - p) / (p - 1)) (μ ^ (1 - p) / (φ * (p - 1))) hα hγ + ⟨μ ^ (2 - p) / (φ * (2 - p)), by positivity⟩] + rw [tweedieMeasure] + congr 1 + · -- atom at 0 + rw [ENNReal.smul_def] + congr 1 + rw [tweedieProbZero, ENNReal.ofReal_eq_coe_nnreal (Real.exp_nonneg _)] + congr 1 + apply NNReal.coe_injective + change Real.exp (-(μ ^ (2 - p) / (φ * (2 - p)))) + = Real.exp (-μ ^ (2 - p) / (φ * (2 - p))) + rw [neg_div] + · -- absolutely continuous part + exact MeasureTheory.withDensity_congr_ae (tweedie_density_series_ae μ hμ hφ hp₁ hp₂) diff --git a/Statlib/Tweedie/CompoundPoissonCore.lean b/Statlib/Tweedie/CompoundPoissonCore.lean new file mode 100644 index 0000000..c3a2ad0 --- /dev/null +++ b/Statlib/Tweedie/CompoundPoissonCore.lean @@ -0,0 +1,74 @@ +/- +Copyright (c) 2026 Bjørn Kjos-Hanssen. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Bjørn Kjos-Hanssen +-/ +import Mathlib +import Statlib.Tweedie.GammaConvolution + +/-! +# Compound Poisson construction (core, Mathlib-only) + +This file isolates the construction of the compound-Poisson law and basic facts about it. It +deliberately depends only on `Mathlib` (and `GammaConvolution`), so that it does not import the +Tweedie development. +-/ +open scoped BigOperators +open scoped Real +open scoped Pointwise +open scoped NNReal ENNReal +set_option maxRecDepth 4000 +set_option synthInstance.maxSize 128 +set_option relaxedAutoImplicit false +set_option autoImplicit false +set_option grind.warning false + + +namespace CompoundPoisson +open MeasureTheory ProbabilityTheory +variable (μ : Measure ℝ) [IsProbabilityMeasure μ] (lam : ℝ≥0) +/-- The base probability space: the first coordinate is the Poisson count `N`, the second is the +i.i.d. sequence of jumps with common law `μ`. -/ +noncomputable def baseMeasure : Measure (ℕ × (ℕ → ℝ)) := + (poissonMeasure lam).prod (Measure.infinitePi (fun _ : ℕ => μ)) +instance : IsProbabilityMeasure (baseMeasure μ lam) := by + unfold baseMeasure; infer_instance +/-- The compound Poisson random variable: the sum of the first `N` jumps. -/ +def jumpSum (p : ℕ × (ℕ → ℝ)) : ℝ := ∑ i ∈ Finset.range p.1, p.2 i +lemma measurable_jumpSum : Measurable (jumpSum : ℕ × (ℕ → ℝ) → ℝ) := by + unfold jumpSum + apply measurable_from_prod_countable_right + intro n; simp only + exact Finset.measurable_sum _ (fun i _ => measurable_pi_apply i) + +/-- The compound Poisson law: the push-forward of the base measure under `jumpSum`. -/ +noncomputable def compoundPoisson : Measure ℝ := + (baseMeasure μ lam).map jumpSum + +instance : IsProbabilityMeasure (compoundPoisson μ lam) := by + unfold compoundPoisson + exact Measure.isProbabilityMeasure_map measurable_jumpSum.aemeasurable +end CompoundPoisson + +open CompoundPoisson GammaConv +open MeasureTheory ProbabilityTheory + +/- +Lintegral against a compound-Poisson law: a Poisson-weighted sum of the convolution powers of +the jump law. +-/ +lemma compoundPoisson_lintegral (G : Measure ℝ) [IsProbabilityMeasure G] (lam : ℝ≥0) + {f : ℝ → ℝ≥0∞} (hf : Measurable f) : + ∫⁻ z, f z ∂(CompoundPoisson.compoundPoisson G lam) + = ∑' n : ℕ, ENNReal.ofReal (Real.exp (-(lam : ℝ)) * (lam : ℝ) ^ n / n.factorial) + * ∫⁻ z, f z ∂(GammaConv.convPow G n) := by + rw [compoundPoisson, MeasureTheory.lintegral_map hf CompoundPoisson.measurable_jumpSum, + baseMeasure, + MeasureTheory.lintegral_prod (fun a => f (jumpSum a)) (hf.comp measurable_jumpSum).aemeasurable, + MeasureTheory.lintegral_countable'] + refine tsum_congr fun n => ?_ + rw [poissonMeasure_singleton, mul_comm] + congr 1 + rw [← infinitePi_range_sum_eq_convPow, + MeasureTheory.lintegral_map hf (Finset.measurable_sum _ fun i _ => measurable_pi_apply i)] + rfl diff --git a/Statlib/Tweedie/GammaConvolution.lean b/Statlib/Tweedie/GammaConvolution.lean new file mode 100644 index 0000000..0392c54 --- /dev/null +++ b/Statlib/Tweedie/GammaConvolution.lean @@ -0,0 +1,314 @@ +/- +Copyright (c) 2026 Bjørn Kjos-Hanssen. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Bjørn Kjos-Hanssen +-/ +import Mathlib + +/-! +# Convolution of Gamma distributions + +The key analytic fact behind the compound-Poisson construction of the Tweedie +distribution: the convolution of two Gamma measures with the **same rate** `r` +is again a Gamma measure, with the shape parameters added: +`Gamma(a, r) ∗ Gamma(b, r) = Gamma(a + b, r)`. + +This is the statement that the sum of two independent Gamma random variables with +the same rate is Gamma distributed. +-/ + +open MeasureTheory ProbabilityTheory Measure Real +open scoped ENNReal NNReal + +namespace GammaConv + +/- +The real Beta integral over `[0,1]`, expressed via the Gamma function. +-/ +lemma real_betaIntegral_eq (a b : ℝ) (ha : 0 < a) (hb : 0 < b) : + ∫ t in (0:ℝ)..1, t ^ (a - 1) * (1 - t) ^ (b - 1) + = Real.Gamma a * Real.Gamma b / Real.Gamma (a + b) := by + rw [mul_comm, eq_div_iff] + · have h_beta : (∫ t in (0 : ℝ)..1, t^(a - 1) * (1 - t)^(b - 1)) + = Complex.betaIntegral (a : ℂ) (b : ℂ) := by + convert intervalIntegral.integral_ofReal.symm using 1 + convert intervalIntegral.integral_congr fun x hx => ?_ using 1 + norm_num [Complex.ofReal_cpow, show 0 ≤ x by aesop, show x ≤ 1 by aesop] + have := @Complex.Gamma_mul_Gamma_eq_betaIntegral (a : ℂ) (b : ℂ) ?_ ?_ <;> norm_cast at * + rw [← Complex.ofReal_inj] + simp_all only [Complex.ofReal_add, mul_comm, Complex.ofReal_mul] + rw [← Complex.Gamma_ofReal, ← Complex.Gamma_ofReal, ← Complex.Gamma_ofReal]; norm_cast at * + aesop + · positivity + +/-- Real Beta integral: `∫_0^z x^(a-1) (z-x)^(b-1) dx = z^(a+b-1) * Γ a * Γ b / Γ (a+b)` +for `a, b > 0` and `z > 0`. Proved by the substitution `x = z t`. -/ +lemma integral_rpow_mul_rpow_sub (a b : ℝ) (ha : 0 < a) (hb : 0 < b) {z : ℝ} (hz : 0 < z) : + ∫ x in (0:ℝ)..z, x ^ (a - 1) * (z - x) ^ (b - 1) + = z ^ (a + b - 1) * (Real.Gamma a * Real.Gamma b / Real.Gamma (a + b)) := by + have h := intervalIntegral.integral_comp_mul_left (a := 0) (b := 1) (c := z) + (fun x => x ^ (a - 1) * (z - x) ^ (b - 1)) hz.ne' + simp only [mul_zero, mul_one] at h + have h_subst : ∫ x in (0:ℝ)..z, x ^ (a - 1) * (z - x) ^ (b - 1) + = ∫ t in (0:ℝ)..1, (z * t) ^ (a - 1) * (z - z * t) ^ (b - 1) * z := by + have hg : ∫ t in (0:ℝ)..1, (z * t) ^ (a - 1) * (z - z * t) ^ (b - 1) * z + = z * ∫ t in (0:ℝ)..1, (z * t) ^ (a - 1) * (z - z * t) ^ (b - 1) := by + rw [← intervalIntegral.integral_const_mul]; apply intervalIntegral.integral_congr + intro x hx; ring + rw [hg, h, smul_eq_mul, ← mul_assoc, mul_inv_cancel₀ hz.ne', one_mul] + have h_simplify : ∫ t in (0:ℝ)..1, (z * t) ^ (a - 1) * (z - z * t) ^ (b - 1) * z + = z^(a - 1) * z^(b - 1) * z * ∫ t in (0:ℝ)..1, t ^ (a - 1) * (1 - t) ^ (b - 1) := by + rw [← intervalIntegral.integral_const_mul] + apply intervalIntegral.integral_congr + intro t ht + rw [Set.uIcc_of_le (by norm_num)] at ht + obtain ⟨ht0, ht1⟩ := ht + simp only + rw [Real.mul_rpow hz.le ht0, show z - z * t = z * (1 - t) by ring, + Real.mul_rpow hz.le (by linarith)] + ring + rw [h_subst, h_simplify, real_betaIntegral_eq a b ha hb] + rw [show z ^ (a - 1) * z ^ (b - 1) * z = z ^ (a + b - 1) by + rw [← Real.rpow_add hz, ← Real.rpow_add_one hz.ne']; ring_nf] + +/- +Convolution identity for gamma densities (real, pointwise) for `z > 0`. +-/ +lemma gamma_density_conv (a b r : ℝ) (ha : 0 < a) (hb : 0 < b) (hr : 0 < r) {z : ℝ} (hz : 0 < z) : + ∫ x in (0:ℝ)..z, gammaPDFReal a r x * gammaPDFReal b r (z - x) + = gammaPDFReal (a + b) r z := by + rw [intervalIntegral.integral_of_le hz.le] + rw [MeasureTheory.integral_Ioc_eq_integral_Ioo] + convert congr_arg (fun x : ℝ => (r ^ a / Real.Gamma a) * (r ^ b / Real.Gamma b) + * Real.exp (- (r * z)) * x) (integral_rpow_mul_rpow_sub a b ha hb hz) using 1 + · rw [intervalIntegral.integral_of_le hz.le, MeasureTheory.integral_Ioc_eq_integral_Ioo] + rw [← MeasureTheory.integral_const_mul] + refine MeasureTheory.setIntegral_congr_fun measurableSet_Ioo fun x hx => ?_ + rw [gammaPDFReal, gammaPDFReal] + rw [if_pos hx.1.le, if_pos (by linarith [hx.2])] + rw [show - (r * z) = - (r * x) + (r * x - r * z) by ring]; rw [Real.exp_add] + ring_nf + · unfold gammaPDFReal; rw [if_pos hz.le] + ring_nf + rw [Real.rpow_add hr] + norm_num [ne_of_gt (Real.Gamma_pos_of_pos ha), ne_of_gt (Real.Gamma_pos_of_pos hb)]; ring + +/- +The pointwise convolution of the gamma densities, at the level of `lintegral` over all of +`ℝ` w.r.t. Lebesgue measure. +-/ +lemma gammaPDF_conv_inner (a b r : ℝ) (ha : 0 < a) (hb : 0 < b) (hr : 0 < r) {z : ℝ} (hz : 0 < z) : + ∫⁻ x, gammaPDF a r x * gammaPDF b r (z - x) ∂volume = gammaPDF (a + b) r z := by + unfold gammaPDF + norm_num [← ENNReal.ofReal_mul, gammaPDFReal] + -- Apply the convolution identity for gamma densities. + have h_conv : ∫ x in Set.Ioo 0 z, (r ^ a / Real.Gamma a * x ^ (a - 1) * Real.exp (-(r * x))) + * (r ^ b / Real.Gamma b * (z - x) ^ (b - 1) * Real.exp (-(r * (z - x)))) + = r ^ (a + b) / Real.Gamma (a + b) * z ^ (a + b - 1) * Real.exp (-(r * z)) := by + convert congr_arg (fun x : ℝ => x * (r ^ a / Real.Gamma a) * (r ^ b / Real.Gamma b) + * Real.exp (- (r * z))) (integral_rpow_mul_rpow_sub a b ha hb hz) using 1 + <;> ring_nf + · rw [← MeasureTheory.integral_Ioc_eq_integral_Ioo, ← intervalIntegral.integral_of_le hz.le] + norm_num [mul_assoc, mul_comm, mul_left_comm, ← Real.exp_add]; ring_nf + exact Or.inl <| Or.inl <| by rw [← intervalIntegral.integral_const_mul]; congr; ext; ring_nf + · rw [Real.rpow_add hr] + ring_nf + norm_num [mul_assoc, mul_comm, mul_left_comm, ne_of_gt (Real.Gamma_pos_of_pos ha), + ne_of_gt (Real.Gamma_pos_of_pos hb)] + rw [← h_conv, if_pos hz.le, MeasureTheory.ofReal_integral_eq_lintegral_ofReal] + · rw [← MeasureTheory.lintegral_indicator] <;> norm_num [Set.indicator] + rw [← MeasureTheory.lintegral_congr_ae] + filter_upwards [ + MeasureTheory.measure_eq_zero_iff_ae_notMem.1 (MeasureTheory.measure_singleton 0), + MeasureTheory.measure_eq_zero_iff_ae_notMem.1 (MeasureTheory.measure_singleton z)] + with x hx₁ hx₂ + split_ifs + · rw [← ENNReal.ofReal_mul (mul_nonneg (mul_nonneg (div_nonneg (Real.rpow_nonneg hr.le _) + (Real.Gamma_nonneg_of_nonneg ha.le)) (Real.rpow_nonneg (by linarith) _)) + (Real.exp_nonneg _))] + · linarith + · linarith + · grind + · grind + · simp_all + · simp_all + · simp_all + · exact (by contrapose! h_conv; rw [MeasureTheory.integral_undef h_conv]; positivity) + · filter_upwards [MeasureTheory.ae_restrict_mem measurableSet_Ioo] with x hx using mul_nonneg + (mul_nonneg (mul_nonneg (by positivity) (Real.rpow_nonneg hx.1.le _)) (Real.exp_nonneg _)) + (mul_nonneg (mul_nonneg (by positivity) (Real.rpow_nonneg (sub_nonneg.2 hx.2.le) _)) + (Real.exp_nonneg _)) + +/- +The convolution of two gamma measures with the same rate `r` is a gamma measure with +the shapes added. +-/ +theorem gammaMeasure_conv (a b r : ℝ) (ha : 0 < a) (hb : 0 < b) (hr : 0 < r) : + (gammaMeasure a r) ∗ (gammaMeasure b r) = gammaMeasure (a + b) r := by + -- By definition of gamma measure, we know that + have h_gamma_def : ∀ c r : ℝ, 0 < c → 0 < r → gammaMeasure c r + = volume.withDensity (fun x => gammaPDF c r x) := by + aesop + -- Apply the equality of measures to rewrite the goal in terms of densities. + suffices h_conv : ∀ {g : ℝ → ENNReal}, Measurable g → ∫⁻ z, g z ∂(volume.withDensity + (fun x => gammaPDF a r x) ∗ volume.withDensity (fun x => gammaPDF b r x)) + = ∫⁻ z, g z * gammaPDF (a + b) r z ∂volume by + ext s hs; specialize @h_conv (Set.indicator s 1) + simp_all only [Set.indicator, Pi.one_apply, ite_mul, one_mul, zero_mul] + convert h_conv (measurable_const.indicator hs) using 1 + · erw [MeasureTheory.lintegral_indicator hs]; aesop + · rw [h_gamma_def _ _ (add_pos ha hb) hr, MeasureTheory.withDensity_apply _ hs] + rw [← MeasureTheory.lintegral_indicator] <;> norm_num [Set.indicator, hs] + intro g hg + have h_conv : ∫⁻ z, g z ∂(volume.withDensity (fun x => gammaPDF a r x) ∗ volume.withDensity + (fun x => gammaPDF b r x)) + = ∫⁻ x, ∫⁻ y, g (x + y) * gammaPDF b r y * gammaPDF a r x ∂volume ∂volume := by + rw [MeasureTheory.Measure.lintegral_conv] + · rw [MeasureTheory.lintegral_withDensity_eq_lintegral_mul] + · simp only [Pi.mul_apply, mul_comm, mul_left_comm] + congr! 1 + ext x + rw [MeasureTheory.lintegral_withDensity_eq_lintegral_mul] + · simp only [Pi.mul_apply, mul_comm] + · rw [← MeasureTheory.lintegral_const_mul'] + · congr; ext; ring + · exact ENNReal.ofReal_ne_top + · exact Measurable.ennreal_ofReal (measurable_gammaPDFReal b r) + · exact hg.comp (measurable_const.add measurable_id') + · exact Measurable.ennreal_ofReal (ProbabilityTheory.measurable_gammaPDFReal a r) + · refine Measurable.lintegral_prod_right ?_ + exact hg.comp (measurable_fst.add measurable_snd) + · exact hg + -- By Fubini's theorem, we can interchange the order of integration. + have h_fubini : ∫⁻ x, ∫⁻ y, g (x + y) * gammaPDF b r y * gammaPDF a r x ∂volume ∂volume + = ∫⁻ y, ∫⁻ x, g y * gammaPDF b r (y - x) * gammaPDF a r x ∂volume ∂volume := by + have h_fubini : ∀ x, ∫⁻ y, g (x + y) * gammaPDF b r y * gammaPDF a r x ∂volume + = ∫⁻ y, g y * gammaPDF b r (y - x) * gammaPDF a r x ∂volume := by + intro x; rw [← MeasureTheory.lintegral_sub_right_eq_self _ x]; congr; ext y; ring_nf + rw [funext h_fubini, MeasureTheory.lintegral_lintegral_swap] + refine AEMeasurable.mul ?_ ?_ + · refine AEMeasurable.mul ?_ ?_ + · exact hg.aemeasurable.comp_aemeasurable (measurable_snd.aemeasurable) + · refine Measurable.aemeasurable ?_ + refine Measurable.ennreal_ofReal ?_ + fun_prop + · refine Measurable.aemeasurable ?_ + exact (measurable_gammaPDFReal a r |> Measurable.comp <| measurable_fst).ennreal_ofReal + -- By the properties of the gamma function, we know that + -- $\int_{0}^{y} \gamma(a, r, x) \gamma(b, r, y-x) \, dx = \gamma(a+b, r, y)$ for $y > 0$. + have h_gamma_conv : ∀ y > 0, ∫⁻ x, gammaPDF a r x * gammaPDF b r (y - x) ∂volume + = gammaPDF (a + b) r y := by + intro y hy; exact (by + convert gammaPDF_conv_inner a b r ha hb hr hy using 1) + rw [h_conv, h_fubini] + rw [MeasureTheory.lintegral_congr_ae] + filter_upwards [ + MeasureTheory.measure_eq_zero_iff_ae_notMem.mp (MeasureTheory.measure_singleton 0)] with y hy + by_cases hy_pos : 0 < y + · rw [← h_gamma_conv y hy_pos, ← MeasureTheory.lintegral_const_mul] + · congr + ext + ring + · exact Measurable.mul (Measurable.ennreal_ofReal (measurable_gammaPDFReal _ _)) + (Measurable.ennreal_ofReal (measurable_gammaPDFReal _ _ |> Measurable.comp + <| measurable_const.sub measurable_id')) + · rw [MeasureTheory.lintegral_congr_ae, MeasureTheory.lintegral_zero] + · simp [gammaPDF_of_neg (show y < 0 from lt_of_le_of_ne (le_of_not_gt hy_pos) hy)] + · filter_upwards [] with x + by_cases hx : 0 < x + · simp_all only [gt_iff_lt, Set.mem_singleton_iff, not_lt, mul_eq_zero] + exact Or.inl <| Or.inr <| by rw [gammaPDF_of_neg (by linarith)] + · simp_all;grind +suggestions + +/-! ## Convolution powers and the law of a sum of i.i.d. variables -/ + +/-- The `n`-fold additive convolution power of a measure `G` (with `convPow G 0 = δ₀`). -/ +noncomputable def convPow (G : Measure ℝ) : ℕ → Measure ℝ + | 0 => Measure.dirac 0 + | (n + 1) => convPow G n ∗ G + +@[simp] lemma convPow_zero (G : Measure ℝ) : convPow G 0 = Measure.dirac 0 := rfl + +lemma convPow_succ (G : Measure ℝ) (n : ℕ) : convPow G (n + 1) = convPow G n ∗ G := rfl + +instance convPow_isProbabilityMeasure (G : Measure ℝ) [IsProbabilityMeasure G] (n : ℕ) : + IsProbabilityMeasure (convPow G n) := by + induction n with + | zero => rw [convPow_zero]; infer_instance + | succ n ih => rw [convPow_succ]; infer_instance + +/- +The `n`-fold convolution power of a gamma measure (for `n ≥ 1`) is a gamma measure with the +shape scaled by `n`. +-/ +lemma convPow_gamma (α γ : ℝ) (hα : 0 < α) (hγ : 0 < γ) (n : ℕ) (hn : 0 < n) : + convPow (gammaMeasure α γ) n = gammaMeasure ((n : ℝ) * α) γ := by + haveI : IsProbabilityMeasure (gammaMeasure α γ) := isProbabilityMeasure_gammaMeasure hα hγ + obtain ⟨m, rfl⟩ := Nat.exists_eq_succ_of_ne_zero hn.ne' + clear hn + induction m with + | zero => + rw [convPow_succ, convPow_zero, Measure.dirac_zero_conv] + norm_num + | succ k ih => + rw [convPow_succ, ih, gammaMeasure_conv _ _ _ (by positivity) hα hγ] + congr 1; push_cast; ring + +/- +The law of the sum of `n` i.i.d. variables (over `Fin n`) is the `n`-fold convolution power. +-/ +lemma finPiSum_eq_convPow (G : Measure ℝ) [IsProbabilityMeasure G] (n : ℕ) : + (Measure.pi (fun _ : Fin n => G)).map (fun v => ∑ i, v i) = convPow G n := by + induction n with + | zero => ext s hs; simp [hs] + | succ n ih => + have h_eq : (Measure.pi (fun _ : Fin (n + 1) => G)).map (fun v => ∑ i : Fin (n + 1), v i) + = (G.prod (Measure.pi (fun _ : Fin n => G))).map (fun p : ℝ × (Fin n → ℝ) + => p.1 + ∑ i : Fin n, p.2 i) := by + rw [← MeasureTheory.measurePreserving_piFinSuccAbove (fun _ : Fin (n + 1) => G) 0 |>.map_eq] + rw [MeasureTheory.Measure.map_map] + · congr with v; simp [Fin.sum_univ_succ] + rfl + · fun_prop + · fun_prop + rw [h_eq, convPow_succ] + have hsum : (fun p : ℝ × (Fin n → ℝ) => p.1 + ∑ i : Fin n, p.2 i) + = (fun q : ℝ × ℝ => q.1 + q.2) ∘ + (Prod.map (id : ℝ → ℝ) (fun w : Fin n → ℝ => ∑ i, w i)) := rfl + rw [hsum, ← MeasureTheory.Measure.map_map (by fun_prop) (by fun_prop), + ← MeasureTheory.Measure.map_prod_map _ _ (by fun_prop) (by fun_prop), + MeasureTheory.Measure.map_id, ih] + exact MeasureTheory.Measure.conv_comm G (convPow G n) + +/- +The law of the partial sum of the first `n` coordinates of an i.i.d. sequence (the infinite +product measure) is the `n`-fold convolution power. +-/ +lemma infinitePi_range_sum_eq_convPow (G : Measure ℝ) [IsProbabilityMeasure G] (n : ℕ) : + (Measure.infinitePi (fun _ : ℕ => G)).map (fun y => ∑ i ∈ Finset.range n, y i) + = convPow G n := by + convert finPiSum_eq_convPow G n using 1 + -- By definition of product measure, we can rewrite the right-hand side. + have h_prod : (Measure.pi fun _ : Fin n => G) = Measure.map + (fun y : ℕ → ℝ => fun i : Fin n => y i) (infinitePi fun _ => G) := by + refine MeasureTheory.Measure.pi_eq ?_ + intro s hs; erw [MeasureTheory.Measure.map_apply] + · convert MeasureTheory.Measure.infinitePi_pi _ _ + any_goals exact Finset.range n + rotate_left + rotate_left + · grind + · use fun i => if hi : i < n then s ⟨ i, hi ⟩ else Set.univ + · aesop + · grind + · rw [Finset.prod_range] + grind + · exact measurable_pi_lambda _ fun _ => measurable_pi_apply _ + · exact MeasurableSet.univ_pi hs + rw [h_prod, MeasureTheory.Measure.map_map] + · congr! 1 + ext; simp [Finset.sum_range] + · exact Finset.measurable_sum _ fun _ _ => measurable_pi_apply _ + · exact measurable_pi_lambda _ fun _ => measurable_pi_apply _ + +end GammaConv diff --git a/Statlib/Tweedie/Tweedie.lean b/Statlib/Tweedie/Tweedie.lean new file mode 100644 index 0000000..526eedb --- /dev/null +++ b/Statlib/Tweedie/Tweedie.lean @@ -0,0 +1,1101 @@ +/- +Copyright (c) 2026 Bjørn Kjos-Hanssen. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Bjørn Kjos-Hanssen +-/ +module +public import Mathlib.Analysis.SpecialFunctions.Gaussian.GaussianIntegral +public import Mathlib.Data.Real.StarOrdered +public import Mathlib.MeasureTheory.Measure.ProbabilityMeasure +public import Mathlib.Probability.Moments.Variance + +/-! +# Tweedie distribution: integral of the (compound-Poisson) density + + For `1 < p < 2` the Tweedie exponential dispersion model with mean `μ` and + dispersion `φ` has a density supported on `(0, ∞)` whose integral is + `1 - exp(-μ^(2-p) / (φ (2-p)))`; the + remaining mass is the atom at `0`. + +-/ + +@[expose] public section + +open MeasureTheory Set Real + +/-- The "a"-factor of the Tweedie density (the infinite series). -/ +noncomputable def a (y φ p : ℝ) := + let α := (2 - p) / (1 - p) + (1/y) * ∑' j : ℕ, (y ^ (- j * α) * (p - 1) ^(α * j)) / + (φ ^ (j * (1 - α)) * (2 - p) ^ j * Nat.factorial j * Gamma (- j * α)) + +/-- The Tweedie density. -/ +noncomputable def tweediePDF (μ φ p : ℝ) := + Set.indicator ({y | 0 < y}) + (fun y => a y φ p * exp ((1 / φ) * ((μ ^ (1 - p) / (1 - p)) * y - (μ ^ (2 - p) / (2 - p))))) + +/-- The Poisson rate `z = μ^(2-p) / (φ (2-p))`; the answer is `1 - exp(-z)`. -/ +noncomputable def twZ (μ φ p : ℝ) : ℝ := μ ^ (2 - p) / (φ * (2 - p)) + +/-- The `j`-th summand of the integrand `a y * exp(...)`, faithful to the definition of `a`. -/ +noncomputable def twG (μ φ p : ℝ) (j : ℕ) (y : ℝ) : ℝ := + let α := (2 - p) / (1 - p) + (1/y) * ((y ^ (- j * α) * (p - 1) ^(α * j)) / + (φ ^ (j * (1 - α)) * (2 - p) ^ j * Nat.factorial j * Gamma (- j * α))) + * exp ((1 / φ) * ((μ ^ (1 - p) / (1 - p)) * y - (μ ^ (2 - p) / (2 - p)))) + +/-- Closed form of `twG j y` for `y > 0`: a constant times `y^(-jα-1) * exp(-(rate)·y)`. -/ +lemma tw_pt (μ : ℝ) {φ p : ℝ} (hp₁ : 1 < p) (hφ : 0 < φ) (j : ℕ) {y : ℝ} (hy : 0 < y) : + twG μ φ p j y + = (Real.exp (-twZ μ φ p) * (p-1)^(((2-p)/(1-p))*(j:ℝ)) + / (φ^((j:ℝ)*(1-(2-p)/(1-p))) * (2-p)^j * (Nat.factorial j) + * Gamma (-(j:ℝ)*((2-p)/(1-p))))) + * (y ^ (-(j:ℝ)*((2-p)/(1-p)) - 1) * exp (-(μ ^ (1 - p) / (φ * (p - 1)) * y))) := by + have h1p : (1 - p) ≠ 0 := by linarith + have hpm1 : (p - 1) ≠ 0 := by linarith + rw [twG] + have hexp : exp ((1 / φ) * ((μ ^ (1 - p) / (1 - p)) * y - (μ ^ (2 - p) / (2 - p)))) + = exp (-twZ μ φ p) * exp (-(μ ^ (1 - p) / (φ * (p - 1)) * y)) := by + rw [← exp_add]; congr 1; rw [twZ]; field_simp; ring + rw [hexp] + have hypow : (1/y) * (y ^ (-(j:ℝ)*((2-p)/(1-p)))) = y ^ (-(j:ℝ)*((2-p)/(1-p)) - 1) := by + rw [Real.rpow_sub hy, rpow_one]; ring + rw [← hypow]; ring + +/-- The integrand is the pointwise tsum of the `twG` family. -/ +lemma tw_pointwise (μ φ p : ℝ) (y : ℝ) : + a y φ p * exp ((1 / φ) * ((μ ^ (1 - p) / (1 - p)) * y - (μ ^ (2 - p) / (2 - p)))) + = ∑' j : ℕ, twG μ φ p j y := by + simp only [a, twG] + rw [tsum_mul_right, tsum_mul_left] + +/-- The zeroth term vanishes identically (`Γ(0) = 0`). -/ +lemma twG_zero (μ φ p y : ℝ) : twG μ φ p 0 y = 0 := by + simp [twG, Gamma_zero] + +/-- The summand is nonnegative on `(0, ∞)`. -/ +lemma twG_nonneg (μ φ p : ℝ) (hp₁ : 1 < p) (hp₂ : p < 2) (hφ : 0 < φ) + (j : ℕ) {y : ℝ} (hy : 0 < y) : 0 ≤ twG μ φ p j y := by + cases j with + | zero => rw [twG_zero] + | succ n => + set j := n + 1 + have hjpos : j > 0 := by simp [j] + have hjposR : 0 < (j:ℝ) := by exact_mod_cast hjpos + have hαneg : (2 - p) / (1 - p) < 0 := div_neg_of_pos_of_neg (by linarith) (by linarith) + have ha0pos : 0 < -(j:ℝ) * ((2-p)/(1-p)) := by + have : -(j:ℝ) < 0 := by simpa using hjposR + exact mul_pos_of_neg_of_neg this hαneg + rw [tw_pt μ hp₁ hφ j hy] + positivity + + +/-- Each `twG j` is integrable on `(0, ∞)`. -/ +lemma tw_integrable_on {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) + (j : ℕ) : IntegrableOn (fun y => twG μ φ p j y) (Set.Ioi 0) := by + rcases Nat.eq_zero_or_pos j with hj0 | hjpos0 + · subst hj0; simp only [twG_zero]; exact integrableOn_zero + · have hjpos : 0 < (j:ℝ) := by exact_mod_cast hjpos0 + have hαneg : (2 - p) / (1 - p) < 0 := div_neg_of_pos_of_neg (by linarith) (by linarith) + have ha0pos : 0 < -(j:ℝ) * ((2-p)/(1-p)) := by + have : -(j:ℝ) < 0 := by simpa using hjpos + exact mul_pos_of_neg_of_neg this hαneg + have hrpos : 0 < μ ^ (1 - p) / (φ * (p - 1)) := by positivity + have hbase : IntegrableOn + (fun y => y ^ (-(j:ℝ)*((2-p)/(1-p)) - 1) + * exp (-(μ ^ (1 - p) / (φ * (p - 1))) * y ^ (1:ℝ))) (Set.Ioi 0) := + integrableOn_rpow_mul_exp_neg_mul_rpow (by linarith [ha0pos]) (le_refl 1) hrpos + have hbase2 : IntegrableOn + (fun y => y ^ (-(j:ℝ)*((2-p)/(1-p)) - 1) + * exp (-(μ ^ (1 - p) / (φ * (p - 1)) * y))) (Set.Ioi 0) := by + apply IntegrableOn.congr_fun hbase _ measurableSet_Ioi + intro y hy; simp only [Real.rpow_one, neg_mul] + have hconst : IntegrableOn + (fun y => (Real.exp (-twZ μ φ p) * (p-1)^(((2-p)/(1-p))*(j:ℝ)) + / (φ^((j:ℝ)*(1-(2-p)/(1-p))) * (2-p)^j * (Nat.factorial j) + * Gamma (-(j:ℝ)*((2-p)/(1-p))))) + * (y ^ (-(j:ℝ)*((2-p)/(1-p)) - 1) + * exp (-(μ ^ (1 - p) / (φ * (p - 1)) * y)))) (Set.Ioi 0) := + hbase2.const_mul _ + apply IntegrableOn.congr_fun hconst _ measurableSet_Ioi + intro y hy; simp only [Set.mem_Ioi] at hy + exact (tw_pt μ hp₁ hφ j hy).symm + +/-- The per-term integral for `j ≥ 1`: it equals `exp(-z) z^j / j!`. -/ +lemma tw_integral_term {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) + {j : ℕ} (hj : 1 ≤ j) : + ∫ y in Set.Ioi (0:ℝ), twG μ φ p j y + = exp (-twZ μ φ p) * (twZ μ φ p) ^ j / (Nat.factorial j) := by + have h1p : (1 - p) ≠ 0 := by linarith + have h2p : (2 - p) ≠ 0 := by linarith + have hp1 : 0 < p - 1 := by linarith + have hjpos : 0 < (j:ℝ) := by exact_mod_cast Nat.lt_of_lt_of_le Nat.zero_lt_one hj + have hαneg : (2 - p) / (1 - p) < 0 := div_neg_of_pos_of_neg (by linarith) (by linarith) + have ha0pos : 0 < -(j:ℝ) * ((2-p)/(1-p)) := by + have : -(j:ℝ) < 0 := by simpa using hjpos + exact mul_pos_of_neg_of_neg this hαneg + have hμpow : 0 < μ ^ (1 - p) := rpow_pos_of_pos hμ _ + have hrpos : 0 < μ ^ (1 - p) / (φ * (p - 1)) := by positivity + rw [setIntegral_congr_fun measurableSet_Ioi + (fun y hy => tw_pt μ hp₁ hφ j hy)] + rw [integral_const_mul] + rw [Real.integral_rpow_mul_exp_neg_mul_Ioi ha0pos hrpos] + set α := (2-p)/(1-p) with hα + set r := μ ^ (1 - p) / (φ * (p - 1)) with hrdef + have hrinvpos : 0 < 1/r := by positivity + have hW : (p-1)^α * (1/r)^(-α) / (φ^(1-α)*(2-p)) = μ^(2-p)/(φ*(2-p)) := by + have h1pneg : (1 - p) < 0 := by linarith + have hkey : (1 - p) * α = 2 - p := by rw [hα]; field_simp + have hrinv : (1/r) = φ * (p-1) / μ^(1-p) := by rw [hrdef, one_div_div] + rw [hrinv] + rw [Real.div_rpow (by positivity) (le_of_lt hμpow)] + rw [Real.mul_rpow (le_of_lt hφ) (le_of_lt hp1)] + rw [← rpow_mul hμ.le] + rw [show (1 - p) * (-α) = -(2-p) by rw [mul_neg, hkey]] + rw [Real.rpow_neg hμ.le (2-p), rpow_neg hp1.le α, rpow_neg hφ.le α] + have hφfac : φ^(1-α) = φ^(1:ℝ) * φ^(-α) := by rw [← rpow_add hφ]; ring_nf + rw [hφfac, rpow_one, rpow_neg hφ.le α] + field_simp + have hcore : (p-1)^(α*(j:ℝ)) * (1/r)^(-(j:ℝ)*α) / (φ^((j:ℝ)*(1-α))*(2-p)^j) + = (μ^(2-p)/(φ*(2-p)))^j := by + have hn1 : 0 ≤ (p-1)^α := rpow_nonneg hp1.le α + have hn2 : 0 ≤ (1/r)^(-α) := rpow_nonneg hrinvpos.le _ + have hn3 : 0 ≤ φ^(1-α) := rpow_nonneg hφ.le _ + have hn4 : 0 ≤ 2-p := by linarith + rw [Real.rpow_mul hp1.le α (j:ℝ)] + rw [show (-(j:ℝ)*α) = (-α)*(j:ℝ) by ring] + rw [Real.rpow_mul hrinvpos.le (-α) (j:ℝ)] + rw [show (j:ℝ)*(1-α) = (1-α)*(j:ℝ) by ring] + rw [Real.rpow_mul hφ.le (1-α) (j:ℝ)] + rw [← rpow_natCast (2-p) j, ← rpow_natCast (μ^(2-p)/(φ*(2-p))) j] + rw [← mul_rpow hn1 hn2] + rw [← mul_rpow hn3 hn4] + rw [← div_rpow (mul_nonneg hn1 hn2) (mul_nonneg hn3 hn4)] + rw [hW] + have hΓ : Gamma (-(j:ℝ)*α) ≠ 0 := ne_of_gt (Real.Gamma_pos_of_pos ha0pos) + have hfac : (Nat.factorial j : ℝ) ≠ 0 := by exact_mod_cast Nat.factorial_ne_zero j + rw [twZ] + rw [show exp (-(μ^(2-p)/(φ*(2-p)))) * (p-1)^(α*(j:ℝ)) + / (φ^((j:ℝ)*(1-α)) * (2-p)^j * (Nat.factorial j) * Gamma (-(j:ℝ)*α)) + * ((1/r)^(-(j:ℝ)*α) * Gamma (-(j:ℝ)*α)) + = exp (-(μ^(2-p)/(φ*(2-p)))) / (Nat.factorial j) + * (Real.Gamma (-(j:ℝ)*α)/Real.Gamma (-(j:ℝ)*α)) + * ((p-1)^(α*(j:ℝ)) * (1/r)^(-(j:ℝ)*α) / (φ^((j:ℝ)*(1-α))*(2-p)^j)) + by field_simp] + rw [div_self hΓ, hcore] + ring + +/-- The per-term integral for `j = 0` is `0`. -/ +lemma tw_integral_zero (μ φ p : ℝ) : + ∫ y in Set.Ioi (0:ℝ), twG μ φ p 0 y = 0 := by + simp [twG_zero] + +/-- Summability of the integral norms (needed to swap `∫` and `∑'`). -/ +lemma tw_summable_norm (μ φ p : ℝ) (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) : + Summable (fun j : ℕ => ∫ y in Set.Ioi (0:ℝ), ‖twG μ φ p j y‖) := by + have hsum2 : Summable (fun j : ℕ => + exp (-twZ μ φ p) * (twZ μ φ p) ^ j / (Nat.factorial j)) := by + have := (Real.summable_pow_div_factorial (twZ μ φ p)).mul_left (Real.exp (-twZ μ φ p)) + simpa [mul_div_assoc] using this + have hnorm_eq (j) : ∫ y in Set.Ioi 0, ‖twG μ φ p j y‖ + = ∫ y in Set.Ioi 0, twG μ φ p j y := by + apply setIntegral_congr_fun measurableSet_Ioi + intro y hy; simp only [Set.mem_Ioi] at hy + exact norm_of_nonneg (twG_nonneg μ φ p hp₁ hp₂ hφ j hy) + apply Summable.of_nonneg_of_le _ _ hsum2 + · intro j; rw [hnorm_eq] + exact setIntegral_nonneg measurableSet_Ioi + fun y => twG_nonneg μ φ p hp₁ hp₂ hφ j + · intro j + rw [hnorm_eq] + rcases Nat.eq_zero_or_pos j with hj0 | hjpos + · subst hj0; rw [tw_integral_zero]; positivity + · rw [tw_integral_term hp₁ hp₂ hμ hφ hjpos] + +/-- The series of per-term values sums to `1 - exp(-z)`. -/ +lemma tw_tsum (μ φ p : ℝ) : + ∑' j : ℕ, (if j = 0 then (0:ℝ) + else exp (-twZ μ φ p) * (twZ μ φ p) ^ j / (Nat.factorial j)) + = 1 - exp (-twZ μ φ p) := by + set z := twZ μ φ p + have hexp : exp z = ∑' n : ℕ, z ^ n / n.factorial := by + rw [Real.exp_eq_exp_ℝ, NormedSpace.exp_eq_tsum_div] + have hsum : Summable (fun n : ℕ => z ^ n / n.factorial) := summable_pow_div_factorial z + have hsum2 : Summable (fun n : ℕ => exp (-z) * z ^ n / (Nat.factorial n)) := by + have := hsum.mul_left (Real.exp (-z)) + simpa [mul_div_assoc] using this + have hsum3 : Summable (fun j : ℕ => + (if j = 0 then exp (-z) * z ^ j / (Nat.factorial j) else 0)) := by + apply summable_of_ne_finset_zero (s := {0}) + intro b hb; simp at hb; simp [hb] + have key : (fun j : ℕ => (if j = 0 then (0:ℝ) else exp (-z) * z ^ j / (Nat.factorial j))) + = (fun j => exp (-z) * z ^ j / (Nat.factorial j) + - (if j = 0 then exp (-z) * z ^ j / (Nat.factorial j) else 0)) := by + funext j; by_cases hj : j = 0 <;> simp [hj] + rw [key, Summable.tsum_sub hsum2 hsum3] + have h1 : ∑' j : ℕ, exp (-z) * z ^ j / (Nat.factorial j) = 1 := by + rw [show (fun j : ℕ => exp (-z) * z ^ j / (Nat.factorial j)) + = (fun j => exp (-z) * (z ^ j / (Nat.factorial j))) by funext j; ring] + rw [tsum_mul_left, ← hexp, ← exp_add]; simp + have h2 : ∑' j : ℕ, (if j = 0 then exp (-z) * z ^ j / (Nat.factorial j) else 0) + = exp (-z) := by + rw [tsum_eq_single 0] + · simp + · intro b hb; simp [hb] + rw [h1, h2] + +/-- The integral of the Tweedie density equals `1 - exp(-μ^(2-p)/(φ(2-p)))`. -/ +lemma tweediePDF_integral {μ φ p} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) + (hφ : 0 < φ) : + ∫ y, tweediePDF μ φ p y = + 1 - exp (-μ ^ (2 - p) / (φ * (2 - p))) := by + rw [tweediePDF] + rw [show {y : ℝ | 0 < y} = Set.Ioi 0 from rfl] + rw [MeasureTheory.integral_indicator measurableSet_Ioi] + rw [setIntegral_congr_fun measurableSet_Ioi (fun y _ => tw_pointwise μ φ p y)] + rw [← integral_tsum_of_summable_integral_norm + (fun j => tw_integrable_on hp₁ hp₂ hμ hφ j) + (tw_summable_norm μ φ p hp₁ hp₂ hμ hφ)] + rw [show (∑' j : ℕ, ∫ y in Set.Ioi (0:ℝ), twG μ φ p j y) + = ∑' j : ℕ, (if j = 0 then (0:ℝ) + else exp (-twZ μ φ p) * (twZ μ φ p) ^ j / (Nat.factorial j)) by + apply tsum_congr + intro j + rcases Nat.eq_zero_or_pos j with hj0 | hjpos + · subst hj0; rw [tw_integral_zero]; simp + · rw [tw_integral_term hp₁ hp₂ hμ hφ hjpos, if_neg (by omega)]] + rw [tw_tsum, twZ, neg_div] + +open NNReal ENNReal Real +noncomputable section + +/-- Probability of zero according to Tweedie distribution. -/ +def tweedieProbZero (μ φ p : ℝ) : ℝ≥0 := + ⟨rexp (-μ ^ (2 - p) / (φ * (2 - p))), exp_nonneg _⟩ + +/-- Nonnegativity of the Tweedie PDF. -/ +lemma tweediePDF_nonneg {y μ φ p : ℝ} (hφ : 0 ≤ φ) + (hp₁ : 1 < p) (hp₂ : p ≤ 2) + : tweediePDF μ φ p y ≥ 0 := by + simp only [tweediePDF, indicator, mem_setOf_eq, a, one_div, neg_mul, ge_iff_le] + split_ifs with g₀ + · apply mul_nonneg + · apply mul_nonneg + · positivity + · refine tsum_nonneg ?_ + intro j + apply mul_nonneg + · positivity + · simp only [mul_inv_rev] + apply mul_nonneg + · rw [inv_nonneg] + refine Gamma_nonneg_of_nonneg ?_ + rw [mul_div, ← neg_div, ← neg_mul, mul_comm, ← mul_div] + apply mul_nonneg + · linarith + · suffices 0 ≤ j / (p - 1) by + convert this using 1 + have : p - 1 = -(1 - p) := by simp + rw [this] + field_simp + positivity + · positivity + · exact exp_nonneg _ + simp + +/-- `ENNReal`-valued Tweedie probability density function. -/ +def tweediePDF' (μ : ℝ) {φ p : ℝ} + (hp₁ : 1 < p) (hp₂ : p ≤ 2) + (hφ : 0 ≤ φ) (y : ℝ) : ℝ≥0∞:= + let nn : NNReal := (⟨tweediePDF μ φ p y, by + by_cases H : y < 0 + · unfold tweediePDF + simp only [indicator, mem_setOf_eq, one_div] + rw [if_neg (by linarith)] + · simp only [not_lt] at H + exact tweediePDF_nonneg hφ hp₁ hp₂⟩ : NNReal) + (nn : ENNReal) + +/-- The Tweedie measure. -/ +def tweedieMeasure (μ : ℝ) {φ p : ℝ} (hφ : 0 ≤ φ) + (hp₁ : 1 < p) (hp₂ : p < 2) : Measure ℝ := + (tweedieProbZero μ φ p) • (Measure.dirac 0) + + (volume.withDensity (tweediePDF' μ hp₁ (by linarith) hφ)) + +/-- The Tweedie measure is a probability measure. -/ +lemma tweedieMeasure_prob (μ : ℝ) (hμ : 0 < μ) {φ p : ℝ} (hφ : 0 < φ) + (hp₁ : 1 < p) (hp₂ : p < 2) : + IsProbabilityMeasure (tweedieMeasure μ (show 0 ≤ φ by linarith) hp₁ hp₂) := by + refine isProbabilityMeasure_iff.mpr ?_ + unfold tweedieMeasure + simp only [Measure.coe_add, Measure.coe_smul, Pi.add_apply, Pi.smul_apply, measure_univ, + ENNReal.smul_one, MeasurableSet.univ, withDensity_apply, Measure.restrict_univ] + unfold tweedieProbZero + have : ∫⁻ (a : ℝ), tweediePDF' μ hp₁ (by linarith) (show 0 ≤ φ by linarith) a = + 1 - tweedieProbZero μ φ p := by + have ht := tweediePDF_integral hp₁ hp₂ hμ hφ + have : tweedieProbZero μ φ p ≤ 1 := by + refine exp_le_one_iff.mpr ?_ + field_simp + simp only [mul_zero, Left.neg_nonpos_iff] + positivity + have : (1 : NNReal) - ⟨rexp (-μ ^ (2 - p) / (φ * (2 - p))), exp_nonneg _⟩ + = ⟨1 - tweedieProbZero μ φ p, by + generalize tweedieProbZero μ φ p = α at * + exact sub_nonneg_of_le this⟩ := by + unfold tweedieProbZero + have (a : ℝ) (ha : a < 1) (ha' : 0 ≤ a) : + (1 : NNReal) - ⟨a,ha'⟩ = ⟨1-a, by linarith⟩ := by + refine (toNNReal_eq_iff_eq_coe ?_).mpr rfl + have : 1 - a > 0 := by linarith + exact Ne.symm (Std.ne_of_lt this) + apply this + refine exp_lt_one_iff.mpr ?_ + ring_nf + simp only [Left.neg_neg_iff] + apply _root_.mul_pos <| rpow_pos_of_pos hμ (2 - p) + · simp only [inv_pos, lt_neg_add_iff_add_lt, add_zero] + nth_rw 1 [mul_comm] + exact (mul_lt_mul_iff_of_pos_left hφ).mpr hp₂ + have : (1 : ENNReal) - ENNReal.ofNNReal ⟨rexp (-μ ^ (2 - p) / (φ * (2 - p))), exp_nonneg _⟩ + = ENNReal.ofNNReal ⟨1 - tweedieProbZero μ φ p, + sub_nonneg_of_le (by simp;tauto)⟩ := by rw [← this]; simp + unfold tweedieProbZero + rw [this] + simp only [tweedieProbZero] + have : rexp (-μ ^ (2 - p) / (φ * (2 - p))) = 1 - ∫ (y : ℝ), tweediePDF μ φ p y := by linarith + simp_rw [this] + have (a : Real) (ha : 0 ≤ 1 - a) + (ha' : 0 ≤ a) + : ofNNReal (⟨(1 : ℝ) - (⟨1 - a, ha⟩ : NNReal), + by simp;tauto⟩ : NNReal) = ofNNReal ⟨a, ha'⟩ := by simp + specialize this (∫ (y : ℝ), tweediePDF μ φ p y) (by + rw [tweediePDF_integral] + all_goals try linarith + simp only [_root_.sub_sub_cancel] + exact exp_nonneg (-μ ^ (2 - p) / (φ * (2 - p)))) + (by + rw [tweediePDF_integral] + all_goals try linarith + simp + field_simp + simp only [mul_zero, Left.neg_nonpos_iff] + apply mul_nonneg + · refine rpow_nonneg ?_ (2 - p) + linarith + · simp + linarith) + symm + unfold tweediePDF' + convert this + · rw [MeasureTheory.lintegral_coe_eq_integral] + · refine (toReal_eq_toReal_iff' ?_ ?_).mp ?_ + · simp + · simp + · refine toReal_ofReal ?_ + refine integral_nonneg ?_ + intro + simp + · suffices Integrable (fun x ↦ tweediePDF μ φ p x) volume by + convert this + apply Integrable.of_integral_ne_zero + rw [tweediePDF_integral hp₁ hp₂ hμ hφ] + have (a) (ha : a ≠ 0) : 1 - rexp a ≠ 0 := by + contrapose! ha + apply exp_injective + rw [exp_zero] + linarith + apply this + simp only [ne_eq, _root_.div_eq_zero_iff, neg_eq_zero, mul_eq_zero, not_or] + constructor + · refine (rpow_ne_zero ?_ ?_).mpr ?_ + all_goals linarith + · constructor + all_goals linarith + rw [this] + unfold tweedieProbZero + have (a : NNReal) (ha : 1 - a.1 > 0) : + a + (1 - a) = 1 := by + apply NNReal.coe_injective + change a.toReal + (1-a).toReal = 1 + have : 1 - a.toReal = (1-a).toReal := by + refine (toNNReal_eq_iff_eq_coe ?_).mp rfl + apply ne_of_gt + refine NNReal.coe_pos.mp ?_ + simp only [val_eq_coe, gt_iff_lt, sub_pos, coe_lt_one, NNReal.coe_pos, + tsub_pos_iff_lt] at ha ⊢ + exact ha + rw [← this] + linarith + have (a : NNReal) (ha : 1 - a.1 > 0) : + ofNNReal a + ofNNReal (1 - a) = 1 := by + rw [← ENNReal.coe_add] + rw [this] + · simp + · tauto + apply this + simp + field_simp + simp only [mul_zero, Left.neg_neg_iff] + apply div_pos + · apply rpow_pos_of_pos + linarith + · linarith + +/-- The Tweedie measure as a probability measure. -/ +def tweedieProbMeasure {μ φ p : ℝ} (hμ : 0 < μ) (hφ : 0 < φ) + (hp₁ : 1 < p) (hp₂ : p < 2) : MeasureTheory.ProbabilityMeasure ℝ := { + val := tweedieMeasure μ (by linarith) hp₁ hp₂ + property := + tweedieMeasure_prob μ hμ hφ hp₁ hp₂ + } + + +/-! ## Expectation of the Tweedie distribution + +We now prove that the mean (expectation) of `tweedieProbMeasure` equals the parameter `μ`. +The measure is a mixture of an atom at `0` (which contributes nothing to the mean) and an +absolutely continuous part with density `tweediePDF`. The mean of the continuous part is +`∫ y, y * tweediePDF μ φ p y`, which we evaluate term-by-term. +-/ + +/- +Closed form of `y * twG j y` for `y > 0`: the same constant as in `tw_pt`, but with +the power `y^(-jα)` (the factor `y` cancels one power of `y`). +-/ +lemma tw_yG_pt (μ : ℝ) {φ p : ℝ} (hp₁ : 1 < p) (hφ : 0 < φ) (j : ℕ) {y : ℝ} (hy : 0 < y) : + y * twG μ φ p j y + = (Real.exp (-twZ μ φ p) * (p-1)^(((2-p)/(1-p))*(j:ℝ)) + / (φ^((j:ℝ)*(1-(2-p)/(1-p))) * (2-p)^j * (Nat.factorial j) + * Real.Gamma (-(j:ℝ)*((2-p)/(1-p))))) + * (y ^ (-(j:ℝ)*((2-p)/(1-p))) * Real.exp (-(μ ^ (1 - p) / (φ * (p - 1)) * y))) := by + convert congr_arg (fun x : ℝ => y * x) (tw_pt μ hp₁ hφ j hy) using 1; ring_nf + have : y ^ (p * (1 - p)⁻¹ * ↑j - (1 - p)⁻¹ * ↑j * 2) + = y * y ^ (-1 + p * (1 - p)⁻¹ * ↑j - (1 - p)⁻¹ * ↑j * 2) := by + have (z : ℝ) : y^z * y = y ^ (z + 1) := by + refine Eq.symm (Real.rpow_add_one ?_ z) + linarith + simp_rw [mul_comm] at this + rw [this] + ring_nf + rw [this] + ring_nf + +/- +Each `y * twG j` is integrable on `(0, ∞)`. +-/ +lemma tw_yG_integrable_on {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) + (j : ℕ) : IntegrableOn (fun y => y * twG μ φ p j y) (Set.Ioi 0) := by + rcases Nat.eq_zero_or_pos j with hj0 | hjpos0 + · simp [hj0, twG_zero] + · have h_integrable : IntegrableOn (fun y => y ^ (-(j:ℝ)*((2-p)/(1-p))) * Real.exp (-(μ ^ (1 - p) + / (φ * (p - 1)) * y))) (Set.Ioi 0) := by + have h_integrable : ∀ {s b : ℝ}, -1 < s → 0 < b → + IntegrableOn (fun y => y ^ s * Real.exp (-b * y)) (Set.Ioi 0) := by + intro s b hs hb + convert (integrableOn_rpow_mul_exp_neg_mul_rpow (show -1 < s by linarith) + (show 1 ≤ (1 : ℝ) by norm_num) hb) using 1; norm_num + convert h_integrable _ _ using 1 + rotate_left + · exact - (j : ℝ) * ((2 - p) / (1 - p)) + · exact μ ^ (1 - p) / (φ * (p - 1)) + · nlinarith [show (j : ℝ) ≥ 1 by + norm_cast, mul_div_cancel₀ (2 - p) (by linarith : (1 - p) ≠ 0)] + · exact div_pos (Real.rpow_pos_of_pos hμ _) (mul_pos hφ (by linarith)) + · norm_num + refine h_integrable.const_mul ?_ |> fun h => h.congr ?_ + · exact (Real.exp (-twZ μ φ p) * (p - 1) ^ (((2 - p) / (1 - p)) * j) / (φ ^ ((j : ℝ) * + (1 - (2 - p) / (1 - p))) * (2 - p) ^ j * (Nat.factorial j) * Real.Gamma (- (j : ℝ) + * ((2 - p) / (1 - p))))) + filter_upwards [MeasureTheory.ae_restrict_mem measurableSet_Ioi] with y hy + using by rw [tw_yG_pt μ hp₁ hφ j hy] + +/- +The per-term mean integral: `∫ y * twG j = K * (j · exp(-z) z^j / j!)`, +where `K = φ(2-p)/μ^(1-p)`. (Valid for all `j`, since both sides vanish at `j = 0`.) +-/ +lemma tw_mean_term {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) (j : ℕ) : + ∫ y in Set.Ioi (0:ℝ), y * twG μ φ p j y + = (φ * (2-p) / μ^(1-p)) + * ((j:ℝ) * Real.exp (-twZ μ φ p) * (twZ μ φ p) ^ j / (Nat.factorial j)) := by + by_cases hj : j = 0 + · simp_all [twG] + simp_all only [twG, one_div, neg_mul] + have h_integral : ∫ y in Set.Ioi (0:ℝ), y * twG μ φ p j y + = ∫ y in Set.Ioi (0:ℝ), (Real.exp (-twZ μ φ p) * (p-1)^(((2-p)/(1-p))*(j:ℝ)) + / (φ^((j:ℝ)*(1-(2-p)/(1-p))) * (2-p)^j * (Nat.factorial j) * Real.Gamma (-(j:ℝ)*((2-p)/(1-p))))) + * (y ^ (-(j:ℝ)*((2-p)/(1-p))) * Real.exp (-(μ ^ (1 - p) / (φ * (p - 1)) * y))) := by + exact MeasureTheory.setIntegral_congr_fun measurableSet_Ioi fun y hy => tw_yG_pt μ hp₁ hφ j hy + -- Apply the integral formula for the gamma function. + have h_gamma : ∫ y in Set.Ioi (0:ℝ), y ^ (-(j:ℝ)*((2-p)/(1-p))) * Real.exp (-(μ ^ (1 - p) + / (φ * (p - 1)) * y)) = (1 / (μ ^ (1 - p) / (φ * (p - 1)))) ^ (-(j:ℝ)*((2-p)/(1-p)) + 1) + * Real.Gamma (-(j:ℝ)*((2-p)/(1-p)) + 1) := by + convert integral_rpow_mul_exp_neg_mul_Ioi _ _ using 1 + · norm_num + · nlinarith [show (j : ℝ) ≥ 1 by exact Nat.one_le_cast.mpr (Nat.pos_of_ne_zero hj), + mul_div_cancel₀ (2 - p) (by linarith : (1 - p) ≠ 0)] + · exact div_pos (Real.rpow_pos_of_pos hμ _) (mul_pos hφ (by linarith)) + convert h_integral using 1 + · unfold twG; congr; ext + ring_nf + · rw [MeasureTheory.integral_const_mul, h_gamma] + rw [Real.Gamma_add_one] + · rw [Real.rpow_add, Real.rpow_one] + · field_simp + rw [eq_comm, neg_div', div_eq_iff] + · unfold twZ; ring_nf + rw [show (-1 + p) = (p - 1) by ring, + show (p * φ * (μ ^ (1 - p)) ⁻¹ - φ * (μ ^ (1 - p)) ⁻¹) + = (p - 1) * φ * (μ ^ (1 - p)) ⁻¹ by ring] + rw [Real.mul_rpow (by nlinarith) (by positivity), + Real.mul_rpow (by nlinarith) (by positivity)]; ring_nf + rw [show (- (p * j * (1 - p) ⁻¹) + j * (1 - p) ⁻¹ * 2 : ℝ) + = - (p * j * (1 - p) ⁻¹ - j * (1 - p) ⁻¹ * 2) by ring, Real.rpow_neg (by linarith)] + ring_nf + rw [show (-1 + p) = (p - 1) by ring, Real.rpow_def_of_pos (by linarith), + Real.rpow_def_of_pos (by positivity), Real.rpow_def_of_pos (by positivity)]; ring_nf + rw [Real.log_inv, Real.log_rpow (by positivity)]; ring_nf + rw [Real.rpow_def_of_pos (by positivity)]; ring_nf + rw [Real.rpow_def_of_pos (by positivity)]; ring_nf + rw [show (- (p * φ) + φ * 2) = φ * (2 - p) by ring, mul_inv] + norm_num [Real.exp_add, Real.exp_neg, Real.exp_nat_mul, Real.exp_log, hμ, hφ, hp₁, hp₂] + ring_nf + norm_num [← Real.exp_nat_mul, ← Real.exp_neg, ← Real.exp_add, ← Real.exp_sub, + hμ.ne', hφ.ne', show (2 - p) ≠ 0 by linarith]; ring_nf + norm_num [mul_assoc, ← Real.exp_add]; ring_nf + grind +revert + · exact mul_ne_zero (mul_ne_zero (by linarith) (pow_ne_zero _ (by linarith))) (ne_of_gt + (Real.Gamma_pos_of_pos ( + by nlinarith [ + show (j : ℝ) ≥ 1 by exact Nat.one_le_cast.mpr (Nat.pos_of_ne_zero hj), + mul_div_cancel₀ ((2 - p) * j) (by linarith : (1 - p) ≠ 0)]))) + · exact one_div_pos.mpr (div_pos (Real.rpow_pos_of_pos hμ _) (mul_pos hφ (by linarith))) + · exact mul_ne_zero (neg_ne_zero.mpr (Nat.cast_ne_zero.mpr hj)) (div_ne_zero (by linarith) + (by linarith)) + +/-- Summability of `fun j => (j:ℝ) * z^j / j!` for any real `z`. -/ +lemma tw_summable_n_pow (z : ℝ) : + Summable (fun j : ℕ => (j:ℝ) * z ^ j / (Nat.factorial j)) := by + refine summable_of_ratio_norm_eventually_le ?_ ?_ (r := 2 / 3) + · norm_num + · norm_num [Nat.factorial_succ, abs_div, abs_mul] + refine ⟨⌈3 * |z|⌉₊ + 1, fun n hn => ?_⟩; rw [mul_div_mul_left _ _ (by positivity)]; ring_nf + nlinarith [show (n : ℝ) ≥ ⌈3 * |z|⌉₊ + 1 by exact_mod_cast hn, Nat.le_ceil (3 * |z|), + show (0 : ℝ) ≤ |z| ^ n * (n.factorial : ℝ) ⁻¹ by positivity] + +/-- The "Poisson mean" identity: `∑' j, j z^j / j! = z · exp z`. -/ +lemma tw_tsum_n_pow (z : ℝ) : + ∑' j : ℕ, (j:ℝ) * z ^ j / (Nat.factorial j) = z * Real.exp z := by + convert (Summable.tsum_eq_zero_add + (show Summable fun j : ℕ => ((j : ℝ) * z ^ j / (j.factorial : ℝ)) from ?_)) using 1 + · norm_num [Nat.factorial_succ, pow_succ', mul_assoc, mul_comm, mul_left_comm, div_eq_mul_inv, + tsum_mul_left] + field_simp + rw [Real.exp_eq_exp_ℝ, NormedSpace.exp_eq_tsum_div] + simp +decide only [mul_div_assoc] + exact Eq.symm _root_.tsum_mul_left + · exact tw_summable_n_pow z + +/- +Summability of the integral norms for the mean (needed to swap `∫` and `∑'`). +-/ +lemma tw_mean_summable_norm (μ φ p : ℝ) (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) : + Summable (fun j : ℕ => ∫ y in Set.Ioi (0:ℝ), ‖y * twG μ φ p j y‖) := by + convert Summable.mul_left (φ * (2 - p) / μ ^ (1 - p) * Real.exp (-twZ μ φ p)) + (tw_summable_n_pow (twZ μ φ p)) using 2 + rw [MeasureTheory.setIntegral_congr_fun measurableSet_Ioi fun y hy => + Real.norm_of_nonneg (mul_nonneg hy.out.le (twG_nonneg μ φ p hp₁ hp₂ hφ _ hy))]; + rw [tw_mean_term hp₁ hp₂ hμ hφ] + ring + +/- +The series of per-term mean values sums to `μ`. +-/ +lemma tw_mean_tsum {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) : + ∑' j : ℕ, (φ * (2-p) / μ^(1-p)) + * ((j:ℝ) * Real.exp (-twZ μ φ p) * (twZ μ φ p) ^ j / (Nat.factorial j)) = μ := by + -- Let `z = twZ μ φ p` and `K = φ*(2-p)/μ^(1-p)`. + -- Each summand is `K * ((j:ℝ) * exp(-z) * z^j / j!)`. + set z := twZ μ φ p + set K := φ * (2 - p) / μ ^ (1 - p) with hK + -- So the sum equals `K * exp(-z) * (z * exp z) = K * z * (exp (-z) * exp z) = K * z` + -- (since `exp (-z) * exp z = 1`, via `Real.exp_neg` and `inv_mul_cancel` or `← Real.exp_add`). + have hsum : ∑' j : ℕ, (j : ℝ) * Real.exp (-z) * z ^ j / (Nat.factorial j) + = z * Real.exp z * Real.exp (-z) := by + have := @tw_tsum_n_pow z + simp_all +decide only [div_eq_mul_inv, mul_comm, mul_left_comm, mul_assoc] + simp +decide only [← mul_assoc, ← this] + exact _root_.tsum_mul_right + convert congr_arg (fun x : ℝ => K * x) hsum using 1 + · exact _root_.tsum_mul_left + · simp +zetaDelta only at * + unfold twZ; norm_num [mul_assoc, ← Real.exp_add]; ring_nf + rw [show (2 - p) = (1 - p) + 1 by ring, Real.rpow_add hμ, Real.rpow_one]; ring_nf + nlinarith [mul_inv_cancel_left₀ (show φ * 2 - φ * p ≠ 0 by nlinarith) μ, mul_inv_cancel₀ + (show μ ^ (1 - p) ≠ 0 by positivity)] + +/- +The mean of the continuous part: `∫ y, tweediePDF μ φ p y * y = μ`. +-/ +lemma tweedie_mean_value {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) : + ∫ y, tweediePDF μ φ p y * y = μ := by + -- First, show `tweediePDF μ φ p y * y = Set.indicator (Set.Ioi 0) + -- (fun y => (a y φ p * Real.exp (...)) * y) y`. + have h_indicator : (∫ y, tweediePDF μ φ p y * y) + = (∫ y in Set.Ioi (0:ℝ), (a y φ p * Real.exp ((1 / φ) * ((μ ^ (1 - p) / (1 - p)) + * y - (μ ^ (2 - p) / (2 - p)))) * y)) := by + rw [← MeasureTheory.integral_indicator] <;> norm_num [Set.indicator, tweediePDF] + rw [h_indicator] + -- On `Ioi 0`, `(a y φ p * Real.exp (...)) * y = y * (a y φ p * Real.exp (...)) + -- = y * ∑' j, twG μ φ p j y = ∑' j, y * twG μ φ p j y` + -- (use `tw_pointwise` then `tsum_mul_left`). Apply `setIntegral_congr_fun measurableSet_Ioi`. + have h_indicator : (∫ y in Set.Ioi (0:ℝ), (a y φ p * Real.exp ((1 / φ) * ((μ ^ (1 - p) / (1 - p)) + * y - (μ ^ (2 - p) / (2 - p)))) * y)) + = (∫ y in Set.Ioi (0:ℝ), ∑' j : ℕ, y * twG μ φ p j y) := by + refine MeasureTheory.setIntegral_congr_fun measurableSet_Ioi fun y hy => ?_ + convert congr_arg (fun x : ℝ => x * y) (tw_pointwise μ φ p y) using 1; ring_nf + exact _root_.tsum_mul_left + rw [h_indicator] + convert tw_mean_tsum hp₁ hp₂ hμ hφ using 1 + rw [← MeasureTheory.integral_tsum_of_summable_integral_norm] + · exact tsum_congr fun j => tw_mean_term hp₁ hp₂ hμ hφ j + · exact fun j => tw_yG_integrable_on hp₁ hp₂ hμ hφ j + · convert tw_mean_summable_norm μ φ p hp₁ hp₂ hμ hφ using 1 + +/-- `tweediePDF μ φ p y * y` is integrable. -/ +lemma tweedie_mean_integrable {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) : + Integrable (fun y => tweediePDF μ φ p y * y) := by + apply Integrable.of_integral_ne_zero + rw [tweedie_mean_value hp₁ hp₂ hμ hφ] + exact ne_of_gt hμ + +/- +The expectation of the Tweedie measure equals `μ`. +-/ +theorem tweedieMeasure_expectation {μ φ p : ℝ} (hμ : 0 < μ) (hφ : 0 < φ) + (hp₁ : 1 < p) (hp₂ : p < 2) : + ∫ y, y ∂(tweedieMeasure μ (show (0:ℝ) ≤ φ by linarith) hp₁ hp₂) = μ := by + -- Let `f y := (tweediePDF μ φ p y).toNNReal`. Then `f y • y = tweediePDF μ φ p y * y`. + set f : ℝ → ℝ≥0 := fun y => (tweediePDF μ φ p y).toNNReal + -- Show that `volume.withDensity (tweediePDF' μ hp₁ (by linarith) _)` + -- equals `volume.withDensity (fun y => ↑(f y))`. + have h_withDensity : volume.withDensity (tweediePDF' μ hp₁ (by linarith) hφ.le) + = volume.withDensity (fun y => f y) := by + refine MeasureTheory.withDensity_congr_ae ?_ + filter_upwards [] with y using (ENNReal.ofReal_eq_coe_nnreal _).symm + -- Show that `Integrable (fun y => y) (volume.withDensity (fun y => ↑(f y)))`. + have h_integrable : Integrable (fun y => y) (volume.withDensity (fun y => f y)) := by + have h_integrable : Integrable (fun y => f y • y) volume := by + convert tweedie_mean_integrable hp₁ hp₂ hμ hφ using 1 + ext y; exact mul_eq_mul_right_iff.mpr (Or.inl <| Real.coe_toNNReal _ <| + tweediePDF_nonneg (by linarith) hp₁ (by linarith)) + rw [MeasureTheory.integrable_withDensity_iff_integrable_smul₀] + · convert h_integrable using 1 + · have h_integrable : Integrable (fun y => tweediePDF μ φ p y) volume := by + contrapose! h_integrable + have := tweediePDF_integral hp₁ hp₂ hμ hφ + rw [MeasureTheory.integral_undef h_integrable] at this + linarith [Real.exp_pos (-μ ^ (2 - p) / (φ * (2 - p))), + Real.exp_lt_one_iff.mpr (show -μ ^ (2 - p) / (φ * (2 - p)) < 0 by + exact div_neg_of_neg_of_pos (neg_neg_of_pos (Real.rpow_pos_of_pos hμ _)) + (mul_pos hφ (by linarith)))] + exact h_integrable.1.aemeasurable.real_toNNReal + convert congr_arg (fun x : ℝ => x + 0) (tweedie_mean_value hp₁ hp₂ hμ hφ) using 1 + · rw [tweedieMeasure, MeasureTheory.integral_add_measure] + · norm_num [h_withDensity, h_integrable] + · convert integral_withDensity_eq_integral_smul₀ _ (fun y => y) using 1 + · simp +zetaDelta only + congr! 1 + ext y; simp +decide [Real.toNNReal_of_nonneg (tweediePDF_nonneg + (by linarith : 0 ≤ φ) hp₁ (by linarith))] + rfl + · have h_integrable : Integrable (fun y => tweediePDF μ φ p y) volume := by + exact (by + apply MeasureTheory.Integrable.of_integral_ne_zero + have := tweediePDF_integral hp₁ hp₂ hμ hφ + linarith [Real.exp_pos (-μ ^ (2 - p) / (φ * (2 - p))), + Real.exp_lt_one_iff.mpr (show -μ ^ (2 - p) / (φ * (2 - p)) < 0 by + exact div_neg_of_neg_of_pos (neg_neg_of_pos (Real.rpow_pos_of_pos hμ _)) + (mul_pos hφ (by linarith)))]) + exact h_integrable.1.aemeasurable.real_toNNReal + · constructor <;> norm_num [MeasureTheory.HasFiniteIntegral] + exact measurable_id.aestronglyMeasurable + · aesop + · ring + +/-- The expectation of the Tweedie probability measure equals `μ`. -/ +theorem tweedieProbMeasure_expectation {μ φ p : ℝ} (hμ : 0 < μ) (hφ : 0 < φ) + (hp₁ : 1 < p) (hp₂ : p < 2) : + ∫ y, y ∂ (tweedieProbMeasure hμ hφ hp₁ hp₂) = μ := + tweedieMeasure_expectation hμ hφ hp₁ hp₂ + + + +/-! ## Variance of the Tweedie distribution + +We prove the famous identity `Var = φ · μ^p`. The argument mirrors the mean +computation: we evaluate the second moment `∫ y², dP` term-by-term and obtain +`μ² + φ·μ^p`, and then the variance is `E[Y²] - (E[Y])² = φ·μ^p`. +-/ + +/- +Closed form of `y² * twG j y` for `y > 0`: the same constant as in `tw_pt`, but with +the power `y^(-jα+1)` (the factor `y²` cancels one power of `y` and adds one). +-/ +lemma tw_y2G_pt (μ : ℝ) {φ p : ℝ} (hp₁ : 1 < p) (hφ : 0 < φ) (j : ℕ) {y : ℝ} (hy : 0 < y) : + y ^ 2 * twG μ φ p j y + = (Real.exp (-twZ μ φ p) * (p-1)^(((2-p)/(1-p))*(j:ℝ)) + / (φ^((j:ℝ)*(1-(2-p)/(1-p))) * (2-p)^j * (Nat.factorial j) + * Real.Gamma (-(j:ℝ)*((2-p)/(1-p))))) + * (y ^ (-(j:ℝ)*((2-p)/(1-p)) + 1) * Real.exp (-(μ ^ (1 - p) / (φ * (p - 1)) * y))) := by + convert congr_arg (fun x : ℝ => y * x) (tw_yG_pt μ hp₁ hφ j hy) using 1 + · ring + · rw [ Real.rpow_add hy, Real.rpow_one ]; ring + +/- +Each `y² * twG j` is integrable on `(0, ∞)`. +-/ +lemma tw_y2G_integrable_on {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) + (j : ℕ) : IntegrableOn (fun y => y ^ 2 * twG μ φ p j y) (Set.Ioi 0) := by + by_cases hj : j = 0 + · unfold twG; aesop + · have h_integrable : IntegrableOn (fun y => y ^ (-(j:ℝ)*((2-p)/(1-p)) + 1) + * Real.exp (-(μ ^ (1 - p) / (φ * (p - 1)) * y))) (Set.Ioi 0) := by + have h_integrable : ∀ {s b : ℝ}, -1 < s → 0 < b → IntegrableOn (fun y => y ^ s + * Real.exp (-b * y)) (Set.Ioi 0) := by + intro s b hs hb + convert (integrableOn_rpow_mul_exp_neg_mul_rpow (show -1 < s by linarith) + (show 1 ≤ (1 : ℝ) by norm_num) hb) using 1; norm_num + convert h_integrable _ _ using 1 <;> norm_num [ neg_div, div_neg ] + · congr! 1 + · nlinarith [ show (j : ℝ) ≥ 1 by exact Nat.one_le_cast.mpr (Nat.pos_of_ne_zero hj), + mul_div_cancel₀ (2 - p) (by linarith : (1 - p) ≠ 0) ] + · exact div_pos (Real.rpow_pos_of_pos hμ _) (mul_pos hφ (by linarith)) + refine h_integrable.const_mul ?_ |> fun h => h.congr ?_ + · exact Real.exp (-twZ μ φ p) * (p - 1) ^ (((2 - p) / (1 - p)) * j) + / (φ ^ ((j : ℝ) * (1 - (2 - p) / (1 - p))) * (2 - p) ^ j * (j.factorial : ℝ) + * Real.Gamma (- (j : ℝ) * ((2 - p) / (1 - p)))) + filter_upwards [ MeasureTheory.ae_restrict_mem measurableSet_Ioi ] with y hy using by + rw [ tw_y2G_pt μ hp₁ hφ j hy ] + +/- +Recurrence relating the second-moment per-term integral to the mean per-term integral: +`∫ y², G_j = ((-jα+1)·(1/b)) · ∫ y, G_j`, where `b = μ^(1-p)/(φ(p-1))`, i.e. `1/b = φ(p-1)/μ^(1-p)`. +This follows from `∫ y^(s+1) e^{-by} = ((s+1)/b) ∫ y^s e^{-by}` (a Gamma recurrence). +-/ +lemma tw_2nd_moment_recurrence {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) + (j : ℕ) : + ∫ y in Set.Ioi (0:ℝ), y ^ 2 * twG μ φ p j y + = ((-(j:ℝ)*((2-p)/(1-p)) + 1) * (φ * (p-1) / μ^(1-p))) + * ∫ y in Set.Ioi (0:ℝ), y * twG μ φ p j y := by + by_cases hj : j = 0 + · unfold twG; aesop + · have h_integral : ∫ y in Set.Ioi (0:ℝ), y ^ 2 * twG μ φ p j y + = (∫ y in Set.Ioi (0:ℝ), y ^ (-(j:ℝ)*((2-p)/(1-p)) + 1) + * Real.exp (-(μ ^ (1 - p) / (φ * (p - 1)) * y))) * (Real.exp (-twZ μ φ p) + * (p-1)^(((2-p)/(1-p))*(j:ℝ)) / (φ^((j:ℝ)*(1-(2-p)/(1-p))) * (2-p)^j * (Nat.factorial j) + * Real.Gamma (-(j:ℝ)*((2-p)/(1-p))))) := by + rw [ ← MeasureTheory.integral_mul_const ] + refine MeasureTheory.setIntegral_congr_fun measurableSet_Ioi fun y hy => ?_ + rw [ tw_y2G_pt μ hp₁ hφ j hy ]; ring + have h_integral : ∫ y in Set.Ioi (0:ℝ), y * twG μ φ p j y + = (∫ y in Set.Ioi (0:ℝ), y ^ (-(j:ℝ)*((2-p)/(1-p))) * Real.exp (-(μ ^ (1 - p) + / (φ * (p - 1)) * y))) * (Real.exp (-twZ μ φ p) * (p-1)^(((2-p)/(1-p))*(j:ℝ)) / (φ^((j:ℝ) + * (1-(2-p)/(1-p))) * (2-p)^j * (Nat.factorial j) * Real.Gamma (-(j:ℝ)*((2-p)/(1-p))))) := by + rw [ ← MeasureTheory.integral_mul_const ] + refine MeasureTheory.setIntegral_congr_fun measurableSet_Ioi fun y hy => ?_ + rw [ tw_yG_pt μ hp₁ hφ j hy ]; ring + rw [ ‹∫ y in Set.Ioi 0, y ^ 2 * twG μ φ p j y = _›, h_integral ] + have h_integral : ∫ y in Set.Ioi (0:ℝ), y ^ (-(j:ℝ)*((2-p)/(1-p)) + 1) * Real.exp (-(μ ^ (1 - p) + / (φ * (p - 1)) * y)) = (1 / (μ ^ (1 - p) / (φ * (p - 1)))) ^ (-(j:ℝ)*((2-p)/(1-p)) + 2) + * Real.Gamma (-(j:ℝ)*((2-p)/(1-p)) + 2) := by + convert integral_rpow_mul_exp_neg_mul_Ioi _ _ using 1 + · exact MeasureTheory.setIntegral_congr_fun measurableSet_Ioi fun x hx => by ring_nf + · nlinarith [ show (j : ℝ) ≥ 1 by exact Nat.one_le_cast.mpr (Nat.pos_of_ne_zero hj), + mul_div_cancel₀ (2 - p) (by linarith : (1 - p) ≠ 0) ] + · exact div_pos (Real.rpow_pos_of_pos hμ _) (mul_pos hφ (by linarith)) + rw [ h_integral, show (∫ y in Set.Ioi 0, y ^ (- (j : ℝ) * ((2 - p) / (1 - p))) + * Real.exp (- (μ ^ (1 - p) / (φ * (p - 1)) * y))) = (1 / (μ ^ (1 - p) + / (φ * (p - 1)))) ^ (- (j : ℝ) * ((2 - p) / (1 - p)) + 1) * Real.Gamma (- (j : ℝ) * ((2 - p) + / (1 - p)) + 1) from ?_ ] + · rw [ show (-j * ((2 - p) / (1 - p)) + 2 : ℝ) + = (-j * ((2 - p) / (1 - p)) + 1) + 1 by ring, Real.rpow_add_one ] <;> norm_num + · rw [ show (- (j * ((2 - p) / (1 - p))) + 1 + 1 : ℝ) + = (- (j * ((2 - p) / (1 - p))) + 1) + 1 by ring, + Real.Gamma_add_one (by nlinarith [ + show (j : ℝ) ≥ 1 by exact Nat.one_le_cast.mpr (Nat.pos_of_ne_zero hj), + mul_div_cancel₀ (2 - p) (by linarith : (1 - p) ≠ 0) + ] + ) + ] + ring + · exact ⟨⟨hφ.ne', by linarith⟩, by positivity⟩ + · convert integral_rpow_mul_exp_neg_mul_Ioi _ _ using 1 + · norm_num + · nlinarith [ show (j : ℝ) ≥ 1 by exact Nat.one_le_cast.mpr (Nat.pos_of_ne_zero hj), + mul_div_cancel₀ (2 - p) (by linarith : (1 - p) ≠ 0) ] + · exact div_pos (Real.rpow_pos_of_pos hμ _) (mul_pos hφ (by linarith)) + +/-- The per-term second-moment integral. (Valid for all `j`, since both sides vanish at `j = 0`.) -/ +lemma tw_2nd_moment_term {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) (j : ℕ) : + ∫ y in Set.Ioi (0:ℝ), y ^ 2 * twG μ φ p j y + = (φ * (2-p) / μ^(1-p)) * (φ * (p-1) / μ^(1-p)) + * (((j:ℝ) + (2-p)/(p-1) * (j:ℝ)^2) * Real.exp (-twZ μ φ p) + * (twZ μ φ p) ^ j / (Nat.factorial j)) := by + rw [tw_2nd_moment_recurrence hp₁ hp₂ hμ hφ, tw_mean_term hp₁ hp₂ hμ hφ] + have h1p : (1 - p) ≠ 0 := by linarith + have hpm1 : (p - 1) ≠ 0 := by linarith + have hμp : μ ^ (1 - p) ≠ 0 := by positivity + field_simp + ring + +/- +Summability of `fun j => (j:ℝ)^2 * z^j / j!` for any real `z`. +-/ +lemma tw_summable_n2_pow (z : ℝ) : + Summable (fun j : ℕ => (j:ℝ)^2 * z ^ j / (Nat.factorial j)) := by + refine summable_of_ratio_norm_eventually_le ?_ ?_ (r := 2 / 3) + · norm_num + · -- For large enough $n$, the term $(n+1) * |z| / n^2$ will be less than $2/3$. + have h_bound : ∃ N : ℕ, ∀ n ≥ N, (n + 1) * |z| / n^2 ≤ 2 / 3 := by + exact ⟨⌈3 * |z|⌉₊ + 1, fun n hn => by rw [ + div_le_iff₀ ] <;> nlinarith [ Nat.le_ceil (3 * |z|), + show (n : ℝ) ≥ ⌈3 * |z|⌉₊ + 1 by exact_mod_cast hn, abs_nonneg z ]⟩ + obtain ⟨N, hN⟩ := h_bound + filter_upwards [ Filter.eventually_ge_atTop N, Filter.eventually_gt_atTop 0 ] with n hn hn' + specialize hN n hn + simp_all +decide only [Nat.cast_add, Nat.cast_one, Nat.factorial_succ, Nat.cast_mul, norm_div, + norm_mul, norm_pow, norm_eq_abs, sq_abs, RCLike.norm_natCast] + convert mul_le_mul_of_nonneg_right hN + (show 0 ≤ (n ^ 2 * |z| ^ n / n.factorial : ℝ) by positivity) using 1; norm_cast; ring_nf + -- Simplifying the right-hand side: + field_simp + ring_nf + push_cast; ring + +/- +The "Poisson second moment" identity: `∑' j, j² z^j / j! = (z²+z) · exp z`. +-/ +lemma tw_tsum_n2_pow (z : ℝ) : + ∑' j : ℕ, (j:ℝ)^2 * z ^ j / (Nat.factorial j) = (z^2 + z) * Real.exp z := by + -- We'll use the fact that $\sum_{j=0}^{\infty} j(j-1) \frac{z^j}{j!} = z^2 e^z$. + have h1 : ∑' j : ℕ, (j * (j - 1) : ℝ) * z ^ j / j.factorial = z^2 * Real.exp z := by + -- Split the sum into two parts: one for $j=0$ and $j=1$, and the rest. + have h_split : ∑' j : ℕ, (j * (j - 1) : ℝ) * z^j / j.factorial + = ∑' j : ℕ, if j ≥ 2 then (j * (j - 1) : ℝ) * z^j / j.factorial else 0 := by + exact tsum_congr fun n => by rcases n with (_ | _ | n) <;> norm_num + -- For $j \geq 2$, we can simplify the term + -- $(j * (j - 1) : ℝ) * z^j / j.factorial$ to $z^2 * z^{j-2} / (j-2)!$. + have h_simplify : ∀ j : ℕ, j ≥ 2 → (j * (j - 1) : ℝ) * z^j / j.factorial + = z^2 * z^(j-2) / (j-2).factorial := by + intro j hj; rcases j with (_ | _ | j) <;> norm_num [ Nat.factorial ] at * + rw [ div_eq_div_iff ] <;> first | positivity | ring! + -- Substitute the simplified terms back into the sum. + have h_sum_simplified : ∑' j : ℕ, (if j ≥ 2 then (j * (j - 1) : ℝ) * z^j / j.factorial else 0) + = ∑' j : ℕ, z^2 * z^j / j.factorial := by + rw [ ← tsum_eq_tsum_of_ne_zero_bij ] + · use fun x => x.val - 2 + · exact fun x y h => by + rcases x with ⟨_ | _ | x, hx⟩ <;> rcases y with ⟨_ | _ | y, hy⟩ <;> cases h <;> trivial + · intro x hx; use ⟨x + 2, by aesop⟩; aesop + · aesop + simp_all +decide only [ge_iff_le, exp_eq_exp_ℝ, NormedSpace.exp_eq_tsum_div] + simp +decide only [mul_div_assoc] + exact _root_.tsum_mul_left + convert congr_arg₂ (· + ·) h1 (tw_tsum_n_pow z) using 1 + · rw [ ← Summable.tsum_add ] + · congr; ext j; ring + · contrapose! h1 + rw [ tsum_eq_zero_of_not_summable h1 ]; norm_num [ Real.exp_ne_zero ] + exact fun h => h1 <| by subst h; exact ⟨_, hasSum_single 0 fun j hj => by aesop⟩ + · convert tw_summable_n_pow z using 1 + · ring + +/- +Summability of the second-moment integral norms (needed to swap `∫` and `∑'`). +-/ +lemma tw_2nd_moment_summable_norm (μ φ p : ℝ) + (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) : + Summable (fun j : ℕ => ∫ y in Set.Ioi (0:ℝ), ‖y ^ 2 * twG μ φ p j y‖) := by + have h_summable : Summable (fun j : ℕ => ∫ y in Set.Ioi (0:ℝ), y ^ 2 * twG μ φ p j y) := by + simp_all +decide only [tw_2nd_moment_term] + refine Summable.mul_left _ ?_ + convert Summable.add (Summable.mul_left (Real.exp (-twZ μ φ p)) + (tw_summable_n_pow (twZ μ φ p))) (Summable.mul_left (Real.exp (-twZ μ φ p) + * (2 - p) / (p - 1)) (tw_summable_n2_pow (twZ μ φ p))) using 2; ring + convert h_summable using 1 + generalize_proofs at * + exact funext fun j => MeasureTheory.setIntegral_congr_fun measurableSet_Ioi fun y hy => by + rw [ Real.norm_of_nonneg (mul_nonneg (sq_nonneg _) (twG_nonneg μ φ p hp₁ hp₂ hφ j hy)) ] + +/- +The series of per-term second-moment values sums to `μ² + φ·μ^p`. +-/ +lemma tw_2nd_moment_tsum {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) : + ∑' j : ℕ, (φ * (2-p) / μ^(1-p)) * (φ * (p-1) / μ^(1-p)) + * (((j:ℝ) + (2-p)/(p-1) * (j:ℝ)^2) * Real.exp (-twZ μ φ p) + * (twZ μ φ p) ^ j / (Nat.factorial j)) + = μ^2 + φ * μ^p := by + -- Factor out common terms and apply the series summations. + have h_series : ∑' j : ℕ, (j + (2 - p) / (p - 1) * j^2) * Real.exp (-twZ μ φ p) + * (twZ μ φ p) ^ j / j.factorial = + Real.exp (-twZ μ φ p) * (twZ μ φ p * Real.exp (twZ μ φ p) + (2 - p) / (p - 1) + * ((twZ μ φ p)^2 + twZ μ φ p) * Real.exp (twZ μ φ p)) := by + have h_series : ∑' j : ℕ, (j : ℝ) * twZ μ φ p ^ j / j.factorial = twZ μ φ p + * Real.exp (twZ μ φ p) ∧ ∑' j : ℕ, (j : ℝ) ^ 2 * twZ μ φ p ^ j / j.factorial + = (twZ μ φ p ^ 2 + twZ μ φ p) * Real.exp (twZ μ φ p) := by + exact ⟨by simpa using tw_tsum_n_pow (twZ μ φ p), + by simpa using tw_tsum_n2_pow (twZ μ φ p)⟩ + convert congr_arg₂ (· + ·) (congr_arg (fun x : ℝ => x * Real.exp (-twZ μ φ p)) h_series.1) + (congr_arg (fun x : ℝ => x * (2 - p) * Real.exp (-twZ μ φ p) / (p - 1)) h_series.2) using 1 + · norm_num [ add_mul, mul_assoc, mul_div_assoc, tsum_mul_left, tsum_mul_right ] + convert congr_arg₂ (· + ·) (tsum_mul_right) (tsum_mul_right) using 1 + · rw [ ← Summable.tsum_add ] + · congr; ext j; ring + · exact Summable.mul_right _ <| by simpa only [ mul_div_assoc ] using tw_summable_n_pow _ + · convert Summable.mul_right ((2 - p) * (Real.exp (-twZ μ φ p) / (p - 1))) + (tw_summable_n2_pow (twZ μ φ p)) using 2; ring + all_goals infer_instance + · ring + convert congr_arg (fun x : ℝ => (φ * (2 - p) / μ ^ (1 - p)) * (φ * (p - 1) / μ ^ (1 - p)) * x) + h_series using 1 + · exact _root_.tsum_mul_left + · unfold twZ; norm_num [ Real.rpow_sub hμ ]; ring_nf + norm_num [ Real.exp_neg, Real.exp_ne_zero ]; ring_nf + field_simp + grind + +/- +The second moment of the continuous part: `∫ y, tweediePDF μ φ p y * y² = μ² + φ·μ^p`. +-/ +lemma tweedie_2nd_moment_value {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) (hμ : 0 < μ) (hφ : 0 < φ) : + ∫ y, tweediePDF μ φ p y * y ^ 2 = μ^2 + φ * μ^p := by + -- Apply the decomposition of the integral into the sum of integrals over the components + -- and use the fact that the integral of `y^2` with respect to the point mass part is 0. + have h_decomp : ∫ y, tweediePDF μ φ p y * y ^ 2 + = ∫ y in Set.Ioi (0:ℝ), ∑' j : ℕ, y ^ 2 * twG μ φ p j y := by + rw [ ← MeasureTheory.integral_indicator ] <;> norm_num [ Set.indicator ] + congr with x + by_cases hx : 0 < x + · simp +decide only [tweediePDF, one_div, mem_setOf_eq, hx, indicator_of_mem, ↓reduceIte] + convert congr_arg (fun y : ℝ => x ^ 2 * y) (tw_pointwise μ φ p x) using 1 + · ring_nf + · exact _root_.tsum_mul_left + · simp +decide [ hx, tweediePDF ] + rw [ h_decomp, ← tw_2nd_moment_tsum hp₁ hp₂ hμ hφ, + ← MeasureTheory.integral_tsum_of_summable_integral_norm ] + · exact tsum_congr fun j => tw_2nd_moment_term hp₁ hp₂ hμ hφ j + · exact fun j => tw_y2G_integrable_on hp₁ hp₂ hμ hφ j + · convert tw_2nd_moment_summable_norm μ φ p hp₁ hp₂ hμ hφ using 1 + +/- +`tweediePDF μ φ p y * y²` is integrable. +-/ +lemma tweedie_2nd_moment_integrable {μ φ p : ℝ} (hp₁ : 1 < p) (hp₂ : p < 2) + (hμ : 0 < μ) (hφ : 0 < φ) : + Integrable (fun y => tweediePDF μ φ p y * y ^ 2) := by + have := tweedie_2nd_moment_value hp₁ hp₂ hμ hφ + exact (by contrapose! this; rw [ MeasureTheory.integral_undef this ]; positivity) + +/- +The second moment of the Tweedie measure equals `μ² + φ·μ^p`. +-/ +theorem tweedieMeasure_2nd_moment {μ φ p : ℝ} (hμ : 0 < μ) (hφ : 0 < φ) + (hp₁ : 1 < p) (hp₂ : p < 2) : + ∫ y, y ^ 2 ∂(tweedieMeasure μ (show (0:ℝ) ≤ φ by linarith) hp₁ hp₂) = μ^2 + φ * μ^p := by + set f : ℝ → ℝ≥0 := fun y => (tweediePDF μ φ p y).toNNReal + have h_withDensity : volume.withDensity (tweediePDF' μ hp₁ (by linarith) hφ.le) = + volume.withDensity (fun y => f y) := by + refine MeasureTheory.withDensity_congr_ae ?_ + filter_upwards [ ] with y using (ENNReal.ofReal_eq_coe_nnreal _).symm + convert congr_arg (fun x : ℝ => x + 0) (tweedie_2nd_moment_value hp₁ (by linarith) hμ hφ) using 1 + · generalize_proofs at * + · unfold tweedieMeasure; norm_num [ h_withDensity ] + rw [ MeasureTheory.integral_add_measure ] <;> norm_num [ tweedieProbZero ] + · convert integral_withDensity_eq_integral_smul₀ _ (fun y => y ^ 2) using 1 + · congr! 1 + generalize_proofs at * + ext; simp [f]; ring_nf + exact mul_eq_mul_right_iff.mpr (Or.inl <| by + rw [ Real.toNNReal_of_nonneg <| tweediePDF_nonneg (by linarith) hp₁ (by linarith) ] + norm_cast) + · have h_integrable : MeasureTheory.Integrable (fun y => tweediePDF μ φ p y) volume := by + exact (by + apply MeasureTheory.Integrable.of_integral_ne_zero + have := tweediePDF_integral hp₁ (by linarith) hμ hφ + linarith [ Real.exp_pos (-μ ^ (2 - p) / (φ * (2 - p))), Real.exp_lt_one_iff.mpr + (show -μ ^ (2 - p) / (φ * (2 - p)) < 0 by + exact div_neg_of_neg_of_pos (neg_neg_of_pos (Real.rpow_pos_of_pos hμ _)) + (mul_pos hφ (by linarith))) ]) + generalize_proofs at * + exact h_integrable.1.aemeasurable.real_toNNReal + · constructor + · exact Continuous.aestronglyMeasurable (continuous_pow 2) + · simp +decide [ HasFiniteIntegral ] + · have h_integrable : Integrable (fun y => f y • y ^ 2) volume := by + convert tweedie_2nd_moment_integrable hp₁ hp₂ hμ hφ using 1 + generalize_proofs at * + ext y + exact mul_eq_mul_right_iff.mpr (Or.inl <| Real.coe_toNNReal _ <| + tweediePDF_nonneg (by linarith) hp₁ (by linarith)) + rw [ MeasureTheory.integrable_withDensity_iff_integrable_smul₀ ] + · exact h_integrable + · have h_integrable : AEMeasurable (fun y => tweediePDF μ φ p y) volume := by + have h_integrable : Integrable (fun y => tweediePDF μ φ p y) volume := by + exact (by + apply MeasureTheory.Integrable.of_integral_ne_zero + have := tweediePDF_integral hp₁ (by linarith) hμ hφ + linarith [ + Real.exp_pos (-μ ^ (2 - p) / (φ * (2 - p))), + Real.exp_lt_one_iff.mpr ( + show -μ ^ (2 - p) / (φ * (2 - p)) < 0 by + exact div_neg_of_neg_of_pos (neg_neg_of_pos (Real.rpow_pos_of_pos hμ _)) + (mul_pos hφ (by linarith))) + ] + ) + generalize_proofs at *; exact h_integrable.1.aemeasurable + generalize_proofs at * + exact AEMeasurable.real_toNNReal h_integrable + · ring + +/-- The second moment of the Tweedie probability measure equals `μ² + φ·μ^p`. -/ +theorem tweedieProbMeasure_2nd_moment {μ φ p : ℝ} (hμ : 0 < μ) (hφ : 0 < φ) + (hp₁ : 1 < p) (hp₂ : p < 2) : + ∫ y, y ^ 2 ∂ (tweedieProbMeasure hμ hφ hp₁ hp₂) = μ^2 + φ * μ^p := + tweedieMeasure_2nd_moment hμ hφ hp₁ hp₂ + +/-- **The variance of the Tweedie distribution is `φ · μ^p`.** +Here the variance is expressed as `E[Y²] - (E[Y])²`. -/ +theorem tweedieProbMeasure_variance {μ φ p : ℝ} (hμ : 0 < μ) (hφ : 0 < φ) + (hp₁ : 1 < p) (hp₂ : p < 2) : + (∫ y, y ^ 2 ∂ (tweedieProbMeasure hμ hφ hp₁ hp₂)) + - (∫ y, y ∂ (tweedieProbMeasure hμ hφ hp₁ hp₂)) ^ 2 = φ * μ^p := by + rw [tweedieProbMeasure_2nd_moment hμ hφ hp₁ hp₂, tweedieProbMeasure_expectation hμ hφ hp₁ hp₂] + ring + +/- +`y ↦ y²` is integrable with respect to the Tweedie measure +(the second moment is finite). +-/ +lemma tweedieMeasure_sq_integrable {μ φ p : ℝ} (hμ : 0 < μ) (hφ : 0 < φ) + (hp₁ : 1 < p) (hp₂ : p < 2) : + Integrable (fun y => y ^ 2) (tweedieMeasure μ (show (0:ℝ) ≤ φ by linarith) hp₁ hp₂) := by + apply MeasureTheory.Integrable.of_integral_ne_zero + have := tweedieMeasure_2nd_moment hμ hφ hp₁ hp₂ + nlinarith [ Real.rpow_pos_of_pos hμ p ] + +/-- **The variance of the Tweedie distribution is `φ · μ^p`**, stated with Mathlib's +`ProbabilityTheory.variance` (`Var[Y] = 𝔼[(Y - 𝔼[Y])²]`). -/ +theorem tweedieMeasure_variance {μ φ p : ℝ} (hμ : 0 < μ) (hφ : 0 < φ) + (hp₁ : 1 < p) (hp₂ : p < 2) : + ProbabilityTheory.variance (fun y => y) + (tweedieMeasure μ (show (0:ℝ) ≤ φ by linarith) hp₁ hp₂) = φ * μ^p := by + haveI : IsProbabilityMeasure (tweedieMeasure μ (show (0:ℝ) ≤ φ by linarith) hp₁ hp₂) := + tweedieMeasure_prob μ hμ hφ hp₁ hp₂ + have hmem : MemLp (fun y => y) 2 (tweedieMeasure μ (show (0:ℝ) ≤ φ by linarith) hp₁ hp₂) := by + rw [MeasureTheory.memLp_two_iff_integrable_sq (by fun_prop)] + exact tweedieMeasure_sq_integrable hμ hφ hp₁ hp₂ + rw [ProbabilityTheory.variance_eq_sub hmem] + simp only [Pi.pow_apply] + rw [tweedieMeasure_2nd_moment hμ hφ hp₁ hp₂, tweedieMeasure_expectation hμ hφ hp₁ hp₂] + ring diff --git a/Statlib/Tweedie/TweedieAux.lean b/Statlib/Tweedie/TweedieAux.lean new file mode 100644 index 0000000..32f5bbb --- /dev/null +++ b/Statlib/Tweedie/TweedieAux.lean @@ -0,0 +1,87 @@ +/- +Copyright (c) 2026 Bjørn Kjos-Hanssen. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Bjørn Kjos-Hanssen +-/ +import Mathlib + +/-! +# Mathlib-only auxiliary lemmas for the Tweedie / compound-Poisson equivalence + +This file collects purely analytic facts (no dependence on the Tweedie development) used in the +proof that the Tweedie measure equals the compound-Poisson construction: + +* `tw_scalar_identity`: the `rpow` bookkeeping identity relating the Poisson/gamma constants to the + Tweedie series coefficients. +* `ae_summable_of_summable_integral_norm`: if the integral norms of a family of functions are + summable, then the family is pointwise summable almost everywhere. +-/ + +open MeasureTheory Real +open scoped ENNReal NNReal + +namespace TweedieAux + +/- +The scalar `rpow` identity matching the Poisson-weighted gamma-rate term with the Tweedie series +coefficient. Here `α = (2-p)/(p-1)` (gamma shape), `γ = μ^(1-p)/(φ(p-1))` (gamma rate), +`λ = μ^(2-p)/(φ(2-p))` (Poisson rate), and `α' = (2-p)/(1-p) = -α`. +-/ +lemma tw_scalar_identity (μ φ p : ℝ) (hμ : 0 < μ) (hφ : 0 < φ) (hp₁ : 1 < p) (hp₂ : p < 2) + (j : ℕ) : + (μ ^ (2 - p) / (φ * (2 - p))) ^ j + * (μ ^ (1 - p) / (φ * (p - 1))) ^ ((j : ℝ) * ((2 - p) / (p - 1))) + = (p - 1) ^ (((2 - p) / (1 - p)) * (j : ℝ)) + / (φ ^ ((j : ℝ) * (1 - (2 - p) / (1 - p))) * (2 - p) ^ j) := by + -- Rewrite the left-hand side with positive bases + have lhs_rewrite : + (μ ^ (2 - p) / (φ * (2 - p))) ^ j * (μ ^ (1 - p) / (φ * (p - 1))) ^ (j * ((2 - p) / (p - 1))) + = (μ ^ ((2 - p) * j) * φ ^ (-(j : ℝ)) * (2 - p) ^ (-(j : ℝ))) * (μ ^ ((1 - p) * (j * ((2 - p) + / (p - 1)))) * φ ^ (-(j * ((2 - p) / (p - 1))) : ℝ) + * (p - 1) ^ (-(j * ((2 - p) / (p - 1))) : ℝ)) := by + congr 1; + · rw [Real.rpow_mul (by positivity), Real.rpow_neg (by linarith), Real.rpow_neg (by linarith)] + ring_nf + norm_cast; norm_num + ring_nf + rw [show - ( p * φ) + φ * 2 = φ * ( 2 - p) by ring, mul_inv]; ring; + · rw [Real.div_rpow (by positivity) (mul_nonneg hφ.le (by linarith))]; + rw [← Real.rpow_mul (by positivity), Real.mul_rpow (by positivity) (by linarith), + Real.rpow_neg (by linarith), Real.rpow_neg (by linarith)]; ring; + convert lhs_rewrite using 1; norm_num [Real.rpow_def_of_pos, *] + ring_nf + rw [← Real.exp_neg, ← Real.exp_add] + ring_nf + rw [← Real.exp_log (by linarith : 0 < 2 - p)]; norm_num [← Real.exp_add, ← Real.exp_nat_mul] + ring_nf + rw [← Real.exp_neg, ← Real.exp_add] + ring_nf + grind + +/- +If the family `F i` is integrable on `s` and the integral norms `∫_s ‖F i‖` are summable, then +for almost every `y ∈ s` the family `i ↦ F i y` is summable. +-/ +lemma ae_summable_of_summable_integral_norm {ι : Type*} [Countable ι] {F : ι → ℝ → ℝ} {s : Set ℝ} + (hF : ∀ i, IntegrableOn (F i) s) + (hsum : Summable (fun i => ∫ y in s, ‖F i y‖)) : + ∀ᵐ y ∂(volume.restrict s), Summable (fun i => F i y) := by + -- The total integral of `∑ ‖F i‖` is finite, hence the sum is finite a.e. + have h_fin : ∫⁻ y in s, ∑' i, ENNReal.ofReal (‖F i y‖) ∂volume < ⊤ := by + rw [MeasureTheory.lintegral_tsum] + · refine lt_of_le_of_lt (ENNReal.tsum_le_tsum fun i => ?_) + (Summable.tsum_ofReal_lt_top hsum) + rw [MeasureTheory.ofReal_integral_eq_lintegral_ofReal + (MeasureTheory.Integrable.norm (hF i)) + (Filter.Eventually.of_forall fun x => norm_nonneg _)] + · exact fun i => ENNReal.continuous_ofReal.measurable.comp_aemeasurable + ((hF i).aemeasurable.norm) + have h_ae : ∀ᵐ y ∂volume.restrict s, ∑' i, ENNReal.ofReal (‖F i y‖) < ⊤ := by + refine MeasureTheory.ae_lt_top' ?_ (ne_of_lt h_fin) + exact AEMeasurable.tsum (fun i => (hF i).aemeasurable.norm.ennreal_ofReal) + filter_upwards [h_ae] with y hy + refine Summable.of_norm ?_ + convert ENNReal.summable_toReal (ne_of_lt hy) with i + rw [ENNReal.toReal_ofReal (norm_nonneg _)] + +end TweedieAux