diff --git a/LMLExplore/Algebra/GeomSum.lean b/LMLExplore/Algebra/GeomSum.lean new file mode 100644 index 00000000..94e718d8 --- /dev/null +++ b/LMLExplore/Algebra/GeomSum.lean @@ -0,0 +1,34 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import Mathlib.Algebra.BigOperators.NatAntidiagonal +public import Mathlib.Algebra.Ring.GeomSum + +/-! +# Geometric sums over an antidiagonal + +The sum of `a ^ i * b ^ j` over the pairs `(i, j)` with `i + j = n` is the two-variable geometric +sum of Mathlib's `geom_sum₂_mul`: multiplied by `a - b`, it is `a ^ (n + 1) - b ^ (n + 1)`. + +## Main statements + +* `Finset.Nat.sum_antidiagonal_pow_mul_pow_mul_sub`: + `(∑ x ∈ antidiagonal n, a ^ x.1 * b ^ x.2) * (a - b) = a ^ (n + 1) - b ^ (n + 1)`. +-/ + +public section + +namespace Finset.Nat + +/-- **The geometric sum over an antidiagonal**, a form of `geom_sum₂_mul`: +`(∑ i + j = n, a ^ i * b ^ j) * (a - b) = a ^ (n + 1) - b ^ (n + 1)`. -/ +theorem sum_antidiagonal_pow_mul_pow_mul_sub {R : Type*} [CommRing R] (a b : R) (n : ℕ) : + (∑ x ∈ antidiagonal n, a ^ x.1 * b ^ x.2) * (a - b) = a ^ (n + 1) - b ^ (n + 1) := by + rw [sum_antidiagonal_eq_sum_range_succ_mk, ← geom_sum₂_mul] + simp + +end Finset.Nat diff --git a/LMLExplore/EValues/Numeraire/Conditional.lean b/LMLExplore/EValues/Numeraire/Conditional.lean new file mode 100644 index 00000000..491c5b39 --- /dev/null +++ b/LMLExplore/EValues/Numeraire/Conditional.lean @@ -0,0 +1,113 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import LMLExplore.EValues.Numeraire + +/-! +# Numeraires of nulls with a common conditional law given a statistic + +Let `T : Ω → 𝓨` be a statistic and `κ : Kernel 𝓨 Ω` a kernel such that `κ y` is carried by the +fiber `T ⁻¹' {y}`. Suppose that `κ` is the conditional law given `T` of every member `P` of a +null `S`, in the sense that `P = κ ∘ₘ P.map T`. For an alternative `Q`, the measure +`κ ∘ₘ Q.map T` has the law of `T` under `Q` and the conditional law `κ` of the null given `T`. + +The likelihood ratio of `Q` to `κ ∘ₘ Q.map T` integrates to one against `κ y` for `Q.map T`-almost +every `y`: its conditional expectation given `T` is one. So it is an e-variable for `S`, with +expectation exactly the mass `P univ` (one for probability measures) under every member `P` of `S` +such that the law of `T` under `P` is absolutely continuous with respect to its law under `Q`. This +is the conditional e-variable of [HGLLA24]. When moreover `κ ∘ₘ Q.map T` belongs to the effective +null of `S` (for example when it is a mixture of members of `S`, by +`ProbabilityTheory.comp_mem_effectiveSet`), this likelihood ratio is the numeraire of `S` under +`Q` by the criterion of [LRR25, Theorem 4.1] (`ProbabilityTheory.isNumeraire_of_ae_eq_rnDeriv`), +and `κ ∘ₘ Q.map T` is the reverse information projection of `Q` on `S` when it is equivalent to +`Q` (`ProbabilityTheory.ripr_eq_of_isEVar_rnDeriv`). + +## Main statements + +* `ProbabilityTheory.ae_lintegral_rnDeriv_comp_map_eq_one`: the likelihood ratio of `Q` to + `κ ∘ₘ Q.map T` has conditional expectation one given `T`. +* `ProbabilityTheory.lintegral_rnDeriv_comp_map`, `ProbabilityTheory.isEVar_rnDeriv_comp_map`: it + is an exact e-variable for every null of which `κ` is the conditional law given `T`. + +## References + +* [HGLLA24] Y. Hao, P. Grünwald, T. Lardy, L. Long and R. Adams, *E-values for k-sample tests + with exponential families*, arXiv:2303.00471. +* [LRR25] M. Larsson, A. Ramdas and J. Ruf, *The numeraire e-variable and reverse information + projection*, Ann. Statist. 53 (2025), arXiv:2402.18810v4, Theorem 4.1. +-/ + +public section + +open MeasureTheory Set + +open scoped ENNReal + +namespace ProbabilityTheory + +variable {Ω 𝓨 : Type*} {mΩ : MeasurableSpace Ω} {m𝓨 : MeasurableSpace 𝓨} {T : Ω → 𝓨} + {κ : Kernel 𝓨 Ω} {S : Set (Measure Ω)} {Q P : Measure Ω} + +/-- **The conditional expectation of the likelihood ratio given the statistic is one.** If `κ y` +is carried by the fiber `T ⁻¹' {y}` for `Q.map T`-almost every `y` and the finite measure `Q` is +absolutely continuous with respect to `κ ∘ₘ Q.map T`, the likelihood ratio of `Q` to +`κ ∘ₘ Q.map T` has integral one against `κ y` for `Q.map T`-almost every `y`. -/ +theorem ae_lintegral_rnDeriv_comp_map_eq_one [IsFiniteMeasure Q] [IsFiniteKernel κ] + (hT : Measurable T) (hκ : ∀ᵐ y ∂Q.map T, ∀ᵐ x ∂κ y, T x = y) (hQ : Q ≪ κ ∘ₘ Q.map T) : + ∀ᵐ y ∂Q.map T, ∫⁻ x, Q.rnDeriv (κ ∘ₘ Q.map T) x ∂κ y = 1 := by + set ν := Q.map T + set X := Q.rnDeriv (κ ∘ₘ ν) + have hX : Measurable X := Measure.measurable_rnDeriv _ _ + have hφ : Measurable fun y ↦ ∫⁻ x, X x ∂κ y := hX.lintegral_kernel + refine (withDensity_eq_iff_of_sigmaFinite (g := 1) hφ.aemeasurable aemeasurable_one).mp ?_ + rw [withDensity_one] + ext A hA + rw [withDensity_apply _ hA, Measure.map_apply hT hA, ← lintegral_indicator hA] + calc ∫⁻ y, A.indicator (fun y ↦ ∫⁻ x, X x ∂κ y) y ∂ν + _ = ∫⁻ y, ∫⁻ x, (T ⁻¹' A).indicator X x ∂κ y ∂ν := by + refine lintegral_congr_ae ?_ + filter_upwards [hκ] with y hκ + by_cases hy : y ∈ A + · rw [indicator_of_mem hy] + refine lintegral_congr_ae ?_ + filter_upwards [hκ] with x hx + rw [indicator_of_mem (show x ∈ T ⁻¹' A by simpa [hx])] + · rw [indicator_of_notMem hy, lintegral_congr_ae (g := 0), Pi.zero_def, lintegral_zero] + filter_upwards [hκ] with x hx + rw [indicator_of_notMem (show x ∉ T ⁻¹' A by simpa [hx]), Pi.zero_apply] + _ = ∫⁻ x, (T ⁻¹' A).indicator X x ∂(κ ∘ₘ ν) := + (Measure.lintegral_bind κ.aemeasurable (hX.indicator (hT hA)).aemeasurable).symm + _ = Q (T ⁻¹' A) := by + rw [lintegral_indicator (hT hA), ← withDensity_apply _ (hT hA), + Measure.withDensity_rnDeriv_eq _ _ hQ] + +/-- **The likelihood ratio to `κ ∘ₘ Q.map T` has expectation one under the null.** Under a +measure `P` whose conditional law given `T` is `κ` and such that the law of `T` under `P` is +absolutely continuous with respect to its law under `Q`, the likelihood ratio of `Q` to +`κ ∘ₘ Q.map T` has expectation the mass of `P`. -/ +theorem lintegral_rnDeriv_comp_map [IsFiniteMeasure Q] [IsFiniteKernel κ] (hT : Measurable T) + (hκ : ∀ᵐ y ∂Q.map T, ∀ᵐ x ∂κ y, T x = y) (hQ : Q ≪ κ ∘ₘ Q.map T) (hP : κ ∘ₘ P.map T = P) + (hPQ : P.map T ≪ Q.map T) : + ∫⁻ x, Q.rnDeriv (κ ∘ₘ Q.map T) x ∂P = P univ := by + conv_lhs => rw [← hP] + rw [Measure.lintegral_bind κ.aemeasurable (Measure.measurable_rnDeriv _ _).aemeasurable, + lintegral_congr_ae (hPQ.ae_le (ae_lintegral_rnDeriv_comp_map_eq_one hT hκ hQ)), + lintegral_one, Measure.map_apply hT .univ, preimage_univ] + +/-- **The conditional e-variable** [HGLLA24]. If `κ` is the conditional law given `T` of every +member of a null `S` of finite measures, under each of which the law of `T` is absolutely +continuous with respect to its law under `Q`, the likelihood ratio of `Q` to `κ ∘ₘ Q.map T` is an +e-variable for `S`, exact by `lintegral_rnDeriv_comp_map`. -/ +theorem isEVar_rnDeriv_comp_map [IsFiniteMeasure Q] [IsFiniteKernel κ] (hT : Measurable T) + (hκ : ∀ᵐ y ∂Q.map T, ∀ᵐ x ∂κ y, T x = y) (hQ : Q ≪ κ ∘ₘ Q.map T) + (hS : ∀ P ∈ S, IsFiniteMeasure P) (hSκ : ∀ P ∈ S, κ ∘ₘ P.map T = P) + (hSQ : ∀ P ∈ S, P.map T ≪ Q.map T) : + IsEVar (Q.rnDeriv (κ ∘ₘ Q.map T)) S := + .of_lintegral_le_measure_univ hS (Measure.measurable_rnDeriv _ _) fun P hP ↦ + (lintegral_rnDeriv_comp_map hT hκ hQ (hSκ P hP) (hSQ P hP)).le + +end ProbabilityTheory diff --git a/LMLExplore/EValues/TwoSample/Exponential.lean b/LMLExplore/EValues/TwoSample/Exponential.lean new file mode 100644 index 00000000..556318c7 --- /dev/null +++ b/LMLExplore/EValues/TwoSample/Exponential.lean @@ -0,0 +1,450 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import LMLExplore.EValues.Numeraire.Conditional +public import LMLExplore.Probability.Distributions.Exponential +public import TauCeti.Probability.Distributions.Uniform.Basic +import LMLExplore.Probability.Distributions.Gamma +import LMLExplore.Probability.Kernel.Composition.MeasureComp +import TauCeti.Probability.Distributions.Beta.Uniform +import TauCeti.Probability.Distributions.Gamma.Beta +import TauCeti.Probability.Distributions.Gamma.Measurability +import TauCeti.Probability.Distributions.Gamma.PDF + +/-! +# The numeraire of two exponential observations against a common rate + +Two independent observations `X₁, X₂` have the laws `Exp(l₁)` and `Exp(l₂)`, and the null +hypothesis is that they have a common, unknown rate: it is the set of the products +`Exp(l) ⊗ Exp(l)`, `l > 0`. This is the two-sample test with a common nuisance parameter +of [HGLLA24] for the exponential family. + +Under every member of the null, the sum `s = X₁ + X₂` has the law `Gamma(2, l)` and, given `s`, +the first observation is uniform on `[0, s]` and the second one is `s - X₁`: the conditional law +given the sum is `ProbabilityTheory.Kernel.segmentUniform s` for every member. For the alternative +`Q = Exp(l₁) ⊗ Exp(l₂)`, let `Qᶜ = Kernel.segmentUniform ∘ₘ Q.map (x.1 + x.2)` be the law under +which the sum has its law under `Q` and, given the sum `s`, the first observation is uniform on +`[0, s]`. Then: + +* the likelihood ratio of `Q` to `Qᶜ`, the conditional e-variable of [HGLLA24], is the numeraire + of the null under `Q`, and `Qᶜ` is the reverse information projection of `Q` on the null; +* for `l₁ ≠ l₂`, `Qᶜ` is the mixture of the members `Exp(1 / m) ⊗ Exp(1 / m)` of the null whose + mean `m` is uniform between `1 / l₁` and `1 / l₂`; +* no finite combination of members of the null is equal to it, against the conjecture of + [HGLLA24, §5.2]. + +These identifications, the mixture and the failure of finite combinations are [EVR26, Note 28]. + +## Main definitions + +* `ProbabilityTheory.Kernel.segmentUniform`: the law of `(u * s, (1 - u) * s)` for `u` uniform + on `[0, 1]`. + +## Main statements + +* `ProbabilityTheory.isNumeraire_rnDeriv_segmentUniform_comp`: the likelihood ratio of `Q` to + `Qᶜ` is the numeraire. +* `ProbabilityTheory.ripr_prod_expMeasure`: `Qᶜ` is the reverse information projection. +* `ProbabilityTheory.segmentUniform_comp_map_add_prod_expMeasure`: `Qᶜ` is the mixture of the + members of the null whose mean is uniform between the two means. +* `ProbabilityTheory.sum_smul_prod_expMeasure_ne_ripr`: no finite combination of members of the + null is the reverse information projection. + +## Implementation notes + +That the law of the sum and the conditional law given the sum are as described under the null is +the beta–gamma algebra of Tau Ceti (`TauCeti.Probability.map_div_add_prod_gammaMeasure`). The law +of the sum under `Q` is Tau Ceti's hypoexponential law +(`ProbabilityTheory.map_add_prod_expMeasure_eq_withDensity`), with density +`l₁ l₂ (exp (- l₁ x) - exp (- l₂ x)) / (l₂ - l₁)` on `(0, ∞)`, which is also the density of the +mixture of the laws `Gamma(2, 1 / m)` over the uniform mean `m` +(`ProbabilityTheory.map_add_prod_expMeasure_eq_bind`). The numeraire is then identified +through the effective null: the likelihood ratio is an e-variable by +`ProbabilityTheory.isEVar_rnDeriv_comp_map`, and `Qᶜ` is a mixture of members of the null, so +that `ProbabilityTheory.isNumeraire_of_ae_eq_rnDeriv` and +`ProbabilityTheory.ripr_eq_of_isEVar_rnDeriv` apply. For the last statement, a finite combination +of the laws `Gamma(2, l)` has a density `∑ c l * l ^ 2 * x * exp (- l x)`, which decays either too +slowly or too fast at infinity compared to `exp (- min l₁ l₂ * x)`. + +## References + +* [HGLLA24] Y. Hao, P. Grünwald, T. Lardy, L. Long and R. Adams, *E-values for k-sample tests + with exponential families*, arXiv:2303.00471, §5.2. +* [LRR25] M. Larsson, A. Ramdas and J. Ruf, *The numeraire e-variable and reverse information + projection*, Ann. Statist. 53 (2025), arXiv:2402.18810v4, Theorem 4.1. +* [EVR26] *Answers to open questions in e-values and anytime-valid inference*, paper project + `RemyDegenne/evalues_research` (2026), Note 28. +-/ + +public section + +open MeasureTheory Set Filter Real TauCeti.Probability + +open scoped ENNReal Topology + +namespace ProbabilityTheory + +/-! ### The uniform law on a segment -/ + +/-- **The uniform law on a segment**: for `s : ℝ`, the law of `(u * s, (1 - u) * s)` for `u` +uniform on `[0, 1]`. For `s ≥ 0`, its first coordinate is uniform on `[0, s]` and the two +coordinates add up to `s`. -/ +noncomputable def Kernel.segmentUniform : Kernel ℝ (ℝ × ℝ) := + (Kernel.id ×ₖ Kernel.const ℝ (uniformMeasure 0 1)).map fun p ↦ (p.2 * p.1, (1 - p.2) * p.1) + +/-- The uniform law on a segment is a Markov kernel. -/ +instance Kernel.isMarkovKernel_segmentUniform : IsMarkovKernel Kernel.segmentUniform := by + unfold Kernel.segmentUniform + exact Kernel.IsMarkovKernel.map _ (by fun_prop) + +/-- The uniform law on the segment of the pairs with sum `s`. -/ +@[simp] +theorem Kernel.segmentUniform_apply (s : ℝ) : + Kernel.segmentUniform s = (uniformMeasure 0 1).map fun u ↦ (u * s, (1 - u) * s) := by + rw [Kernel.segmentUniform, Kernel.map_apply _ (by fun_prop), Kernel.prod_apply, Kernel.id_apply, + Kernel.const_apply, Measure.dirac_prod, Measure.map_map (by fun_prop) (by fun_prop)] + rfl + +/-- The mixture of the uniform laws on the segments against a law `ν` of the sum is the image of +`uniformMeasure 0 1 ⊗ ν` under `(u, s) ↦ (u * s, (1 - u) * s)`. -/ +theorem Kernel.segmentUniform_comp (ν : Measure ℝ) [SFinite ν] : + Kernel.segmentUniform ∘ₘ ν = + ((uniformMeasure 0 1).prod ν).map fun p ↦ (p.1 * p.2, (1 - p.1) * p.2) := by + rw [Kernel.segmentUniform, ← Measure.map_comp _ _ (by fun_prop), ← Measure.compProd_eq_comp_prod, + Measure.compProd_const, ← Measure.prod_swap, Measure.map_map (by fun_prop) measurable_swap] + rfl + +/-- Under the uniform law on the segment of the pairs with sum `s`, the coordinates add up to +`s`. -/ +theorem Kernel.ae_add_segmentUniform (s : ℝ) : ∀ᵐ x ∂Kernel.segmentUniform s, x.1 + x.2 = s := by + rw [Kernel.segmentUniform_apply, + ae_map_iff (by fun_prop) (measurableSet_eq_fun (by fun_prop) measurable_const)] + exact ae_of_all _ fun u ↦ by ring + +/-! ### The null of two exponential observations with a common rate -/ + +section Null + +variable {l l₁ l₂ : ℝ} + +/-- **The conditional law given the sum under the null**: two independent `Exp(l)` variables are +`(u * s, (1 - u) * s)` for a sum `s` of law `Gamma(2, l)` and an independent uniform `u`. -/ +private lemma segmentUniform_comp_gammaMeasure (hl : 0 < l) : + Kernel.segmentUniform ∘ₘ gammaMeasure 2 l = (expMeasure l).prod (expMeasure l) := by + have h := map_div_add_prod_gammaMeasure one_pos one_pos hl + rw [betaMeasure_one_one, one_add_one_eq_two] at h + rw [Kernel.segmentUniform_comp, ← h, Measure.map_map (by fun_prop) (by fun_prop)] + conv_rhs => rw [← Measure.map_id (μ := (expMeasure l).prod (expMeasure l))] + refine Measure.map_congr ?_ + filter_upwards [ae_mem_prod_Ioi_gammaMeasure (a := 1) (r := l) (b := 1) (s := l)] with z hz + have hz0 : z.1 + z.2 ≠ 0 := (add_pos hz.1 hz.2).ne' + refine Prod.ext ?_ ?_ + · simp [field] + · simp [field] + +/-- Every member of the null is its own `Qᶜ`. -/ +private lemma segmentUniform_comp_map_add_prod_expMeasure_self (hl : 0 < l) : + Kernel.segmentUniform ∘ₘ ((expMeasure l).prod (expMeasure l)).map (fun x ↦ x.1 + x.2) = + (expMeasure l).prod (expMeasure l) := by + rw [map_add_prod_expMeasure_self hl, segmentUniform_comp_gammaMeasure hl] + +/-- `Exp(l₁) ⊗ Exp(l₂)` and its `Qᶜ` are equivalent: both are equivalent to `Exp(1) ⊗ Exp(1)`, +which is its own `Qᶜ`. -/ +private lemma absolutelyContinuous_segmentUniform_comp (h₁ : 0 < l₁) (h₂ : 0 < l₂) : + (expMeasure l₁).prod (expMeasure l₂) ≪ + Kernel.segmentUniform ∘ₘ ((expMeasure l₁).prod (expMeasure l₂)).map (fun x ↦ x.1 + x.2) ∧ + Kernel.segmentUniform ∘ₘ ((expMeasure l₁).prod (expMeasure l₂)).map (fun x ↦ x.1 + x.2) ≪ + (expMeasure l₁).prod (expMeasure l₂) := by + have h1 := segmentUniform_comp_map_add_prod_expMeasure_self (l := 1) one_pos + refine ⟨(prod_expMeasure_absolutelyContinuous l₁ l₂ one_pos one_pos).trans ?_, ?_⟩ + · rw [← h1] + exact ((prod_expMeasure_absolutelyContinuous 1 1 h₁ h₂).map (by fun_prop)).comp_right _ + · refine (((prod_expMeasure_absolutelyContinuous l₁ l₂ one_pos one_pos).map + (by fun_prop)).comp_right _).trans ?_ + rw [h1] + exact prod_expMeasure_absolutelyContinuous 1 1 h₁ h₂ + +/-- The members of the null are probability measures. -/ +private lemma isProbabilityMeasure_of_mem_null {P : Measure (ℝ × ℝ)} + (hP : P ∈ Set.range fun l : Ioi (0 : ℝ) ↦ (expMeasure l).prod (expMeasure l)) : + IsProbabilityMeasure P := by + obtain ⟨⟨l, hl⟩, rfl⟩ := hP + have := isProbabilityMeasure_expMeasure hl + infer_instance + +/-- Every member of the null has the conditional law `Kernel.segmentUniform` given the sum. -/ +private lemma segmentUniform_comp_map_add_of_mem {P : Measure (ℝ × ℝ)} + (hP : P ∈ Set.range fun l : Ioi (0 : ℝ) ↦ (expMeasure l).prod (expMeasure l)) : + Kernel.segmentUniform ∘ₘ P.map (fun x ↦ x.1 + x.2) = P := by + obtain ⟨⟨l, hl⟩, rfl⟩ := hP + exact segmentUniform_comp_map_add_prod_expMeasure_self hl + +/-- Under every member of the null, the law of the sum is absolutely continuous with respect to +its law under `Exp(l₁) ⊗ Exp(l₂)`. -/ +private lemma map_add_absolutelyContinuous_of_mem (h₁ : 0 < l₁) (h₂ : 0 < l₂) + {P : Measure (ℝ × ℝ)} + (hP : P ∈ Set.range fun l : Ioi (0 : ℝ) ↦ (expMeasure l).prod (expMeasure l)) : + P.map (fun x ↦ x.1 + x.2) ≪ + ((expMeasure l₁).prod (expMeasure l₂)).map (fun x ↦ x.1 + x.2) := by + obtain ⟨⟨l, hl⟩, rfl⟩ := hP + exact (prod_expMeasure_absolutelyContinuous l l h₁ h₂).map (by fun_prop) + +/-! ### The law of the sum and the mixture -/ + +/-- The density of the sum is nonnegative on `[0, ∞)`. -/ +private lemma sumDensity_nonneg (h₁ : 0 < l₁) (h₂ : 0 < l₂) {x : ℝ} (hx : 0 ≤ x) : + 0 ≤ l₁ * l₂ * (exp (-(l₁ * x)) - exp (-(l₂ * x))) / (l₂ - l₁) := by + rcases le_total l₁ l₂ with h | h + · exact div_nonneg (mul_nonneg (by positivity) (sub_nonneg.mpr (exp_le_exp.mpr + (by nlinarith)))) (sub_nonneg.mpr h) + · exact div_nonneg_of_nonpos (mul_nonpos_of_nonneg_of_nonpos (by positivity) + (sub_nonpos.mpr (exp_le_exp.mpr (by nlinarith)))) (sub_nonpos.mpr h) + +/-- The laws `Gamma(2, 1 / m)`, as a kernel. -/ +private noncomputable def gammaTwoInvKernel : Kernel ℝ ℝ where + toFun m := gammaMeasure 2 m⁻¹ + measurable' := measurable_gammaMeasure.comp (f := fun m : ℝ ↦ ((2 : ℝ), m⁻¹)) (by fun_prop) + +/-- The kernel `gammaTwoInvKernel` maps `m` to `Gamma(2, 1 / m)`. -/ +private lemma gammaTwoInvKernel_apply (m : ℝ) : gammaTwoInvKernel m = gammaMeasure 2 m⁻¹ := rfl + +/-- The mixture of the laws `Gamma(2, 1 / m)` over a law of `m`, as a composition with +`gammaTwoInvKernel`. -/ +private lemma bind_gammaMeasure_eq_comp (μ : Measure ℝ) : + (μ.bind fun m ↦ gammaMeasure 2 m⁻¹) = gammaTwoInvKernel ∘ₘ μ := by + simp_rw [← gammaTwoInvKernel_apply] + +/-- `Qᶜ` is in the effective null. -/ +private lemma segmentUniform_comp_mem_effectiveSet (h₁ : 0 < l₁) (h₂ : 0 < l₂) : + Kernel.segmentUniform ∘ₘ ((expMeasure l₁).prod (expMeasure l₂)).map (fun x ↦ x.1 + x.2) ∈ + effectiveSet (Set.range fun l : Ioi (0 : ℝ) ↦ (expMeasure l).prod (expMeasure l)) := by + have hS : ∀ P ∈ Set.range fun l : Ioi (0 : ℝ) ↦ (expMeasure l).prod (expMeasure l), + P univ ≤ 1 := fun P hP ↦ (isProbabilityMeasure_of_mem_null hP).measure_univ.le + rcases eq_or_ne l₁ l₂ with rfl | h + · rw [segmentUniform_comp_map_add_prod_expMeasure_self h₁] + exact subset_effectiveSet hS ⟨⟨l₁, h₁⟩, rfl⟩ + have hα : 0 < min l₁⁻¹ l₂⁻¹ := lt_min (inv_pos.mpr h₁) (inv_pos.mpr h₂) + have := isProbabilityMeasure_uniformMeasure + (min_lt_max.mpr (inv_injective.ne h) : min l₁⁻¹ l₂⁻¹ < max l₁⁻¹ l₂⁻¹) + rw [map_add_prod_expMeasure_eq_bind h₁ h₂ h, bind_gammaMeasure_eq_comp, Measure.comp_assoc] + refine comp_mem_effectiveSet prob_le_one ?_ + filter_upwards [ae_mem_Ioc_uniformMeasure] with m hm + have hm0 : 0 < m⁻¹ := inv_pos.mpr (hα.trans hm.1) + rw [Kernel.comp_apply, gammaTwoInvKernel_apply, segmentUniform_comp_gammaMeasure hm0] + exact subset_effectiveSet hS ⟨⟨m⁻¹, hm0⟩, rfl⟩ + +/-- The likelihood ratio of `Exp(l₁) ⊗ Exp(l₂)` to its `Qᶜ` is an e-variable for the null. -/ +private lemma isEVar_rnDeriv_segmentUniform_comp (h₁ : 0 < l₁) (h₂ : 0 < l₂) : + IsEVar (((expMeasure l₁).prod (expMeasure l₂)).rnDeriv + (Kernel.segmentUniform ∘ₘ ((expMeasure l₁).prod (expMeasure l₂)).map fun x ↦ x.1 + x.2)) + (Set.range fun l : Ioi (0 : ℝ) ↦ (expMeasure l).prod (expMeasure l)) := by + have := isProbabilityMeasure_expMeasure h₁ + have := isProbabilityMeasure_expMeasure h₂ + exact isEVar_rnDeriv_comp_map (by fun_prop) (ae_of_all _ Kernel.ae_add_segmentUniform) + (absolutelyContinuous_segmentUniform_comp h₁ h₂).1 + (fun P hP ↦ have := isProbabilityMeasure_of_mem_null hP; inferInstance) + (fun _ ↦ segmentUniform_comp_map_add_of_mem) + (fun _ ↦ map_add_absolutelyContinuous_of_mem h₁ h₂) + +end Null + +/-! ### The numeraire and the reverse information projection -/ + +section Main + +variable {l₁ l₂ : ℝ} + +/-- **The numeraire of two exponential observations against a common rate.** For positive rates +`l₁, l₂`, the likelihood ratio of `Exp(l₁) ⊗ Exp(l₂)` to the law under which the sum of the two +observations has its law under `Exp(l₁) ⊗ Exp(l₂)` and, given the sum `s`, the first observation +is uniform on `[0, s]`, is the numeraire of the null of the products `Exp(l) ⊗ Exp(l)`, `l > 0`: +it is an e-variable, and no e-variable has a larger expected logarithm under +`Exp(l₁) ⊗ Exp(l₂)` (`ProbabilityTheory.IsNumeraire.eintegral_log_le`) [EVR26, Note 28]. -/ +theorem isNumeraire_rnDeriv_segmentUniform_comp (h₁ : 0 < l₁) (h₂ : 0 < l₂) : + IsNumeraire (((expMeasure l₁).prod (expMeasure l₂)).rnDeriv + (Kernel.segmentUniform ∘ₘ ((expMeasure l₁).prod (expMeasure l₂)).map fun x ↦ x.1 + x.2)) + (Set.range fun l : Ioi (0 : ℝ) ↦ (expMeasure l).prod (expMeasure l)) + ((expMeasure l₁).prod (expMeasure l₂)) := by + have := isProbabilityMeasure_expMeasure h₁ + have := isProbabilityMeasure_expMeasure h₂ + exact isNumeraire_of_ae_eq_rnDeriv (isEVar_rnDeriv_segmentUniform_comp h₁ h₂) + (segmentUniform_comp_mem_effectiveSet h₁ h₂) (absolutelyContinuous_segmentUniform_comp h₁ h₂).1 + .rfl + +/-- **The reverse information projection of two exponential observations** on the null of the +products `Exp(l) ⊗ Exp(l)`, `l > 0`, is the law under which the sum of the two observations has +its law under `Exp(l₁) ⊗ Exp(l₂)` and, given the sum `s`, the first observation is uniform on +`[0, s]` [EVR26, Note 28]. -/ +theorem ripr_prod_expMeasure (h₁ : 0 < l₁) (h₂ : 0 < l₂) : + ripr ((expMeasure l₁).prod (expMeasure l₂)) + (Set.range fun l : Ioi (0 : ℝ) ↦ (expMeasure l).prod (expMeasure l)) = + Kernel.segmentUniform ∘ₘ ((expMeasure l₁).prod (expMeasure l₂)).map fun x ↦ x.1 + x.2 := by + have := isProbabilityMeasure_expMeasure h₁ + have := isProbabilityMeasure_expMeasure h₂ + exact ripr_eq_of_isEVar_rnDeriv + (fun P hP ↦ have := isProbabilityMeasure_of_mem_null hP; inferInstance) + (segmentUniform_comp_mem_effectiveSet h₁ h₂) (absolutelyContinuous_segmentUniform_comp h₁ h₂).1 + (absolutelyContinuous_segmentUniform_comp h₁ h₂).2 (isEVar_rnDeriv_segmentUniform_comp h₁ h₂) + +/-- **The reverse information projection is a continuous mixture of members of the null.** For +distinct positive rates `l₁, l₂`, the law under which the sum of the two observations has its +law under `Exp(l₁) ⊗ Exp(l₂)` and, given the sum `s`, the first observation is uniform on +`[0, s]`, is the mixture of the members `Exp(1 / m) ⊗ Exp(1 / m)` of the null whose mean `m` is +uniform between `1 / l₁` and `1 / l₂` [EVR26, Note 28]. -/ +theorem segmentUniform_comp_map_add_prod_expMeasure (h₁ : 0 < l₁) (h₂ : 0 < l₂) (h : l₁ ≠ l₂) : + Kernel.segmentUniform ∘ₘ ((expMeasure l₁).prod (expMeasure l₂)).map (fun x ↦ x.1 + x.2) = + (uniformMeasure (min l₁⁻¹ l₂⁻¹) (max l₁⁻¹ l₂⁻¹)).bind + fun m ↦ (expMeasure m⁻¹).prod (expMeasure m⁻¹) := by + have hα : 0 < min l₁⁻¹ l₂⁻¹ := lt_min (inv_pos.mpr h₁) (inv_pos.mpr h₂) + rw [map_add_prod_expMeasure_eq_bind h₁ h₂ h, Measure.bind_bind (f := fun m ↦ gammaMeasure 2 m⁻¹) + gammaTwoInvKernel.measurable.aemeasurable Kernel.segmentUniform.aemeasurable] + refine Measure.bind_congr_right ?_ + filter_upwards [ae_mem_Ioc_uniformMeasure] with m hm + exact segmentUniform_comp_gammaMeasure (inv_pos.mpr (hα.trans hm.1)) + +/-- A nonnegative finite combination of the functions `x ↦ l ^ 2 * x * exp (- l x)`, `l > 0`, is +not `C (exp (- a x) - exp (- b x))` on `(0, ∞)` for `a < b` and `C > 0`: at infinity, the terms +with `l ≤ a` decay too slowly and the other ones too fast. -/ +private lemma not_forall_sum_eq_mul_sub {a b C : ℝ} (hab : a < b) (hC : 0 < C) (s : Finset ℝ) + (c : ℝ → ℝ) (hc : ∀ l ∈ s, 0 ≤ c l) (hs : ∀ l ∈ s, 0 < l) : + ¬ ∀ x > 0, ∑ l ∈ s, c l * (l ^ 2 * x * exp (-(l * x))) = + C * (exp (-(a * x)) - exp (-(b * x))) := by + intro h + have hterm : ∀ x > 0, ∀ l ∈ s, 0 ≤ c l * (l ^ 2 * x * exp (-(l * x))) := fun x hx l hl ↦ by + have := hc l hl + positivity + -- the coefficients of the rates `l ≤ a` vanish + have h0 : ∀ l ∈ s, l ≤ a → c l = 0 := by + intro l hl hla + by_contra hcl + have hcl' : 0 < c l * l ^ 2 := mul_pos ((hc l hl).lt_of_ne' hcl) (pow_pos (hs l hl) 2) + set x := 2 * C / (c l * l ^ 2) with hx_def + have hx : 0 < x := by positivity + have h1 := Finset.single_le_sum (hterm x hx) hl + rw [h x hx] at h1 + have h2 : C * exp (-(a * x)) ≤ C * exp (-(l * x)) := + mul_le_mul_of_nonneg_left (exp_le_exp.mpr (by nlinarith)) hC.le + have h3 : 0 < C * exp (-(b * x)) := mul_pos hC (exp_pos _) + have h4 : 0 < C * exp (-(l * x)) := mul_pos hC (exp_pos _) + have h5 : c l * (l ^ 2 * x * exp (-(l * x))) = 2 * (C * exp (-(l * x))) := by + have := (hs l hl).ne' + rw [hx_def] + field_simp + linarith + -- multiply by `exp (a x)` and let `x` tend to infinity + have hlim : Tendsto (fun x ↦ ∑ l ∈ s, c l * l ^ 2 * (x * exp (-((l - a) * x)))) atTop + (𝓝 0) := by + rw [show (0 : ℝ) = ∑ _l ∈ s, 0 by simp] + refine tendsto_finsetSum _ fun l hl ↦ ?_ + by_cases hla : l ≤ a + · simp [h0 l hl hla] + have hε : 0 < l - a := sub_pos.mpr (not_le.mp hla) + have := ((tendsto_pow_mul_exp_neg_atTop_nhds_zero 1).comp + (Tendsto.const_mul_atTop hε tendsto_id)).const_mul (c l * l ^ 2 / (l - a)) + rw [mul_zero] at this + refine this.congr fun x ↦ ?_ + simp only [Function.comp_apply, id, pow_one] + field_simp + have hlim' : Tendsto (fun x ↦ C * (1 - exp (-((b - a) * x)))) atTop (𝓝 C) := by + have : Tendsto (fun x ↦ exp (-((b - a) * x))) atTop (𝓝 0) := + tendsto_exp_neg_atTop_nhds_zero.comp (Tendsto.const_mul_atTop (sub_pos.mpr hab) tendsto_id) + simpa using (this.const_sub 1).const_mul C + refine hC.ne (tendsto_nhds_unique (hlim.congr' ?_) hlim') + filter_upwards [eventually_gt_atTop 0] with x hx + have key : ∀ l, exp (-((l - a) * x)) = exp (-(l * x)) * exp (a * x) := fun l ↦ by + rw [← exp_add] + ring_nf + have ha : exp (-(a * x)) * exp (a * x) = 1 := by + rw [← exp_add] + simp + calc ∑ l ∈ s, c l * l ^ 2 * (x * exp (-((l - a) * x))) + _ = (∑ l ∈ s, c l * (l ^ 2 * x * exp (-(l * x)))) * exp (a * x) := by + rw [Finset.sum_mul] + refine Finset.sum_congr rfl fun l _ ↦ ?_ + rw [key] + ring + _ = C * (1 - exp (-((b - a) * x))) := by + rw [h x hx, key b] + linear_combination C * ha + +/-- **No finite combination of members of the null is the reverse information projection.** For +distinct positive rates `l₁, l₂`, no finite combination of the laws `Exp(l) ⊗ Exp(l)` is the +reverse information projection of `Exp(l₁) ⊗ Exp(l₂)` on the null of the products +`Exp(l) ⊗ Exp(l)`, `l > 0`. In particular it is not a mixture of two members of the null, against +the conjecture of [HGLLA24, §5.2]; this refutation is [EVR26, Note 28]. -/ +theorem sum_smul_prod_expMeasure_ne_ripr (h₁ : 0 < l₁) (h₂ : 0 < l₂) (h : l₁ ≠ l₂) + (s : Finset ℝ) (w : ℝ → ℝ≥0∞) : + ∑ l ∈ s, w l • (expMeasure l).prod (expMeasure l) ≠ + ripr ((expMeasure l₁).prod (expMeasure l₂)) + (Set.range fun l : Ioi (0 : ℝ) ↦ (expMeasure l).prod (expMeasure l)) := by + have := isProbabilityMeasure_expMeasure h₁ + have := isProbabilityMeasure_expMeasure h₂ + rw [ripr_prod_expMeasure h₁ h₂] + intro heq + set s' := s.filter (0 < ·) with hs' + have hs'pos : ∀ l ∈ s', 0 < l := fun l hl ↦ (Finset.mem_filter.mp hl).2 + -- the laws of the sum agree + have hsum : ∀ A, MeasurableSet A → ∑ l ∈ s', w l * gammaMeasure 2 l A = + ((expMeasure l₁).prod (expMeasure l₂)).map (fun x ↦ x.1 + x.2) A := by + intro A hA + have h1 := congrArg (fun μ : Measure (ℝ × ℝ) ↦ μ ((fun x ↦ x.1 + x.2) ⁻¹' A)) heq + simp only [FunLike.coe_sum, Finset.sum_apply, smul_apply, smul_eq_mul] at h1 + rw [← Measure.map_apply (by fun_prop) hA, + map_comp_eq_of_ae_eq (by fun_prop) (ae_of_all _ Kernel.ae_add_segmentUniform)] at h1 + rw [← h1, hs', Finset.sum_filter] + refine Finset.sum_congr rfl fun l _ ↦ ?_ + split_ifs with hl + · rw [← map_add_prod_expMeasure_self hl, Measure.map_apply (by fun_prop) hA] + · simp [expMeasure_of_nonpos (not_lt.mp hl)] + -- so the weights are finite + have hw : ∀ l ∈ s', w l ≠ ∞ := by + have hmass := hsum univ MeasurableSet.univ + rw [measure_univ, Finset.sum_congr rfl fun l hl ↦ by + have := isProbabilityMeasure_gammaMeasure two_pos (hs'pos l hl) + rw [measure_univ, mul_one]] at hmass + exact fun l hl ↦ ne_top_of_le_ne_top ENNReal.one_ne_top + (hmass ▸ Finset.single_le_sum (fun _ _ ↦ zero_le) hl) + -- and the densities of the laws of the sum agree on `(0, ∞)` + have hae : (fun x ↦ ∑ l ∈ s', w l * gammaPDF 2 l x) =ᵐ[volume.restrict (Ioi 0)] + fun x ↦ ENNReal.ofReal (l₁ * l₂ * (exp (-(l₁ * x)) - exp (-(l₂ * x))) / (l₂ - l₁)) := by + refine ae_eq_of_forall_setLIntegral_eq_of_sigmaFinite + (Finset.measurable_fun_sum _ fun l _ ↦ (measurable_gammaPDF 2 l).const_mul _) + (by fun_prop) fun A hA _ ↦ ?_ + rw [lintegral_finsetSum _ fun l _ ↦ (measurable_gammaPDF 2 l).const_mul _, + ← withDensity_apply _ hA, ← map_add_prod_expMeasure_eq_withDensity h₁ h₂ h, ← hsum A hA] + refine Finset.sum_congr rfl fun l _ ↦ ?_ + rw [lintegral_const_mul _ (measurable_gammaPDF 2 l), ← withDensity_apply _ hA, + ← gammaMeasure_eq_withDensity_restrict_Ioi] + have hreal : ∀ x ∈ Ioi (0 : ℝ), + ∑ l ∈ s', (w l).toReal * (l ^ 2 * x * exp (-(l * x))) = + l₁ * l₂ * (exp (-(l₁ * x)) - exp (-(l₂ * x))) / (l₂ - l₁) := by + refine Measure.eqOn_open_of_ae_eq (μ := volume) ?_ isOpen_Ioi (by fun_prop) (by fun_prop) + filter_upwards [hae, ae_restrict_mem measurableSet_Ioi] with x hx (hx0 : 0 < x) + rw [← ENNReal.ofReal_eq_ofReal_iff (Finset.sum_nonneg fun l hl ↦ + have := hs'pos l hl; by positivity) (sumDensity_nonneg h₁ h₂ hx0.le), ← hx, + ENNReal.ofReal_sum_of_nonneg fun l hl ↦ have := hs'pos l hl; by positivity] + refine Finset.sum_congr rfl fun l hl ↦ ?_ + rw [ENNReal.ofReal_mul ENNReal.toReal_nonneg, ENNReal.ofReal_toReal (hw l hl), + gammaPDF_two hx0.le] + -- the density of the law of the sum is symmetric in the two rates: order them + obtain ⟨a, b, hab, ha, hdens⟩ : ∃ a b, a < b ∧ 0 < a ∧ ∀ x, + l₁ * l₂ * (exp (-(l₁ * x)) - exp (-(l₂ * x))) / (l₂ - l₁) = + a * b / (b - a) * (exp (-(a * x)) - exp (-(b * x))) := by + have h12 : l₁ - l₂ ≠ 0 := sub_ne_zero.mpr h + have h21 : l₂ - l₁ ≠ 0 := sub_ne_zero.mpr h.symm + rcases h.lt_or_gt with h' | h' + · exact ⟨l₁, l₂, h', h₁, fun x ↦ by ring⟩ + · refine ⟨l₂, l₁, h', h₂, fun x ↦ ?_⟩ + field_simp + ring + exact not_forall_sum_eq_mul_sub hab (div_pos (mul_pos ha (ha.trans hab)) (sub_pos.mpr hab)) s' + (fun l ↦ (w l).toReal) (fun _ _ ↦ ENNReal.toReal_nonneg) hs'pos fun x hx ↦ by + rw [hreal x hx, hdens] + +end Main + +end ProbabilityTheory diff --git a/LMLExplore/EValues/TwoSample/Geometric.lean b/LMLExplore/EValues/TwoSample/Geometric.lean new file mode 100644 index 00000000..ca077712 --- /dev/null +++ b/LMLExplore/EValues/TwoSample/Geometric.lean @@ -0,0 +1,344 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import LMLExplore.EValues.Numeraire.Conditional +public import Mathlib.Probability.Distributions.Geometric +public import Mathlib.Probability.UniformOn +import LMLExplore.Algebra.GeomSum +import LMLExplore.MeasureTheory.Measure.Singletons +import LMLExplore.Probability.Distributions.Geometric +import LMLExplore.Probability.Distributions.Uniform +import TauCeti.Probability.Distributions.Geometric.Measurability +import TauCeti.Probability.Distributions.Uniform.Basic + +/-! +# The numeraire of two geometric observations against a common success probability + +Two independent observations `X₁, X₂` in `ℕ` have the geometric laws `Geom(p₁)` and `Geom(p₂)` +of Mathlib's `ProbabilityTheory.geometricMeasure`, and the null hypothesis is that they have a +common, unknown success probability: it is the set of the products `Geom(p) ⊗ Geom(p)`. Under +every member of the null the mass of a pair `(j, k)` depends only on `j + k`, so, given the sum +`s`, the first observation is uniform on `{0, …, s}`: the conditional law given the sum is +`ProbabilityTheory.Kernel.antidiagonalUniform s` for every member. + +For the alternative `Q = Geom(p₁) ⊗ Geom(p₂)`, the likelihood ratio of `Q` to the law +`Kernel.antidiagonalUniform ∘ₘ Q.map (x.1 + x.2)`, under which the sum has its law under `Q` and, +given the sum `s`, the first observation is uniform on `{0, …, s}`, is the numeraire of the null +under `Q` (`ProbabilityTheory.isNumeraire_rnDeriv_antidiagonalUniform_comp`), as computed in +[EVR26, Note 28]. As for two exponential observations +(`ProbabilityTheory.isNumeraire_rnDeriv_segmentUniform_comp`), for distinct positive success +probabilities `p₁ ≠ p₂` that law is the mixture of the members `Geom(1 / m) ⊗ Geom(1 / m)` of the +null over `m` uniform between `1 / p₁` and `1 / p₂`, that is, over the mean `m - 1` uniform between +the two means; when `p₁ = p₂` it is `Q` itself. In both cases it belongs to the effective null. + +## Main definitions + +* `ProbabilityTheory.Kernel.antidiagonalUniform`: the uniform law on the pairs `(j, k)` of natural + numbers with `j + k = s`. + +## Main statements + +* `ProbabilityTheory.Kernel.antidiagonalUniform_comp_singleton`, + `ProbabilityTheory.Kernel.antidiagonalUniform_comp_map_add_of_singleton`: the masses of the + points under `Kernel.antidiagonalUniform ∘ₘ ν`, and a law on pairs whose masses depend only on + the sum has the conditional law `Kernel.antidiagonalUniform` given the sum. +* `ProbabilityTheory.isNumeraire_rnDeriv_antidiagonalUniform_comp`: the numeraire of two geometric + observations against a common success probability [EVR26, Note 28]. + +## Implementation notes + +Mathlib's `geometricMeasure 0` is the point mass at `0`, as is `geometricMeasure 1`, so the +statement holds for every pair of success probabilities except those for which `Q` is the point +mass at `(0, 0)`. The mixture is computed on the masses of the points, through the antiderivative +`(1 - 1 / m) ^ (s + 1) / (s + 1)` of `(1 - 1 / m) ^ s / m ^ 2`. + +## References + +* [HGLLA24] Y. Hao, P. Grünwald, T. Lardy, L. Long and R. Adams, *E-values for k-sample tests + with exponential families*, arXiv:2303.00471. +* [LRR25] M. Larsson, A. Ramdas and J. Ruf, *The numeraire e-variable and reverse information + projection*, Ann. Statist. 53 (2025), arXiv:2402.18810v4, Theorem 4.1. +* [EVR26] *Answers to open questions in e-values and anytime-valid inference*, paper project + `RemyDegenne/evalues_research` (2026), Note 28. +-/ + +public section + +open MeasureTheory Set Finset TauCeti.Probability + +open scoped ENNReal + +namespace ProbabilityTheory + +/-! ### The uniform law on an antidiagonal -/ + +/-- **The uniform law on an antidiagonal**: for `s : ℕ`, the uniform law on the `s + 1` pairs +`(j, k)` of natural numbers with `j + k = s`. Its first coordinate is uniform on `{0, …, s}`. -/ +noncomputable def Kernel.antidiagonalUniform : Kernel ℕ (ℕ × ℕ) := + Kernel.ofFunOfCountable fun s ↦ uniformOn (antidiagonal s : Set (ℕ × ℕ)) + +/-- The uniform law on the pairs with sum `s`. -/ +@[simp] +theorem Kernel.antidiagonalUniform_apply (s : ℕ) : + Kernel.antidiagonalUniform s = uniformOn (antidiagonal s : Set (ℕ × ℕ)) := by + unfold Kernel.antidiagonalUniform + rfl + +/-- The uniform law on an antidiagonal is a Markov kernel. -/ +instance Kernel.isMarkovKernel_antidiagonalUniform : IsMarkovKernel Kernel.antidiagonalUniform := + ⟨fun s ↦ by + rw [Kernel.antidiagonalUniform_apply] + exact isProbabilityMeasure_uniformOn (finite_toSet _) ⟨(0, s), by simp⟩⟩ + +/-- The mass of a point under the uniform law on the pairs with sum `s`. -/ +theorem Kernel.antidiagonalUniform_singleton (s : ℕ) (x : ℕ × ℕ) : + Kernel.antidiagonalUniform s {x} = if x.1 + x.2 = s then ((s : ℝ≥0∞) + 1)⁻¹ else 0 := by + classical + rw [Kernel.antidiagonalUniform_apply, ← coe_singleton, uniformOn_apply_finset, + Nat.card_antidiagonal] + split_ifs with hx + · rw [Finset.inter_singleton_of_mem (mem_antidiagonal.mpr hx)] + simp + · rw [Finset.inter_singleton_of_notMem (by simpa using hx)] + simp + +/-- Under the uniform law on the pairs with sum `s`, the coordinates add up to `s`. -/ +theorem Kernel.ae_add_antidiagonalUniform (s : ℕ) : + ∀ᵐ x ∂Kernel.antidiagonalUniform s, x.1 + x.2 = s := by + rw [ae_iff, Kernel.antidiagonalUniform_apply, ← uniformOn_inter_self (finite_toSet _)] + convert measure_empty (μ := uniformOn (antidiagonal s : Set (ℕ × ℕ))) + ext x + simp + +/-! ### Laws on pairs whose masses depend on the sum -/ + +section Pairs + +variable {μ : Measure (ℕ × ℕ)} + +/-- The mass of a point under `Kernel.antidiagonalUniform ∘ₘ ν`. -/ +theorem Kernel.antidiagonalUniform_comp_singleton (ν : Measure ℕ) (x : ℕ × ℕ) : + (Kernel.antidiagonalUniform ∘ₘ ν) {x} = + ν {x.1 + x.2} * (((x.1 + x.2 : ℕ) : ℝ≥0∞) + 1)⁻¹ := by + rw [Measure.bind_apply (measurableSet_singleton x) (Kernel.aemeasurable _), lintegral_countable'] + simp_rw [Kernel.antidiagonalUniform_singleton] + rw [tsum_eq_single (x.1 + x.2) fun s hs ↦ by simp [Ne.symm hs]] + simp [mul_comm] + +/-- A law on pairs whose mass at a point depends only on the sum has the conditional law +`Kernel.antidiagonalUniform` given the sum. -/ +theorem Kernel.antidiagonalUniform_comp_map_add_of_singleton {g : ℕ → ℝ≥0∞} + (h : ∀ x, μ {x} = g (x.1 + x.2)) : + Kernel.antidiagonalUniform ∘ₘ μ.map (fun x ↦ x.1 + x.2) = μ := by + refine Measure.ext_iff_singleton.mpr fun x ↦ ?_ + rw [Kernel.antidiagonalUniform_comp_singleton, Measure.map_add_singleton_eq_sum_antidiagonal, + Finset.sum_congr rfl fun y hy ↦ by rw [h y, mem_antidiagonal.mp hy], Finset.sum_const, + Nat.card_antidiagonal, nsmul_eq_mul, h x, mul_comm _ (g _), mul_assoc] + push_cast + rw [ENNReal.mul_inv_cancel (by positivity) (by simp), mul_one] + +/-- Every law on pairs is absolutely continuous with respect to the law with its law of the sum +and the conditional law `Kernel.antidiagonalUniform` given the sum. -/ +theorem Kernel.absolutelyContinuous_antidiagonalUniform_comp (μ : Measure (ℕ × ℕ)) : + μ ≪ Kernel.antidiagonalUniform ∘ₘ μ.map (fun x ↦ x.1 + x.2) := by + refine Measure.AbsolutelyContinuous.mk fun t _ ht ↦ ?_ + rw [← biUnion_of_singleton t, measure_biUnion_null_iff t.to_countable] at ht ⊢ + intro x hx + have hx' := ht x hx + rw [Kernel.antidiagonalUniform_comp_singleton, mul_eq_zero, + Measure.map_add_singleton_eq_sum_antidiagonal, Finset.sum_eq_zero_iff] at hx' + exact hx'.resolve_right (by simp) x (mem_antidiagonal.mpr rfl) + +end Pairs + +/-! ### Geometric laws -/ + +section Geometric + +variable {p₁ p₂ : unitInterval} + +/-- The mass of a point under `Geom(p) ⊗ Geom(p)` depends only on the sum. -/ +private lemma exists_prod_geometricMeasure_self_singleton (p : unitInterval) : + ∃ g : ℕ → ℝ≥0∞, ∀ x : ℕ × ℕ, + (geometricMeasure p).prod (geometricMeasure p) {x} = g (x.1 + x.2) := by + obtain ⟨q, hq, hqp⟩ := exists_ne_zero_geometricMeasure_eq p + refine ⟨fun s ↦ ENNReal.ofReal ((1 - q) ^ s * q ^ 2), fun x ↦ ?_⟩ + rw [← hqp, prod_geometricMeasure_singleton hq hq] + congr 1 + ring + +/-- If one of the two geometric laws is not a point mass, the sum charges every point. -/ +private lemma map_add_prod_geometricMeasure_singleton_pos + (h : (p₁ ≠ 0 ∧ p₁ ≠ 1) ∨ (p₂ ≠ 0 ∧ p₂ ≠ 1)) (s : ℕ) : + 0 < ((geometricMeasure p₁).prod (geometricMeasure p₂)).map (fun x ↦ x.1 + x.2) {s} := by + rw [Measure.map_add_singleton_eq_sum_antidiagonal] + rcases h with ⟨h0, h1⟩ | ⟨h0, h1⟩ + · refine (lt_of_lt_of_le ?_ (single_le_sum (fun _ _ ↦ zero_le) (mem_antidiagonal.mpr + (add_zero s) : (s, 0) ∈ antidiagonal s))) + rw [← singleton_prod_singleton, Measure.prod_prod] + refine ENNReal.mul_pos ?_ (geometricMeasure_singleton_zero_pos p₂).ne' + rw [geometricMeasure_singleton h0] + exact (ENNReal.ofReal_pos.mpr (geometricMeasure_pos h0 h1 s)).ne' + · refine (lt_of_lt_of_le ?_ (single_le_sum (fun _ _ ↦ zero_le) (mem_antidiagonal.mpr + (zero_add s) : (0, s) ∈ antidiagonal s))) + rw [← singleton_prod_singleton, Measure.prod_prod] + refine ENNReal.mul_pos (geometricMeasure_singleton_zero_pos p₁).ne' ?_ + rw [geometricMeasure_singleton h0] + exact (ENNReal.ofReal_pos.mpr (geometricMeasure_pos h0 h1 s)).ne' + +/-- The geometric laws, as a kernel. -/ +private noncomputable def geomKernel : Kernel unitInterval ℕ := + ⟨geometricMeasure, measurable_geometricMeasure⟩ + +/-- The geometric laws are probability measures. -/ +private instance : IsMarkovKernel geomKernel := ⟨fun _ ↦ isProbabilityMeasure_geometricMeasure⟩ + +/-- The geometric kernel maps `p` to `Geom(p)`. -/ +private lemma geomKernel_apply (p : unitInterval) : geomKernel p = geometricMeasure p := rfl + +/-- The laws `Geom(1 / m) ⊗ Geom(1 / m)`, as a kernel (the success probability is clipped to +`[0, 1]`). -/ +private noncomputable def geomPairInvKernel : Kernel ℝ (ℕ × ℕ) := + (geomKernel ×ₖ geomKernel).comap (fun m : ℝ ↦ projIcc 0 1 zero_le_one m⁻¹) (by fun_prop) + +/-- The kernel `geomPairInvKernel` maps `m` to `Geom(1 / m) ⊗ Geom(1 / m)`, with the success +probability clipped to `[0, 1]`. -/ +private lemma geomPairInvKernel_apply (m : ℝ) : + geomPairInvKernel m = (geometricMeasure (projIcc 0 1 zero_le_one m⁻¹)).prod + (geometricMeasure (projIcc 0 1 zero_le_one m⁻¹)) := by + rw [geomPairInvKernel, Kernel.comap_apply, Kernel.prod_apply, geomKernel_apply] + +/-- The mass of a point under the mixture of the laws `Geom(1 / m) ⊗ Geom(1 / m)` over `m` +uniform between two distinct `a, b ≥ 1`. -/ +private lemma geomPairInvKernel_comp_singleton {a b : ℝ} (ha : 1 ≤ a) (hb : 1 ≤ b) (hab : a ≠ b) + (x : ℕ × ℕ) : + (geomPairInvKernel ∘ₘ uniformMeasure (min a b) (max a b)) {x} = + ENNReal.ofReal ((b - a)⁻¹ * (((1 - b⁻¹) ^ (x.1 + x.2 + 1) - (1 - a⁻¹) ^ (x.1 + x.2 + 1)) / + (((x.1 + x.2 : ℕ) : ℝ) + 1))) := by + set s := x.1 + x.2 + have hone : ∀ m ∈ uIcc a b, 1 ≤ m := fun m hm ↦ (le_inf ha hb).trans hm.1 + have hderiv : ∀ m ∈ uIcc a b, HasDerivAt (fun m : ℝ ↦ (1 - m⁻¹) ^ (s + 1) / (s + 1)) + ((1 - m⁻¹) ^ s * (m ^ 2)⁻¹) m := fun m hm ↦ + ((((hasDerivAt_inv (zero_lt_one.trans_le (hone m hm)).ne').const_sub 1).pow (s + 1)).div_const + _).congr_deriv (by field_simp; push_cast; ring) + have hcont : ContinuousOn (fun m : ℝ ↦ (1 - m⁻¹) ^ s * (m ^ 2)⁻¹) (uIcc a b) := by + refine ContinuousOn.mul + (by fun_prop (disch := exact fun m hm ↦ (zero_lt_one.trans_le (hone m hm)).ne')) ?_ + fun_prop (disch := exact fun m hm ↦ pow_ne_zero 2 (zero_lt_one.trans_le (hone m hm)).ne') + have hmass : ∀ᵐ m ∂uniformMeasure (min a b) (max a b), geomPairInvKernel m {x} = + ENNReal.ofReal ((1 - m⁻¹) ^ s * (m ^ 2)⁻¹) := by + filter_upwards [ae_mem_Ioc_uniformMeasure] with m hm + have hm1 : 1 ≤ m := (le_min ha hb).trans hm.1.le + have hinv : m⁻¹ ∈ Icc (0 : ℝ) 1 := ⟨by positivity, inv_le_one_of_one_le₀ hm1⟩ + have hne : projIcc (0 : ℝ) 1 zero_le_one m⁻¹ ≠ 0 := by + rw [ne_eq, ← Subtype.coe_inj, projIcc_of_mem _ hinv] + exact (inv_pos.mpr (zero_lt_one.trans_le hm1)).ne' + rw [geomPairInvKernel_apply, prod_geometricMeasure_singleton hne hne, projIcc_of_mem _ hinv] + congr 1 + simp only [s] + ring + rw [Measure.bind_apply (measurableSet_singleton x) (Kernel.aemeasurable _), + lintegral_congr_ae hmass, + lintegral_ofReal_uniformMeasure_min_max_of_hasDerivAt hab hderiv hcont fun m hm ↦ ?_, sub_div] + have : 0 ≤ 1 - m⁻¹ := sub_nonneg.mpr (inv_le_one_of_one_le₀ (hone m (uIoc_subset_uIcc hm))) + positivity + +/-- **The conditional law given the sum is a mixture of members of the null**: for distinct +positive success probabilities, the law with the law of the sum under `Geom(p₁) ⊗ Geom(p₂)` and +the conditional law `Kernel.antidiagonalUniform` given the sum is the mixture of the laws +`Geom(1 / m) ⊗ Geom(1 / m)` over `m` uniform between `1 / p₁` and `1 / p₂`. -/ +private lemma antidiagonalUniform_comp_map_add_prod_geometricMeasure (h₁ : p₁ ≠ 0) + (h₂ : p₂ ≠ 0) (h : p₁ ≠ p₂) : + Kernel.antidiagonalUniform ∘ₘ + ((geometricMeasure p₁).prod (geometricMeasure p₂)).map (fun x ↦ x.1 + x.2) = + geomPairInvKernel ∘ₘ uniformMeasure (min (p₁ : ℝ)⁻¹ (p₂ : ℝ)⁻¹) + (max (p₁ : ℝ)⁻¹ (p₂ : ℝ)⁻¹) := by + have h₁' : (0 : ℝ) < p₁ := unitInterval.pos_iff_ne_zero.mpr h₁ + have h₂' : (0 : ℝ) < p₂ := unitInterval.pos_iff_ne_zero.mpr h₂ + have hne : (p₁ : ℝ) ≠ p₂ := Subtype.coe_injective.ne h + have h12 : (p₁ : ℝ) - p₂ ≠ 0 := sub_ne_zero.mpr hne + have h21 : (p₂ : ℝ) - p₁ ≠ 0 := sub_ne_zero.mpr hne.symm + refine Measure.ext_iff_singleton.mpr fun x ↦ ?_ + set s := x.1 + x.2 + have hc : ((s : ℝ≥0∞) + 1)⁻¹ = ENNReal.ofReal ((s : ℝ) + 1)⁻¹ := by + rw [ENNReal.ofReal_inv_of_pos (by positivity), ENNReal.ofReal_add (by positivity) zero_le_one, + ENNReal.ofReal_natCast, ENNReal.ofReal_one] + have hnn : ∀ y ∈ antidiagonal s, + 0 ≤ (1 - (p₁ : ℝ)) ^ y.1 * p₁ * ((1 - p₂) ^ y.2 * p₂) := fun y _ ↦ + mul_nonneg (geometricMeasure_nonneg _ _) (geometricMeasure_nonneg _ _) + -- the law of the sum at `s`, a geometric sum over the antidiagonal + have hsum : ∑ y ∈ antidiagonal s, (1 - (p₁ : ℝ)) ^ y.1 * p₁ * ((1 - p₂) ^ y.2 * p₂) = + p₁ * p₂ * ((1 - p₁) ^ (s + 1) - (1 - p₂) ^ (s + 1)) / (p₂ - p₁) := by + rw [eq_div_iff h21, ← Nat.sum_antidiagonal_pow_mul_pow_mul_sub, sum_mul, sum_mul, mul_sum] + refine sum_congr rfl fun y _ ↦ ?_ + ring + rw [Kernel.antidiagonalUniform_comp_singleton, Measure.map_add_singleton_eq_sum_antidiagonal, + sum_congr rfl fun y _ ↦ prod_geometricMeasure_singleton h₁ h₂ y, + ← ENNReal.ofReal_sum_of_nonneg hnn, hc, ← ENNReal.ofReal_mul (sum_nonneg hnn), hsum, + geomPairInvKernel_comp_singleton (one_le_inv₀ h₁' |>.mpr p₁.2.2) + (one_le_inv₀ h₂' |>.mpr p₂.2.2) (inv_injective.ne hne), inv_inv, inv_inv] + congr 1 + field_simp + ring + +/-- The law with the law of the sum under `Geom(p₁) ⊗ Geom(p₂)` and the conditional law +`Kernel.antidiagonalUniform` given the sum is in the effective null. -/ +private lemma antidiagonalUniform_comp_mem_effectiveSet (p₁ p₂ : unitInterval) : + Kernel.antidiagonalUniform ∘ₘ + ((geometricMeasure p₁).prod (geometricMeasure p₂)).map (fun x ↦ x.1 + x.2) ∈ + effectiveSet (Set.range fun p : unitInterval ↦ + (geometricMeasure p).prod (geometricMeasure p)) := by + have hS : ∀ P ∈ Set.range fun p : unitInterval ↦ (geometricMeasure p).prod (geometricMeasure p), + P univ ≤ 1 := by + rintro _ ⟨p, rfl⟩ + exact prob_le_one + obtain ⟨q₁, hq₁, hq₁p⟩ := exists_ne_zero_geometricMeasure_eq p₁ + obtain ⟨q₂, hq₂, hq₂p⟩ := exists_ne_zero_geometricMeasure_eq p₂ + rw [← hq₁p, ← hq₂p] + rcases eq_or_ne q₁ q₂ with rfl | h + · obtain ⟨g, hg⟩ := exists_prod_geometricMeasure_self_singleton q₁ + rw [Kernel.antidiagonalUniform_comp_map_add_of_singleton hg] + exact subset_effectiveSet hS ⟨q₁, rfl⟩ + have := isProbabilityMeasure_uniformMeasure + (min_lt_max.mpr (inv_injective.ne (Subtype.coe_injective.ne h)) : + min (q₁ : ℝ)⁻¹ (q₂ : ℝ)⁻¹ < max (q₁ : ℝ)⁻¹ (q₂ : ℝ)⁻¹) + rw [antidiagonalUniform_comp_map_add_prod_geometricMeasure hq₁ hq₂ h] + refine comp_mem_effectiveSet prob_le_one (ae_of_all _ fun m ↦ subset_effectiveSet hS ?_) + exact ⟨projIcc 0 1 zero_le_one m⁻¹, (geomPairInvKernel_apply m).symm⟩ + +end Geometric + +/-! ### The numeraire -/ + +/-- **The numeraire of two geometric observations against a common success probability.** If one +of the laws `Geom(p₁)` and `Geom(p₂)` is not a point mass, the likelihood ratio of +`Geom(p₁) ⊗ Geom(p₂)` to the law under which the sum of the two observations has its law under +`Geom(p₁) ⊗ Geom(p₂)` and, given the sum `s`, the first observation is uniform on `{0, …, s}`, is +the numeraire of the null of the products `Geom(p) ⊗ Geom(p)`: it is an e-variable, and no +e-variable has a larger expected logarithm under `Geom(p₁) ⊗ Geom(p₂)` +(`ProbabilityTheory.IsNumeraire.eintegral_log_le`). This is the geometric instance of +[EVR26, Note 28]. -/ +theorem isNumeraire_rnDeriv_antidiagonalUniform_comp {p₁ p₂ : unitInterval} + (h : (p₁ ≠ 0 ∧ p₁ ≠ 1) ∨ (p₂ ≠ 0 ∧ p₂ ≠ 1)) : + IsNumeraire (((geometricMeasure p₁).prod (geometricMeasure p₂)).rnDeriv + (Kernel.antidiagonalUniform ∘ₘ + ((geometricMeasure p₁).prod (geometricMeasure p₂)).map fun x ↦ x.1 + x.2)) + (Set.range fun p : unitInterval ↦ (geometricMeasure p).prod (geometricMeasure p)) + ((geometricMeasure p₁).prod (geometricMeasure p₂)) := by + refine isNumeraire_of_ae_eq_rnDeriv (isEVar_rnDeriv_comp_map (measurable_of_countable _) + (ae_of_all _ Kernel.ae_add_antidiagonalUniform) + (Kernel.absolutelyContinuous_antidiagonalUniform_comp _) (by rintro _ ⟨p, rfl⟩; infer_instance) + ?_ ?_) (antidiagonalUniform_comp_mem_effectiveSet p₁ p₂) + (Kernel.absolutelyContinuous_antidiagonalUniform_comp _) .rfl + · rintro _ ⟨p, rfl⟩ + obtain ⟨g, hg⟩ := exists_prod_geometricMeasure_self_singleton p + exact Kernel.antidiagonalUniform_comp_map_add_of_singleton hg + · rintro _ ⟨p, rfl⟩ + exact Measure.absolutelyContinuous_of_forall_singleton_pos + (map_add_prod_geometricMeasure_singleton_pos h) + +end ProbabilityTheory diff --git a/LMLExplore/MeasureTheory/Measure/Singletons.lean b/LMLExplore/MeasureTheory/Measure/Singletons.lean new file mode 100644 index 00000000..e04503a3 --- /dev/null +++ b/LMLExplore/MeasureTheory/Measure/Singletons.lean @@ -0,0 +1,52 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import Mathlib.MeasureTheory.Measure.AbsolutelyContinuous +public import Mathlib.MeasureTheory.Measure.Map +public import Mathlib.Algebra.Order.Antidiag.Prod + +/-! +# Measures through their masses on points + +Two facts about measures read off their masses on points: a measure that charges every point +dominates every measure, and on a countable additive monoid with an antidiagonal (such as `ℕ`), the +law of the sum of the coordinates of a measure on pairs gives a point `s` the total mass of the +pairs with sum `s`. + +## Main statements + +* `MeasureTheory.Measure.absolutelyContinuous_of_forall_singleton_pos`: every measure is absolutely + continuous with respect to a measure that charges every point. +* `MeasureTheory.Measure.map_add_singleton_eq_sum_antidiagonal`: the mass of `{s}` under the law of + `x.1 + x.2` is `∑ x ∈ antidiagonal s, μ {x}`. +-/ + +public section + +open Set Finset + +namespace MeasureTheory.Measure + +/-- Every measure is absolutely continuous with respect to a measure that charges every point. -/ +theorem absolutelyContinuous_of_forall_singleton_pos {α : Type*} {_ : MeasurableSpace α} + {μ ν : Measure α} (h : ∀ x, 0 < ν {x}) : μ ≪ ν := + AbsolutelyContinuous.mk fun t _ ht ↦ by + rcases t.eq_empty_or_nonempty with rfl | ⟨x, hx⟩ + · exact measure_empty + · exact absurd ht ((h x).trans_le (measure_mono (singleton_subset_iff.mpr hx))).ne' + +/-- On a countable additive monoid with an antidiagonal, the mass of a point `s` under the law of +the sum of the coordinates is the sum of the masses of the pairs with sum `s`. -/ +theorem map_add_singleton_eq_sum_antidiagonal {M : Type*} [AddMonoid M] [HasAntidiagonal M] + [MeasurableSpace M] [Countable M] [MeasurableSingletonClass M] (μ : Measure (M × M)) (s : M) : + μ.map (fun x ↦ x.1 + x.2) {s} = ∑ x ∈ antidiagonal s, μ {x} := by + rw [map_apply (measurable_of_countable _) (measurableSet_singleton s), sum_measure_singleton] + congr 1 + ext x + simp + +end MeasureTheory.Measure diff --git a/LMLExplore/Probability/Distributions/Exponential.lean b/LMLExplore/Probability/Distributions/Exponential.lean new file mode 100644 index 00000000..478b8b44 --- /dev/null +++ b/LMLExplore/Probability/Distributions/Exponential.lean @@ -0,0 +1,173 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import Mathlib.Probability.Distributions.Exponential +public import TauCeti.Probability.Distributions.Hypoexponential.Basic +public import TauCeti.Probability.Distributions.Uniform.Basic +import LMLExplore.Probability.Distributions.Gamma +import LMLExplore.Probability.Distributions.Uniform +import TauCeti.MeasureTheory.Measure.Mixture.Basic +import TauCeti.Probability.Distributions.Beta.Uniform +import TauCeti.Probability.Distributions.Gamma.Beta +import TauCeti.Probability.Distributions.Gamma.Measurability +import TauCeti.Probability.Distributions.Hypoexponential.PDF + +/-! +# Exponential laws and the sum of two independent exponential variables + +Facts about Mathlib's exponential laws `expMeasure l`: they are finite measures, and those with +positive rates are equivalent. The sum of two independent exponential variables of rates `l₁` and +`l₂` has Tau Ceti's hypoexponential law, which is `Gamma(2, l)` when `l₁ = l₂ = l` and has the +density `l₁ l₂ (exp (- l₁ x) - exp (- l₂ x)) / (l₂ - l₁)` on `(0, ∞)` when `l₁ ≠ l₂`. In the +second case it is the mixture of the laws `Gamma(2, 1 / m)` over the mean `m` uniform between +`1 / l₁` and `1 / l₂`. + +## Main statements + +* `ProbabilityTheory.isFiniteMeasure_expMeasure`: exponential laws are finite measures. +* `ProbabilityTheory.expMeasure_absolutelyContinuous`, + `ProbabilityTheory.prod_expMeasure_absolutelyContinuous`: exponential laws, and their products, + with positive rates are equivalent. +* `ProbabilityTheory.map_add_prod_expMeasure_self`: `Exp(l) + Exp(l)` has the law `Gamma(2, l)`. +* `ProbabilityTheory.map_add_prod_expMeasure_eq_hypoexponentialMeasure`: `Exp(l₁) + Exp(l₂)` has + the hypoexponential law of rates `l₁, l₂`. +* `ProbabilityTheory.map_add_prod_expMeasure_eq_withDensity`: its density for `l₁ ≠ l₂`. +* `ProbabilityTheory.map_add_prod_expMeasure_eq_bind`: for `l₁ ≠ l₂`, it is the mixture of the + laws `Gamma(2, 1 / m)` over `m` uniform between `1 / l₁` and `1 / l₂`. +-/ + +public section + +open MeasureTheory Set Real TauCeti.Probability + +namespace ProbabilityTheory + +variable {l l₁ l₂ : ℝ} + +/-- Exponential laws are finite measures: probability measures for a positive rate, and the zero +measure otherwise. -/ +instance isFiniteMeasure_expMeasure (l : ℝ) : IsFiniteMeasure (expMeasure l) := by + rcases le_or_gt l 0 with hl | hl + · rw [expMeasure_of_nonpos hl] + infer_instance + · have := isProbabilityMeasure_expMeasure hl + infer_instance + +/-- Exponential laws with positive rates are equivalent. -/ +theorem expMeasure_absolutelyContinuous (l : ℝ) {l' : ℝ} (hl' : 0 < l') : + expMeasure l ≪ expMeasure l' := by + rw [expMeasure_eq_withDensity_restrict_Ioi l, expMeasure_eq_withDensity_restrict_Ioi l'] + refine (withDensity_absolutelyContinuous _ _).trans (withDensity_absolutelyContinuous' + (measurable_exponentialPDFReal l').ennreal_ofReal.aemeasurable ?_) + filter_upwards [ae_restrict_mem measurableSet_Ioi] with x hx + exact (ENNReal.ofReal_pos.mpr (exponentialPDFReal_pos hl' hx)).ne' + +/-- Products of exponential laws with positive rates are equivalent. -/ +theorem prod_expMeasure_absolutelyContinuous (l₁ l₂ : ℝ) {l₁' l₂' : ℝ} (h₁ : 0 < l₁') + (h₂ : 0 < l₂') : + (expMeasure l₁).prod (expMeasure l₂) ≪ (expMeasure l₁').prod (expMeasure l₂') := + (expMeasure_absolutelyContinuous l₁ h₁).prod (expMeasure_absolutelyContinuous l₂ h₂) + +/-- The sum of two independent `Exp(l)` variables has the law `Gamma(2, l)`. -/ +theorem map_add_prod_expMeasure_self (hl : 0 < l) : + ((expMeasure l).prod (expMeasure l)).map (fun x ↦ x.1 + x.2) = gammaMeasure 2 l := by + have : IsProbabilityMeasure (betaMeasure 1 1) := by rw [betaMeasure_one_one]; infer_instance + have h := congrArg (Measure.map Prod.snd) (map_div_add_prod_gammaMeasure one_pos one_pos hl) + rw [Measure.map_map measurable_snd (by fun_prop), Measure.map_snd_prod, measure_univ, + one_smul, one_add_one_eq_two] at h + exact h + +/-- **The sum of two independent exponential variables is hypoexponential**: for positive rates +`l₁, l₂`, the law of the sum under `Exp(l₁) ⊗ Exp(l₂)` is `hypoexponentialMeasure ![l₁, l₂]`. -/ +theorem map_add_prod_expMeasure_eq_hypoexponentialMeasure (h₁ : 0 < l₁) (h₂ : 0 < l₂) : + ((expMeasure l₁).prod (expMeasure l₂)).map (fun x ↦ x.1 + x.2) = + hypoexponentialMeasure ![l₁, l₂] := by + have hr : ∀ i, 0 < ![l₁, l₂] i := by simp [Fin.forall_fin_two, h₁, h₂] + have hpi : (fun i ↦ expMeasure (![l₁, l₂] i)) = ![expMeasure l₁, expMeasure l₂] := by + ext i : 1 + fin_cases i <;> rfl + rw [hypoexponentialMeasure_eq_map hr, hpi, + ← (measurePreserving_finTwoArrow_vec (expMeasure l₁) (expMeasure l₂)).map_eq, + Measure.map_map (by fun_prop) (MeasurableEquiv.measurable _)] + congr 1 + ext x + simp [Fin.sum_univ_two, MeasurableEquiv.finTwoArrow] + +/-- **The density of the sum of two exponential variables** with distinct positive rates +`l₁, l₂` is `l₁ l₂ (exp (- l₁ x) - exp (- l₂ x)) / (l₂ - l₁)` on `(0, ∞)`. -/ +theorem map_add_prod_expMeasure_eq_withDensity (h₁ : 0 < l₁) (h₂ : 0 < l₂) (h : l₁ ≠ l₂) : + ((expMeasure l₁).prod (expMeasure l₂)).map (fun x ↦ x.1 + x.2) = + (volume.restrict (Ioi 0)).withDensity fun x ↦ + ENNReal.ofReal (l₁ * l₂ * (exp (-(l₁ * x)) - exp (-(l₂ * x))) / (l₂ - l₁)) := by + have hr : ∀ i, 0 < ![l₁, l₂] i := by simp [Fin.forall_fin_two, h₁, h₂] + have hinj : Function.Injective ![l₁, l₂] := by + intro i j hij + fin_cases i <;> fin_cases j <;> simp_all [eq_comm] + rw [map_add_prod_expMeasure_eq_hypoexponentialMeasure h₁ h₂, + ← Measure.restrict_eq_self_of_ae_mem (s := Ioi 0) + (μ := hypoexponentialMeasure ![l₁, l₂]) (ae_pos_hypoexponentialMeasure _), + hypoexponentialMeasure_eq_withDensity, restrict_withDensity measurableSet_Ioi] + refine withDensity_congr_ae ((ae_restrict_iff' measurableSet_Ioi).mpr + (ae_of_all _ fun x (hx : 0 < x) ↦ ?_)) + rw [hypoexponentialPDF_eq_ofReal, hypoexponentialPDFReal_of_injective hr hinj hx] + have h01 : (Finset.univ.erase (0 : Fin 2)) = {1} := by decide + have h10 : (Finset.univ.erase (1 : Fin 2)) = {0} := by decide + have h12 : l₂ - l₁ ≠ 0 := sub_ne_zero.mpr h.symm + have h21 : l₁ - l₂ ≠ 0 := sub_ne_zero.mpr h + simp only [Fin.sum_univ_two, h01, h10, Finset.prod_singleton, Matrix.cons_val_zero, + Matrix.cons_val_one] + congr 1 + field_simp + ring + +/-- The density of the mixture of the laws `Gamma(2, 1 / m)` over the mean `m` uniform between +`1 / l₁` and `1 / l₂`, at `x > 0`. -/ +private lemma lintegral_gammaPDF_uniformMeasure (h₁ : 0 < l₁) (h₂ : 0 < l₂) (h : l₁ ≠ l₂) + {x : ℝ} (hx : 0 < x) : + ∫⁻ m, gammaPDF 2 m⁻¹ x ∂uniformMeasure (min l₁⁻¹ l₂⁻¹) (max l₁⁻¹ l₂⁻¹) = + ENNReal.ofReal (l₁ * l₂ * (exp (-(l₁ * x)) - exp (-(l₂ * x))) / (l₂ - l₁)) := by + have hpos : ∀ m ∈ uIcc l₁⁻¹ l₂⁻¹, 0 < m := fun m hm ↦ + (lt_inf_iff.mpr ⟨inv_pos.mpr h₁, inv_pos.mpr h₂⟩).trans_le hm.1 + have hderiv : ∀ m ∈ uIcc l₁⁻¹ l₂⁻¹, HasDerivAt (fun m ↦ exp (-(m⁻¹ * x))) + (exp (-(m⁻¹ * x)) * ((m ^ 2)⁻¹ * x)) m := + fun m hm ↦ (((hasDerivAt_inv (hpos m hm).ne').mul_const x).neg.exp).congr_deriv + (by simp only [Pi.neg_apply]; ring) + have hcont : ContinuousOn (fun m ↦ exp (-(m⁻¹ * x)) * ((m ^ 2)⁻¹ * x)) (uIcc l₁⁻¹ l₂⁻¹) := by + refine ContinuousOn.mul (by fun_prop (disch := exact fun m hm ↦ (hpos m hm).ne')) ?_ + fun_prop (disch := exact fun m hm ↦ pow_ne_zero 2 (hpos m hm).ne') + have hpdf : ∀ᵐ m ∂uniformMeasure (min l₁⁻¹ l₂⁻¹) (max l₁⁻¹ l₂⁻¹), gammaPDF 2 m⁻¹ x = + ENNReal.ofReal (exp (-(m⁻¹ * x)) * ((m ^ 2)⁻¹ * x)) := by + filter_upwards [ae_mem_Ioc_uniformMeasure] with m hm + rw [gammaPDF_two hx.le] + congr 1 + ring + rw [lintegral_congr_ae hpdf, lintegral_ofReal_uniformMeasure_min_max_of_hasDerivAt + (inv_injective.ne h) hderiv hcont fun m hm ↦ ?_] + · have h12 : l₂ - l₁ ≠ 0 := sub_ne_zero.mpr h.symm + have h21 : l₁ - l₂ ≠ 0 := sub_ne_zero.mpr h + simp only [inv_inv] + congr 1 + field_simp + ring + · have := hpos m (uIoc_subset_uIcc hm) + positivity + +/-- **The sum of two exponential variables is a mixture of gamma laws**: for distinct positive +rates, the law of the sum is the mixture of the laws `Gamma(2, 1 / m)` over the mean `m` uniform +between `1 / l₁` and `1 / l₂`. -/ +theorem map_add_prod_expMeasure_eq_bind (h₁ : 0 < l₁) (h₂ : 0 < l₂) (h : l₁ ≠ l₂) : + ((expMeasure l₁).prod (expMeasure l₂)).map (fun x ↦ x.1 + x.2) = + (uniformMeasure (min l₁⁻¹ l₂⁻¹) (max l₁⁻¹ l₂⁻¹)).bind fun m ↦ gammaMeasure 2 m⁻¹ := by + simp_rw [gammaMeasure_eq_withDensity_restrict_Ioi] + have hf : Measurable (Function.uncurry fun m x ↦ gammaPDF 2 m⁻¹ x) := + measurable_uncurry_gammaPDF.comp (f := fun p : ℝ × ℝ ↦ (((2 : ℝ), p.1⁻¹), p.2)) (by fun_prop) + rw [map_add_prod_expMeasure_eq_withDensity h₁ h₂ h, Measure.bind_withDensity _ _ hf] + refine withDensity_congr_ae ((ae_restrict_iff' measurableSet_Ioi).mpr + (ae_of_all _ fun x (hx : 0 < x) ↦ ?_)) + exact (lintegral_gammaPDF_uniformMeasure h₁ h₂ h hx).symm + +end ProbabilityTheory diff --git a/LMLExplore/Probability/Distributions/Gamma.lean b/LMLExplore/Probability/Distributions/Gamma.lean new file mode 100644 index 00000000..c9d3868e --- /dev/null +++ b/LMLExplore/Probability/Distributions/Gamma.lean @@ -0,0 +1,33 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import Mathlib.Probability.Distributions.Gamma +import Mathlib.Analysis.SpecialFunctions.Gamma.BohrMollerup + +/-! +# The density of the gamma law of shape two + +The density of Mathlib's gamma law `gammaMeasure 2 r` at `x ≥ 0` is `r ^ 2 * x * exp (- r x)`. + +## Main statements + +* `ProbabilityTheory.gammaPDF_two`: the density of `Gamma(2, r)` on `[0, ∞)`. +-/ + +public section + +open Real + +namespace ProbabilityTheory + +/-- The density of `Gamma(2, r)` on `[0, ∞)` is `r ^ 2 * x * exp (- r x)`. -/ +theorem gammaPDF_two {r x : ℝ} (hx : 0 ≤ x) : + gammaPDF 2 r x = ENNReal.ofReal (r ^ 2 * x * exp (-(r * x))) := by + rw [gammaPDF_of_nonneg hx, Real.Gamma_two, show (2 : ℝ) - 1 = 1 by norm_num, Real.rpow_one, + Real.rpow_two, div_one] + +end ProbabilityTheory diff --git a/LMLExplore/Probability/Distributions/Geometric.lean b/LMLExplore/Probability/Distributions/Geometric.lean new file mode 100644 index 00000000..f80079e3 --- /dev/null +++ b/LMLExplore/Probability/Distributions/Geometric.lean @@ -0,0 +1,57 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import Mathlib.Probability.Distributions.Geometric +import TauCeti.Probability.Distributions.Geometric.Basic + +/-! +# Geometric laws with a positive success probability + +Mathlib's `geometricMeasure 0` is the point mass at `0`, as is `geometricMeasure 1` +(Tau Ceti's `TauCeti.Probability.geometricMeasure_zero` and +`TauCeti.Probability.geometricMeasure_one`), so every geometric law is one with a positive success +probability, for which Mathlib's `geometricMeasure_singleton` gives the masses of the points. + +## Main statements + +* `ProbabilityTheory.exists_ne_zero_geometricMeasure_eq`: every geometric law is one with a + positive success probability. +* `ProbabilityTheory.geometricMeasure_singleton_zero_pos`: every geometric law charges `0`. +* `ProbabilityTheory.prod_geometricMeasure_singleton`: the mass of a point under a product of two + geometric laws with positive success probabilities. +-/ + +public section + +open MeasureTheory Set + +namespace ProbabilityTheory + +/-- Every geometric law is one with a positive success probability. -/ +theorem exists_ne_zero_geometricMeasure_eq (p : unitInterval) : + ∃ q : unitInterval, q ≠ 0 ∧ geometricMeasure q = geometricMeasure p := by + by_cases hp : p = 0 + · exact ⟨1, one_ne_zero, by rw [hp, TauCeti.Probability.geometricMeasure_one, + TauCeti.Probability.geometricMeasure_zero]⟩ + · exact ⟨p, hp, rfl⟩ + +/-- The point `0` has positive mass under every geometric law. -/ +theorem geometricMeasure_singleton_zero_pos (p : unitInterval) : 0 < geometricMeasure p {0} := by + obtain ⟨q, hq, hqp⟩ := exists_ne_zero_geometricMeasure_eq p + rw [← hqp, geometricMeasure_singleton hq] + simpa using unitInterval.pos_iff_ne_zero.mpr hq + +/-- The mass of a point under a product of geometric laws with positive success +probabilities. -/ +theorem prod_geometricMeasure_singleton {p₁ p₂ : unitInterval} (h₁ : p₁ ≠ 0) (h₂ : p₂ ≠ 0) + (x : ℕ × ℕ) : + (geometricMeasure p₁).prod (geometricMeasure p₂) {x} = + ENNReal.ofReal ((1 - p₁) ^ x.1 * p₁ * ((1 - p₂) ^ x.2 * p₂)) := by + rw [← singleton_prod_singleton, Measure.prod_prod, geometricMeasure_singleton h₁, + geometricMeasure_singleton h₂, ← ENNReal.ofReal_mul (geometricMeasure_nonneg _ _)] + +end ProbabilityTheory diff --git a/LMLExplore/Probability/Distributions/Uniform.lean b/LMLExplore/Probability/Distributions/Uniform.lean new file mode 100644 index 00000000..d9edf34c --- /dev/null +++ b/LMLExplore/Probability/Distributions/Uniform.lean @@ -0,0 +1,73 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import Mathlib.MeasureTheory.Integral.IntervalIntegral.FundThmCalculus +public import TauCeti.Probability.Distributions.Uniform.Basic + +/-! +# Integrating a derivative against the uniform law + +Against Tau Ceti's uniform law `uniformMeasure a b` on `(a, b]`, the lower Lebesgue integral of a +nonnegative continuous derivative `f'` of `f` is the mean slope `(b - a)⁻¹ * (f b - f a)`, by the +fundamental theorem of calculus. The second form takes the uniform law between two distinct +points `a` and `b` in either order, `uniformMeasure (min a b) (max a b)`; the mean slope is +symmetric in `a` and `b`. + +## Main statements + +* `ProbabilityTheory.lintegral_ofReal_uniformMeasure_of_hasDerivAt`: for `a < b`, + `∫⁻ x, ENNReal.ofReal (f' x) ∂uniformMeasure a b = ENNReal.ofReal ((b - a)⁻¹ * (f b - f a))`. +* `ProbabilityTheory.lintegral_ofReal_uniformMeasure_min_max_of_hasDerivAt`: the same for the + uniform law between `a ≠ b`. +-/ + +public section + +open MeasureTheory Set TauCeti.Probability + +namespace ProbabilityTheory + +variable {a b : ℝ} {f f' : ℝ → ℝ} + +/-- **The integral of a derivative against the uniform law is the mean slope.** For `a < b` and a +function `f` with a continuous derivative `f'` on `[a, b]`, nonnegative on `(a, b]`, the integral +of `f'` against the uniform law on `(a, b]` is `(b - a)⁻¹ * (f b - f a)`. -/ +theorem lintegral_ofReal_uniformMeasure_of_hasDerivAt (hab : a < b) + (hf : ∀ x ∈ Icc a b, HasDerivAt f (f' x) x) (hf' : ContinuousOn f' (Icc a b)) + (hf'0 : ∀ x ∈ Ioc a b, 0 ≤ f' x) : + ∫⁻ x, ENNReal.ofReal (f' x) ∂uniformMeasure a b = + ENNReal.ofReal ((b - a)⁻¹ * (f b - f a)) := by + rw [← uIcc_of_le hab.le] at hf hf' + rw [uniformMeasure_eq_smul, lintegral_smul_measure, ← ofReal_integral_eq_lintegral_ofReal, + ← intervalIntegral.integral_of_le hab.le, + intervalIntegral.integral_eq_sub_of_hasDerivAt hf hf'.intervalIntegrable, smul_eq_mul, + ← ENNReal.ofReal_inv_of_pos (sub_pos.mpr hab), + ← ENNReal.ofReal_mul (inv_nonneg.mpr (sub_pos.mpr hab).le)] + · rw [uIcc_of_le hab.le] at hf' + exact hf'.integrableOn_Icc.mono_set Ioc_subset_Icc_self + · exact ae_restrict_of_forall_mem measurableSet_Ioc hf'0 + +/-- **The integral of a derivative against the uniform law between two points.** For `a ≠ b` and +a function `f` with a continuous derivative `f'` on `uIcc a b`, nonnegative on `uIoc a b`, the +integral of `f'` against the uniform law between `a` and `b` is `(b - a)⁻¹ * (f b - f a)`. -/ +theorem lintegral_ofReal_uniformMeasure_min_max_of_hasDerivAt (hab : a ≠ b) + (hf : ∀ x ∈ uIcc a b, HasDerivAt f (f' x) x) (hf' : ContinuousOn f' (uIcc a b)) + (hf'0 : ∀ x ∈ uIoc a b, 0 ≤ f' x) : + ∫⁻ x, ENNReal.ofReal (f' x) ∂uniformMeasure (min a b) (max a b) = + ENNReal.ofReal ((b - a)⁻¹ * (f b - f a)) := by + rcases hab.lt_or_gt with h | h + · rw [min_eq_left h.le, max_eq_right h.le] + rw [uIcc_of_le h.le] at hf hf' + rw [uIoc_of_le h.le] at hf'0 + exact lintegral_ofReal_uniformMeasure_of_hasDerivAt h hf hf' hf'0 + · rw [min_eq_right h.le, max_eq_left h.le, ← neg_sub a b, inv_neg, ← neg_sub (f a) (f b), + neg_mul_neg] + rw [uIcc_of_ge h.le] at hf hf' + rw [uIoc_of_ge h.le] at hf'0 + exact lintegral_ofReal_uniformMeasure_of_hasDerivAt h hf hf' hf'0 + +end ProbabilityTheory diff --git a/LMLExplore/Probability/Kernel/Composition/MeasureComp.lean b/LMLExplore/Probability/Kernel/Composition/MeasureComp.lean new file mode 100644 index 00000000..c901e496 --- /dev/null +++ b/LMLExplore/Probability/Kernel/Composition/MeasureComp.lean @@ -0,0 +1,50 @@ +/- +Copyright (c) 2026 LMLExplore contributors. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: LMLExplore contributors +-/ +module + +public import Mathlib.Probability.Kernel.Composition.MeasureComp + +/-! +# The law of a statistic under a composition with a kernel carried by its fibers + +If a Markov kernel `κ : Kernel 𝓨 Ω` is carried by the fibers of a statistic `T : Ω → 𝓨`, that is, +`κ y` is carried by `T ⁻¹' {y}`, then the law of `T` under `κ ∘ₘ ν` is `ν`: the composition with +`κ` draws a point of the fiber of each `y`. + +## Main statements + +* `ProbabilityTheory.map_comp_eq_of_ae_eq`: `(κ ∘ₘ ν).map T = ν`. +-/ + +public section + +open MeasureTheory Set + +namespace ProbabilityTheory + +variable {Ω 𝓨 : Type*} {mΩ : MeasurableSpace Ω} {m𝓨 : MeasurableSpace 𝓨} {T : Ω → 𝓨} + {κ : Kernel 𝓨 Ω} + +/-- If `κ y` is a probability measure carried by the fiber `T ⁻¹' {y}` for `ν`-almost every `y`, +the law of `T` under `κ ∘ₘ ν` is `ν`. -/ +theorem map_comp_eq_of_ae_eq [IsMarkovKernel κ] (hT : Measurable T) {ν : Measure 𝓨} + (hκ : ∀ᵐ y ∂ν, ∀ᵐ x ∂κ y, T x = y) : (κ ∘ₘ ν).map T = ν := by + ext A hA + rw [Measure.map_apply hT hA, Measure.bind_apply (hT hA) κ.aemeasurable, + ← lintegral_indicator_one hA] + refine lintegral_congr_ae ?_ + filter_upwards [hκ] with y hκ + by_cases hy : y ∈ A + · rw [indicator_of_mem hy, Pi.one_apply, ← measure_univ (μ := κ y)] + refine measure_congr ?_ + filter_upwards [hκ] with x hx + simp [hx, hy] + · rw [indicator_of_notMem hy] + refine measure_eq_zero_iff_ae_notMem.mpr ?_ + filter_upwards [hκ] with x hx + simpa [hx] using hy + +end ProbabilityTheory