diff --git a/Mathlib.lean b/Mathlib.lean index bb572a5a061e4f..cdd90a93d8d70a 100644 --- a/Mathlib.lean +++ b/Mathlib.lean @@ -6151,6 +6151,7 @@ public import Mathlib.Probability.Combinatorics.BinomialRandomGraph.Defs public import Mathlib.Probability.CondVar public import Mathlib.Probability.ConditionalExpectation public import Mathlib.Probability.ConditionalProbability +public import Mathlib.Probability.Decision.HypothesisTesting.NeymanPearson public import Mathlib.Probability.Decision.Risk.Basic public import Mathlib.Probability.Decision.Risk.Defs public import Mathlib.Probability.Density diff --git a/Mathlib/Probability/Decision/HypothesisTesting/NeymanPearson.lean b/Mathlib/Probability/Decision/HypothesisTesting/NeymanPearson.lean new file mode 100644 index 00000000000000..a45aabe12f63fa --- /dev/null +++ b/Mathlib/Probability/Decision/HypothesisTesting/NeymanPearson.lean @@ -0,0 +1,178 @@ +/- +Copyright (c) 2026 Zixiao Wang. All rights reserved. +Released under Apache 2.0 license as described in the file LICENSE. +Authors: Zixiao Wang, Rajarshi Mukherjee +-/ +module + +public import Mathlib.MeasureTheory.Measure.Decomposition.Lebesgue +public import Mathlib.MeasureTheory.Measure.Decomposition.RadonNikodym + +/-! +# The Neyman–Pearson lemma (sum-of-errors form) + +For two finite measures `P` and `Q` on a measurable space `Ω`, the Neyman–Pearson lemma +characterizes the randomized tests `T : Ω → ℝ` with values in `[0, 1]` that minimize +the sum of error masses +`∫ T dP + t · ∫ (1 - T) dQ` +where `t > 0` is a fixed cost weight on Type II errors (equivalently, the test compares +the likelihood ratio against the threshold `t⁻¹`). + +## Main definitions + +* `ProbabilityTheory.NeymanPearson.npLikelihoodRatio P Q` — the likelihood ratio + `(dQ/d(P+Q)) / (dP/d(P+Q))` of `Q` against `P`. +* `ProbabilityTheory.NeymanPearson.npTest P Q t η` — the Neyman–Pearson test: +`1` if the likelihood ratio is `> t⁻¹`, `η` if it equals `t⁻¹`, and `0` otherwise; + forced to `1` on `{dP/d(P+Q) = 0}`. +* `ProbabilityTheory.NeymanPearson.npSumError P Q t T` — the sum-of-errors cost + `∫ T dP + t · ∫ (1 - T) dQ`. + +## Main result + +* `ProbabilityTheory.NeymanPearson.neyman_pearson` — `npTest P Q t η` minimizes + `npSumError P Q t` over the set of measurable `[0, 1]`-valued tests, for any + `t > 0` and `0 ≤ η ≤ 1`. + +## References +* Ingster, Y. and Suslina, I.A., 2012. Nonparametric goodness-of-fit testing under Gaussian models (Vol. 169). Springer Science & Business Media. +-/ + +@[expose] public section + +open MeasureTheory +open scoped ENNReal + +namespace ProbabilityTheory.NeymanPearson + +variable {Ω : Type*} [MeasurableSpace Ω] (P Q : Measure Ω) + +/-- The likelihood ratio of `Q` against `P` with respect to the dominating measure +`P + Q`, i.e. `(dQ/d(P+Q)) / (dP/d(P+Q))`. On the `P`-null set where `dP/d(P+Q)` +vanishes, this evaluates to `0` by Lean's `a / 0 = 0` convention; the value there +is immaterial to the optimality statement `neyman_pearson` because `npTest` is +hard-coded to `1` on that set. -/ +noncomputable def npLikelihoodRatio : Ω → ℝ := fun ω ↦ + (Q.rnDeriv (P + Q) ω).toReal / (P.rnDeriv (P + Q) ω).toReal + +/-- The Neyman–Pearson randomized test against the threshold `t⁻¹`: returns `1` if the +likelihood ratio is `> t⁻¹`, `η` at the boundary `= t⁻¹`, and `0` otherwise; returns `1` +on the `P`-null set where the ratio is undefined. The test is optimal for +`npSumError P Q t` when `t > 0`. Boundedness in `[0, 1]` requires `0 ≤ η ≤ 1` +(see `npTest_nonneg` and `npTest_le_one`). -/ +noncomputable def npTest (t η : ℝ) : Ω → ℝ := fun ω ↦ + if (P.rnDeriv (P + Q) ω).toReal = 0 then 1 + else if npLikelihoodRatio P Q ω > t⁻¹ then 1 + else if npLikelihoodRatio P Q ω = t⁻¹ then η + else 0 + +/-- The sum-of-errors cost functional `∫ T dP + t · ∫ (1 - T) dQ`. -/ +noncomputable def npSumError (t : ℝ) (T : Ω → ℝ) : ℝ := + ∫ ω, T ω ∂P + t * ∫ ω, (1 - T ω) ∂Q + +/-- The likelihood ratio is measurable. -/ +@[fun_prop] +lemma measurable_npLikelihoodRatio : Measurable (npLikelihoodRatio P Q) := + (Measure.measurable_rnDeriv Q (P + Q)).ennreal_toReal.div + (Measure.measurable_rnDeriv P (P + Q)).ennreal_toReal + +/-- The Neyman–Pearson test is measurable. -/ +@[fun_prop] +lemma measurable_npTest (t η : ℝ) : Measurable (npTest P Q t η) := by + have hh := measurable_npLikelihoodRatio P Q + have hp : Measurable (fun ω ↦ (P.rnDeriv (P + Q) ω).toReal) := by fun_prop + exact Measurable.ite (measurableSet_eq_fun hp measurable_const) measurable_const <| + Measurable.ite (measurableSet_lt measurable_const hh) measurable_const <| + Measurable.ite (measurableSet_eq_fun hh measurable_const) measurable_const measurable_const + +/-- The Neyman–Pearson test is non-negative when `0 ≤ η`. -/ +lemma npTest_nonneg {t η : ℝ} (hη0 : 0 ≤ η) (ω : Ω) : 0 ≤ npTest P Q t η ω := by + simp only [npTest] + split_ifs <;> linarith + +/-- The Neyman–Pearson test is bounded by `1` when `η ≤ 1`. -/ +lemma npTest_le_one {t η : ℝ} (hη1 : η ≤ 1) (ω : Ω) : npTest P Q t η ω ≤ 1 := by + simp only [npTest] + split_ifs <;> linarith + +/-- **Neyman–Pearson lemma (sum-of-errors form).** For finite measures `P` and `Q` on +a measurable space `Ω`, the randomized test `npTest P Q t η` minimizes +`npSumError P Q t T = ∫ T dP + t · ∫ (1 - T) dQ` over all measurable `[0, 1]`-valued +tests `T`, for any threshold `t > 0` and any `0 ≤ η ≤ 1`. -/ +theorem neyman_pearson + [IsFiniteMeasure P] [IsFiniteMeasure Q] + (t η : ℝ) (ht : 0 < t) (hη : η ∈ Set.Icc (0 : ℝ) 1) : + IsMinOn (npSumError P Q t) + { T : Ω → ℝ | Measurable T ∧ (∀ ω, 0 ≤ T ω) ∧ (∀ ω, T ω ≤ 1) } + (npTest P Q t η) := by + intro T hT + obtain ⟨hTm, hT0, hT1⟩ := hT + obtain ⟨hη0, hη1⟩ := hη + set τ := npTest P Q t η + change npSumError P Q t τ ≤ npSumError P Q t T + have hPμ : P ≪ P + Q := Measure.AbsolutelyContinuous.add_right .rfl Q + have hQμ : Q ≪ P + Q := Measure.AbsolutelyContinuous.add_right' .rfl P + have hτm : Measurable τ := measurable_npTest P Q t η + have hτ0 : ∀ ω, 0 ≤ τ ω := npTest_nonneg P Q hη0 + have hτ1 : ∀ ω, τ ω ≤ 1 := npTest_le_one P Q hη1 + have h_int_μ : ∀ F : Ω → ℝ, Measurable F → (∀ ω, 0 ≤ F ω) → (∀ ω, F ω ≤ 1) + → Integrable (fun ω ↦ (P.rnDeriv (P + Q) ω).toReal * F ω) (P + Q) ∧ + Integrable (fun ω ↦ (Q.rnDeriv (P + Q) ω).toReal * (1 - F ω)) (P + Q) := by + intro F hFm hF0 hF1 + have hF_P : Integrable F P := .of_mem_Icc 0 1 hFm.aemeasurable + (ae_of_all _ fun ω ↦ Set.mem_Icc.mpr ⟨hF0 ω, hF1 ω⟩) + have h1F_m : Measurable (fun ω ↦ 1 - F ω) := by fun_prop + have h1F_Q : Integrable (fun ω ↦ 1 - F ω) Q := .of_mem_Icc 0 1 h1F_m.aemeasurable + (ae_of_all _ fun ω ↦ Set.mem_Icc.mpr ⟨by linarith [hF1 ω], by linarith [hF0 ω]⟩) + exact ⟨(integrable_toReal_rnDeriv_mul_iff hPμ).mpr hF_P, + (integrable_toReal_rnDeriv_mul_iff hQμ).mpr h1F_Q⟩ + obtain ⟨hpτ, hq1τ⟩ := h_int_μ τ hτm hτ0 hτ1 + obtain ⟨hpT, hq1T⟩ := h_int_μ T hTm hT0 hT1 + have h_change : ∀ F : Ω → ℝ, + Integrable (fun ω ↦ (P.rnDeriv (P + Q) ω).toReal * F ω) (P + Q) → + Integrable (fun ω ↦ (Q.rnDeriv (P + Q) ω).toReal * (1 - F ω)) (P + Q) → + npSumError P Q t F = + ∫ ω, (P.rnDeriv (P + Q) ω).toReal * F ω + + t * ((Q.rnDeriv (P + Q) ω).toReal * (1 - F ω)) ∂(P + Q) := fun F hpF hq1F ↦ by + change ∫ ω, F ω ∂P + t * ∫ ω, (1 - F ω) ∂Q = _ + rw [← integral_toReal_rnDeriv_mul hPμ, ← integral_toReal_rnDeriv_mul hQμ, + ← integral_const_mul, ← integral_add hpF (hq1F.const_mul t)] + rw [h_change τ hpτ hq1τ, h_change T hpT hq1T] + refine integral_mono (hpτ.add (hq1τ.const_mul t)) (hpT.add (hq1T.const_mul t)) ?_ + intro ω + dsimp only + by_cases hp_zero : (P.rnDeriv (P + Q) ω).toReal = 0 + · have hτ_one : τ ω = 1 := by + change npTest P Q t η ω = 1 + simp only [npTest, if_pos hp_zero] + rw [hp_zero, hτ_one] + nlinarith [mul_nonneg ht.le (mul_nonneg + (ENNReal.toReal_nonneg : (0:ℝ) ≤ (Q.rnDeriv (P + Q) ω).toReal) + (sub_nonneg.mpr (hT1 ω)))] + · have hp_pos : 0 < (P.rnDeriv (P + Q) ω).toReal := + lt_of_le_of_ne ENNReal.toReal_nonneg (Ne.symm hp_zero) + have h_ratio : npLikelihoodRatio P Q ω = + (Q.rnDeriv (P + Q) ω).toReal / (P.rnDeriv (P + Q) ω).toReal := rfl + rcases lt_trichotomy (npLikelihoodRatio P Q ω) t⁻¹ with h_lt | h_eq | h_gt + · have hτ_zero : τ ω = 0 := by + change npTest P Q t η ω = 0 + simp only [npTest, if_neg hp_zero, if_neg (not_lt.mpr h_lt.le), if_neg (ne_of_lt h_lt)] + rw [h_ratio, div_lt_iff₀ hp_pos] at h_lt + have htq_lt_p : t * (Q.rnDeriv (P + Q) ω).toReal < (P.rnDeriv (P + Q) ω).toReal := by + have := mul_lt_mul_of_pos_left h_lt ht + rwa [← mul_assoc, mul_inv_cancel₀ ht.ne', one_mul] at this + rw [hτ_zero]; nlinarith [hT0 ω, hT1 ω, htq_lt_p] + · rw [h_ratio, div_eq_iff hp_zero] at h_eq + have hp_eq : (P.rnDeriv (P + Q) ω).toReal = t * (Q.rnDeriv (P + Q) ω).toReal := by + rw [h_eq, ← mul_assoc, mul_inv_cancel₀ ht.ne', one_mul] + rw [hp_eq]; ring_nf; exact le_rfl + · have hτ_one : τ ω = 1 := by + change npTest P Q t η ω = 1 + simp only [npTest, if_neg hp_zero, if_pos h_gt] + rw [h_ratio, lt_div_iff₀ hp_pos] at h_gt + have hp_lt_tq : (P.rnDeriv (P + Q) ω).toReal < t * (Q.rnDeriv (P + Q) ω).toReal := by + have := mul_lt_mul_of_pos_left h_gt ht + rwa [← mul_assoc, mul_inv_cancel₀ ht.ne', one_mul] at this + rw [hτ_one]; nlinarith [hT0 ω, hT1 ω, hp_lt_tq] + +end ProbabilityTheory.NeymanPearson