prospicio_prob/
tweedie.rs

1//! The Tweedie distribution with power `1 < p < 2`: compound Poisson with
2//! gamma severities.
3
4use prospicio_core::{Error, Result};
5use prospicio_math::roots::bisect;
6use prospicio_math::special::ln_gamma;
7
8use crate::distribution::{Distribution, check_probability};
9use crate::gamma::Gamma;
10use crate::severity::Severity;
11
12/// Tweedie distribution with mean `μ`, dispersion `φ` and power
13/// `1 < p < 2`: variance `φ μ^p`, a point mass at 0 and a continuous
14/// density above it.
15///
16/// It is the compound Poisson sum `Y = X_1 + … + X_N` with
17///
18/// ```text
19/// N ~ Poisson(λ),        λ = μ^(2-p) / (φ (2 - p))
20/// X ~ Gamma(α, θ),       α = (2 - p) / (p - 1),   θ = φ (p - 1) μ^(p-1)
21/// ```
22///
23/// the GLM family for pure premium (losses per exposure), where claim
24/// counts and severities are not modelled separately. Every quantity is a
25/// Poisson-weighted sum over the number of claims `n`, of the matching
26/// quantity of `Gamma(nα, θ)`, summed until the remaining terms are
27/// negligible:
28///
29/// - `P(Y = 0) = e^(-λ)`;
30/// - the distribution function, survival function and layer moments from
31///   the gamma's, each tail summed directly so both keep their precision;
32/// - the density (Dunn & Smyth's series), in log space.
33///
34/// # Example
35///
36/// ```
37/// use prospicio_prob::{Distribution, Severity, Tweedie};
38///
39/// let y = Tweedie::new(500.0, 40.0, 1.6).unwrap();
40/// assert!((y.mean() - 500.0).abs() < 1e-9);
41/// assert!((y.variance() - 40.0 * 500f64.powf(1.6)).abs() < 1e-6);
42/// // P(Y = 0) = e^(-λ).
43/// assert!((y.cdf(0.0) - (-y.lambda()).exp()).abs() < 1e-15);
44/// assert!((y.lev(800.0) + y.stop_loss(800.0) - 500.0).abs() < 1e-9);
45/// ```
46#[derive(Debug, Clone, Copy, PartialEq)]
47pub struct Tweedie {
48    mean: f64,
49    dispersion: f64,
50    power: f64,
51    lambda: f64,
52    severity: Gamma,
53}
54
55impl Tweedie {
56    /// Tweedie with mean `μ > 0`, dispersion `φ > 0` and power `p` in
57    /// `(1, 2)`.
58    pub fn new(mean: f64, dispersion: f64, power: f64) -> Result<Self> {
59        positive("mean", mean)?;
60        positive("dispersion", dispersion)?;
61        if !(power > 1.0 && power < 2.0) {
62            return Err(Error::InvalidParameter {
63                name: "power",
64                value: power,
65                reason: "must be in (1, 2)",
66            });
67        }
68        let lambda = mean.powf(2.0 - power) / (dispersion * (2.0 - power));
69        let shape = (2.0 - power) / (power - 1.0);
70        let scale = dispersion * (power - 1.0) * mean.powf(power - 1.0);
71        Ok(Self {
72            mean,
73            dispersion,
74            power,
75            lambda,
76            severity: Gamma::new(shape, scale)?,
77        })
78    }
79
80    /// The Tweedie equal to a Poisson(`lambda`) number of
81    /// Gamma(`shape`, `scale`) losses: power `(α + 2) / (α + 1)`, mean
82    /// `λαθ`.
83    ///
84    /// ```
85    /// use prospicio_prob::Tweedie;
86    ///
87    /// let y = Tweedie::from_poisson_gamma(3.0, 2.0, 100.0).unwrap();
88    /// assert!((y.power() - 4.0 / 3.0).abs() < 1e-15);
89    /// assert!((y.lambda() - 3.0).abs() < 1e-12);
90    /// ```
91    pub fn from_poisson_gamma(lambda: f64, shape: f64, scale: f64) -> Result<Self> {
92        positive("lambda", lambda)?;
93        Gamma::new(shape, scale)?;
94        let power = (shape + 2.0) / (shape + 1.0);
95        let mean = lambda * shape * scale;
96        let dispersion = scale / ((power - 1.0) * mean.powf(power - 1.0));
97        Self::new(mean, dispersion, power)
98    }
99
100    /// Mean `μ`.
101    pub fn mean_param(&self) -> f64 {
102        self.mean
103    }
104
105    /// Dispersion `φ`.
106    pub fn dispersion(&self) -> f64 {
107        self.dispersion
108    }
109
110    /// Power `p`.
111    pub fn power(&self) -> f64 {
112        self.power
113    }
114
115    /// Poisson mean `λ` of the number of losses.
116    pub fn lambda(&self) -> f64 {
117        self.lambda
118    }
119
120    /// The gamma distribution of each loss.
121    pub fn severity(&self) -> Gamma {
122        self.severity
123    }
124
125    /// Log density at `y > 0` (the continuous part); at `y = 0`, the log
126    /// of the point mass `-λ`. `-inf` below 0.
127    ///
128    /// ```
129    /// use prospicio_prob::Tweedie;
130    ///
131    /// // λ = 1 with unit exponential losses:
132    /// // f(y) = e^(-1-y) Σ_n y^(n-1) / (n! (n-1)!).
133    /// let y = Tweedie::from_poisson_gamma(1.0, 1.0, 1.0).unwrap();
134    /// let (mut series, mut term) = (0.0, 1.0); // term = 2^(n-1) / (n! (n-1)!)
135    /// for n in 1..40 {
136    ///     series += term;
137    ///     term *= 2.0 / (f64::from(n + 1) * f64::from(n));
138    /// }
139    /// assert!((y.ln_pdf(2.0) - (-3.0 + f64::ln(series))).abs() < 1e-13);
140    /// ```
141    pub fn ln_pdf(&self, y: f64) -> f64 {
142        if y < 0.0 || y.is_nan() {
143            return f64::NEG_INFINITY;
144        }
145        if y == 0.0 {
146            return -self.lambda;
147        }
148        // ln of term n: ln Pois(n; λ) + ln Gamma(nα, θ).pdf(y). Terms rise to
149        // a maximum and fall; sum them relative to the largest.
150        let (alpha, theta) = (self.severity.shape(), self.severity.scale());
151        let ln_term = |n: f64| {
152            n * self.lambda.ln() - self.lambda - ln_gamma(n + 1.0) + (n * alpha - 1.0) * y.ln()
153                - y / theta
154                - ln_gamma(n * alpha)
155                - n * alpha * theta.ln()
156        };
157        // The terms peak near n where the Poisson and gamma pulls balance.
158        let peak = {
159            let mut best = (1.0, ln_term(1.0));
160            let guess = (y / (alpha * theta)).max(self.lambda).max(1.0);
161            for n in [guess.floor(), guess.ceil(), self.lambda.floor().max(1.0)] {
162                let v = ln_term(n.max(1.0));
163                if v > best.1 {
164                    best = (n.max(1.0), v);
165                }
166            }
167            // Walk uphill to the exact peak.
168            let (mut n, mut v) = best;
169            loop {
170                let up = ln_term(n + 1.0);
171                if up > v {
172                    (n, v) = (n + 1.0, up);
173                    continue;
174                }
175                if n > 1.0 {
176                    let down = ln_term(n - 1.0);
177                    if down > v {
178                        (n, v) = (n - 1.0, down);
179                        continue;
180                    }
181                }
182                break (n, v);
183            }
184        };
185        let (n0, v0) = peak;
186        let mut sum = 1.0;
187        let mut n = n0 + 1.0;
188        loop {
189            let r = (ln_term(n) - v0).exp();
190            sum += r;
191            if r <= 1e-17 * sum {
192                break;
193            }
194            n += 1.0;
195        }
196        let mut n = n0 - 1.0;
197        while n >= 1.0 {
198            let r = (ln_term(n) - v0).exp();
199            sum += r;
200            if r <= 1e-17 * sum {
201                break;
202            }
203            n -= 1.0;
204        }
205        v0 + sum.ln()
206    }
207
208    /// `Σ_{n≥1} Pois(n; λ) g(Gamma(nα, θ))` for a `g` that grows with `n`
209    /// at most polynomially. Sums up from `n = 1` past both the Poisson
210    /// mode and `beyond` (where `g` levels off), then until a term is
211    /// negligible.
212    fn poisson_sum(&self, beyond: f64, mut g: impl FnMut(&Gamma) -> f64) -> f64 {
213        let (alpha, theta) = (self.severity.shape(), self.severity.scale());
214        let ln_lambda = self.lambda.ln();
215        let floor = self.lambda.max(beyond);
216        let mut sum = 0.0;
217        for n in 1..10_000_000u32 {
218            let nf = f64::from(n);
219            let w = (nf * ln_lambda - self.lambda - ln_gamma(nf + 1.0)).exp();
220            let term = if w == 0.0 {
221                0.0
222            } else {
223                w * g(&Gamma::new(nf * alpha, theta).expect("valid shape and scale"))
224            };
225            sum += term;
226            if nf > floor && (term <= 1e-17 * sum || w == 0.0) {
227                break;
228            }
229        }
230        sum
231    }
232
233    /// The `n` beyond which `Gamma(nα, θ)` sits mostly above `y`.
234    fn claims_to_reach(&self, y: f64) -> f64 {
235        let k = y / (self.severity.shape() * self.severity.scale());
236        k + 10.0 * k.sqrt() + 10.0
237    }
238}
239
240impl Distribution for Tweedie {
241    fn mean(&self) -> f64 {
242        self.mean
243    }
244
245    fn variance(&self) -> f64 {
246        self.dispersion * self.mean.powf(self.power)
247    }
248
249    fn cdf(&self, y: f64) -> f64 {
250        if y < 0.0 {
251            return 0.0;
252        }
253        let zero = (-self.lambda).exp();
254        if y == 0.0 {
255            return zero;
256        }
257        (zero + self.poisson_sum(0.0, |g| g.cdf(y))).min(1.0)
258    }
259
260    fn survival(&self, y: f64) -> f64 {
261        if y < 0.0 {
262            return 1.0;
263        }
264        if y == 0.0 {
265            return -(-self.lambda).exp_m1();
266        }
267        self.poisson_sum(self.claims_to_reach(y), |g| g.survival(y))
268            .min(1.0)
269    }
270
271    /// 0 when `p` is within the point mass; otherwise by bisection on the
272    /// distribution or survival function, whichever is the smaller tail.
273    fn quantile(&self, p: f64) -> Result<f64> {
274        check_probability(p)?;
275        if p <= (-self.lambda).exp() {
276            return Ok(0.0);
277        }
278        if p == 1.0 {
279            return Ok(f64::INFINITY);
280        }
281        let below = |y: f64| {
282            if p <= 0.5 {
283                self.cdf(y) < p
284            } else {
285                self.survival(y) > 1.0 - p
286            }
287        };
288        let mut hi = self.mean + self.variance().sqrt();
289        while below(hi) {
290            hi *= 2.0;
291        }
292        Ok(bisect(0.0, hi, below))
293    }
294}
295
296impl Severity for Tweedie {
297    fn lev(&self, limit: f64) -> f64 {
298        if limit <= 0.0 {
299            return limit;
300        }
301        if limit == f64::INFINITY {
302            return self.mean;
303        }
304        self.poisson_sum(0.0, |g| g.lev(limit))
305    }
306
307    fn stop_loss(&self, retention: f64) -> f64 {
308        if retention <= 0.0 {
309            return self.mean - retention;
310        }
311        if retention == f64::INFINITY {
312            return 0.0;
313        }
314        self.poisson_sum(self.claims_to_reach(retention), |g| g.stop_loss(retention))
315    }
316
317    fn layer(&self, limit: f64, attachment: f64) -> f64 {
318        let a = attachment.max(0.0);
319        self.poisson_sum(self.claims_to_reach(a + limit.min(1e300)), |g| {
320            g.layer(limit, a)
321        })
322    }
323
324    fn layer_second_moment(&self, limit: f64, attachment: f64) -> f64 {
325        let a = attachment.max(0.0);
326        self.poisson_sum(self.claims_to_reach(a + limit.min(1e300)), |g| {
327            g.layer_second_moment(limit, a)
328        })
329    }
330}
331
332fn positive(name: &'static str, value: f64) -> Result<()> {
333    if value.is_finite() && value > 0.0 {
334        Ok(())
335    } else {
336        Err(Error::InvalidParameter {
337            name,
338            value,
339            reason: "must be finite and positive",
340        })
341    }
342}
343
344#[cfg(test)]
345mod tests {
346    use super::*;
347
348    #[test]
349    fn parameterizations_agree() {
350        let y = Tweedie::new(500.0, 40.0, 1.6).unwrap();
351        let z = Tweedie::from_poisson_gamma(y.lambda(), y.severity().shape(), y.severity().scale())
352            .unwrap();
353        assert!((z.mean() / 500.0 - 1.0).abs() < 1e-14);
354        assert!((z.dispersion() / 40.0 - 1.0).abs() < 1e-12);
355        assert!((z.power() - 1.6).abs() < 1e-14);
356        assert!(Tweedie::new(1.0, 1.0, 2.0).is_err());
357        assert!(Tweedie::new(1.0, 1.0, 1.0).is_err());
358    }
359
360    #[test]
361    fn moments_and_layers_from_the_series() {
362        let y = Tweedie::new(500.0, 40.0, 1.6).unwrap();
363        // The unlimited layer is the whole distribution.
364        assert!((y.layer(f64::INFINITY, 0.0) / y.mean() - 1.0).abs() < 1e-12);
365        let m2 = y.layer_second_moment(f64::INFINITY, 0.0);
366        let want = y.variance() + y.mean() * y.mean();
367        assert!((m2 / want - 1.0).abs() < 1e-12, "{m2} {want}");
368        for d in [10.0, 500.0, 5_000.0] {
369            assert!((y.lev(d) + y.stop_loss(d) - 500.0).abs() < 1e-9);
370            assert!((y.cdf(d) + y.survival(d) - 1.0).abs() < 1e-14);
371        }
372    }
373
374    #[test]
375    fn density_integrates_to_the_distribution_function() {
376        let y = Tweedie::new(10.0, 2.0, 1.4).unwrap();
377        // ∫_a^b f by composite Gauss–Legendre equals F(b) - F(a).
378        let (a, b) = (0.5, 30.0);
379        let panels = 400;
380        let h = (b - a) / f64::from(panels);
381        let integral: f64 = (0..panels)
382            .map(|i| {
383                let lo = a + h * f64::from(i);
384                prospicio_math::integrate::gauss_legendre(
385                    |x| Ok::<_, ()>(y.ln_pdf(x).exp()),
386                    lo,
387                    lo + h,
388                )
389                .unwrap()
390            })
391            .sum();
392        assert!(
393            (integral - (y.cdf(b) - y.cdf(a))).abs() < 1e-12,
394            "{integral}"
395        );
396    }
397
398    #[test]
399    fn quantiles_invert_and_respect_the_atom() {
400        let y = Tweedie::new(10.0, 2.0, 1.4).unwrap();
401        let p0 = y.cdf(0.0);
402        assert_eq!(y.quantile(p0 * 0.5), Ok(0.0));
403        for p in [0.5, 0.9, 0.999999] {
404            let x = y.quantile(p).unwrap();
405            assert!((y.cdf(x) - p).abs() < 1e-12, "{p}");
406        }
407    }
408}
prospicio_prob/tweedie.rs

prospicio_prob/
tweedie.rs
