prospicio_prob/
loglogistic.rs

1//! The loglogistic (Fisk) distribution.
2
3use std::f64::consts::PI;
4
5use prospicio_core::{Error, Result};
6use prospicio_math::special::{beta_inc, ln_gamma};
7
8use crate::distribution::{Distribution, check_probability};
9use crate::severity::Severity;
10
11/// Loglogistic distribution with shape `α` and scale `θ`:
12/// `F(x) = (x/θ)^α / (1 + (x/θ)^α)`, as SciPy's `fisk(c=α, scale=θ)` and
13/// actuar's `dllogis(shape, scale)`. Its median is `θ`.
14///
15/// It is also the growth curve of Clark's LDF method, `G(t) = t^ω / (t^ω + θ^ω)`,
16/// so a fitted curve is this distribution's [`cdf`](Distribution::cdf).
17///
18/// The tail is Pareto-like: `E[X^j]` is finite only for `j < α`, so the
19/// mean is infinite for `α <= 1` and the variance for `α <= 2`. Limited and
20/// layer moments are finite for every `α`. With `v = F(u)`,
21/// `E[min(X, u)^j] = θ^j B(1 + j/α, 1 - j/α; v) + u^j (1 - v)`, using the
22/// unnormalized incomplete beta, which for `j >= α` has a non-positive
23/// second argument and is reached by recurrence.
24///
25/// ```
26/// use prospicio_prob::{Distribution, Loglogistic, Severity};
27///
28/// // Shape 1: F(x) = x / (x + θ) and E[min(X, u)] = θ ln(1 + u/θ).
29/// let d = Loglogistic::new(1.0, 2.0).unwrap();
30/// assert!((d.cdf(3.0) - 0.6).abs() < 1e-15);
31/// assert!((d.lev(3.0) - 2.0 * 2.5f64.ln()).abs() < 1e-13);
32/// assert!(d.mean().is_infinite());
33/// ```
34#[derive(Debug, Clone, Copy, PartialEq)]
35pub struct Loglogistic {
36    shape: f64,
37    scale: f64,
38}
39
40impl Loglogistic {
41    /// Loglogistic with shape `α > 0` and scale `θ > 0`.
42    pub fn new(shape: f64, scale: f64) -> Result<Self> {
43        for (name, v) in [("shape", shape), ("scale", scale)] {
44            if !(v.is_finite() && v > 0.0) {
45                return Err(Error::InvalidParameter {
46                    name,
47                    value: v,
48                    reason: "must be finite and positive",
49                });
50            }
51        }
52        Ok(Self { shape, scale })
53    }
54
55    /// Shape `α`.
56    pub fn shape(&self) -> f64 {
57        self.shape
58    }
59
60    /// Scale `θ`, the median.
61    pub fn scale(&self) -> f64 {
62        self.scale
63    }
64
65    /// `E[X^j] = θ^j (πj/α) / sin(πj/α)`, infinite for `j >= α`.
66    fn raw_moment(&self, j: f64) -> f64 {
67        if j >= self.shape {
68            return f64::INFINITY;
69        }
70        let r = PI * j / self.shape;
71        self.scale.powf(j) * r / r.sin()
72    }
73
74    /// `(F(u), S(u))`, each computed without cancellation.
75    fn probs(&self, u: f64) -> (f64, f64) {
76        let z = (u / self.scale).powf(self.shape);
77        (z / (1.0 + z), 1.0 / (1.0 + z))
78    }
79
80    /// `E[min(X, u)^j]` for `u > 0`.
81    fn limited_moment(&self, j: i32, u: f64) -> f64 {
82        if u == f64::INFINITY {
83            return self.raw_moment(f64::from(j));
84        }
85        let (v, s) = self.probs(u);
86        let r = f64::from(j) / self.shape;
87        self.scale.powi(j) * beta_lower(1.0 + r, 1.0 - r, v) + u.powi(j) * s
88    }
89
90    /// `E[X^j; X > u] - u^j S(u)`, the tail part of
91    /// `E[X^j] - E[min(X, u)^j]`, for `j < α` and `u > 0`; it uses
92    /// `E[X^j; X > u] = θ^j B(1 - j/α, 1 + j/α; S(u))`.
93    fn tail_moment(&self, j: i32, u: f64) -> f64 {
94        if u == f64::INFINITY {
95            return 0.0;
96        }
97        let (_, s) = self.probs(u);
98        let r = f64::from(j) / self.shape;
99        (self.scale.powi(j) * beta_lower(1.0 - r, 1.0 + r, s) - u.powi(j) * s).max(0.0)
100    }
101
102    /// `∫_a^b x^(j-1) S(x) dx` for `0 <= a < b <= ∞`: limited-moment
103    /// differences below `c = θ 2^(1/α)` (where `S <= 1/3` begins), and the
104    /// tail series above it, so a far layer does not cancel.
105    fn partial(&self, j: i32, a: f64, b: f64) -> f64 {
106        let below = |u: f64| {
107            if u <= 0.0 {
108                0.0
109            } else {
110                self.limited_moment(j, u) / f64::from(j)
111            }
112        };
113        let c = self.scale * 2f64.powf(1.0 / self.shape);
114        if b <= c {
115            return below(b) - below(a);
116        }
117        if a >= c {
118            return self.tail_series(j, a, b);
119        }
120        below(c) - below(a) + self.tail_series(j, c, b)
121    }
122
123    /// `∫_a^b x^(j-1) S(x) dx` for `θ 2^(1/α) <= a < b`. With
124    /// `r = (θ/x)^α <= 1/2`, it is `(θ^j/α) ∫ r^(-j/α) / (1 + r) dr` over
125    /// `[r_b, r_a]`, integrated term by term in the geometric series of
126    /// `1 / (1 + r)`; infinite when `b = ∞` and `j >= α`.
127    fn tail_series(&self, j: i32, a: f64, b: f64) -> f64 {
128        let ra = (self.scale / a).powf(self.shape);
129        let rb = if b == f64::INFINITY {
130            0.0
131        } else {
132            (self.scale / b).powf(self.shape)
133        };
134        let ln_ratio = (rb / ra).ln();
135        let r = f64::from(j) / self.shape;
136        let mut sum = 0.0;
137        for k in 0..200 {
138            // ∫_{r_b}^{r_a} r^(e - 1) dr with e = k + 1 - j/α.
139            let e = f64::from(k) + 1.0 - r;
140            let term = if e.abs() < 1e-12 {
141                -ln_ratio
142            } else if rb == 0.0 {
143                if e > 0.0 {
144                    ra.powf(e) / e
145                } else {
146                    f64::INFINITY
147                }
148            } else {
149                -ra.powf(e) * (e * ln_ratio).exp_m1() / e
150            };
151            if !term.is_finite() {
152                return f64::INFINITY;
153            }
154            sum += if k % 2 == 0 { term } else { -term };
155            if k > 0 && term.abs() < 1e-17 * sum.abs() {
156                break;
157            }
158        }
159        self.scale.powi(j) / self.shape * sum
160    }
161}
162
163/// The unnormalized lower incomplete beta `B(a, b; x) = ∫_0^x t^(a-1) (1-t)^(b-1) dt`
164/// for `a > 0`, any `b`, and `0 <= x < 1`.
165///
166/// For `b > 0` it is `I_x(a, b) B(a, b)`. Otherwise it steps `b` up to a
167/// positive value and back down by
168/// `B(a, b; x) = ((a + b) B(a, b + 1; x) - x^a (1 - x)^b) / b`. A zero `b`
169/// on the way only arises here with integer `a` (`a + b = 2`), where
170/// `B(a, 0; x) = -ln(1 - x) - Σ_{k<a} x^k / k`.
171fn beta_lower(a: f64, b: f64, x: f64) -> f64 {
172    if x <= 0.0 {
173        return 0.0;
174    }
175    if b > 0.0 {
176        let ln_b = ln_gamma(a) + ln_gamma(b) - ln_gamma(a + b);
177        return beta_inc(a, b, x) * ln_b.exp();
178    }
179    // The smallest c = b + n that is positive, or zero.
180    let steps = (-b).floor();
181    let c = b + steps;
182    let (mut value, mut c) = if c == 0.0 {
183        (beta_zero(a, x), 0.0)
184    } else {
185        let c = c + 1.0;
186        let ln_b = ln_gamma(a) + ln_gamma(c) - ln_gamma(a + c);
187        (beta_inc(a, c, x) * ln_b.exp(), c)
188    };
189    while c > b {
190        c -= 1.0;
191        value = ((a + c) * value - x.powf(a) * (c * (-x).ln_1p()).exp()) / c;
192    }
193    value
194}
195
196/// `B(a, 0; x)` for integer `a >= 1`: a power series for small `x`, where
197/// the closed form would cancel.
198fn beta_zero(a: f64, x: f64) -> f64 {
199    let m = a.round() as i32;
200    if x < 0.5 {
201        let mut sum = 0.0;
202        let mut term = x.powi(m);
203        for n in 0..2000 {
204            let add = term / f64::from(m + n);
205            sum += add;
206            if add < 1e-17 * sum {
207                break;
208            }
209            term *= x;
210        }
211        sum
212    } else {
213        let partial: f64 = (1..m).map(|k| x.powi(k) / f64::from(k)).sum();
214        -(-x).ln_1p() - partial
215    }
216}
217
218impl Distribution for Loglogistic {
219    fn mean(&self) -> f64 {
220        self.raw_moment(1.0)
221    }
222
223    fn variance(&self) -> f64 {
224        if self.shape <= 2.0 {
225            return f64::INFINITY;
226        }
227        let m = self.mean();
228        self.raw_moment(2.0) - m * m
229    }
230
231    fn cdf(&self, x: f64) -> f64 {
232        if x <= 0.0 {
233            return 0.0;
234        }
235        self.probs(x).0
236    }
237
238    fn survival(&self, x: f64) -> f64 {
239        if x <= 0.0 {
240            return 1.0;
241        }
242        self.probs(x).1
243    }
244
245    /// `θ (p / (1 - p))^(1/α)`.
246    fn quantile(&self, p: f64) -> Result<f64> {
247        check_probability(p)?;
248        Ok(self.scale * (p / (1.0 - p)).powf(1.0 / self.shape))
249    }
250}
251
252impl Severity for Loglogistic {
253    fn lev(&self, limit: f64) -> f64 {
254        if limit <= 0.0 {
255            return limit;
256        }
257        self.limited_moment(1, limit)
258    }
259
260    /// From the tail, so it does not cancel against the mean; infinite for
261    /// `α <= 1`.
262    fn stop_loss(&self, retention: f64) -> f64 {
263        if retention <= 0.0 {
264            return self.mean() - retention;
265        }
266        if self.shape <= 1.0 {
267            return f64::INFINITY;
268        }
269        self.tail_moment(1, retention)
270    }
271
272    /// `∫_a^(a+l) S(x) dx`, by the tail series in the tail.
273    fn layer(&self, limit: f64, attachment: f64) -> f64 {
274        let a = attachment.max(0.0);
275        self.partial(1, a, a + limit)
276    }
277
278    /// `2 ∫_a^b (x - a) S(x) dx` with `b = a + limit`.
279    fn layer_second_moment(&self, limit: f64, attachment: f64) -> f64 {
280        let a = attachment.max(0.0);
281        let b = a + limit;
282        let p1 = self.partial(1, a, b);
283        if p1 == 0.0 {
284            return 0.0;
285        }
286        (2.0 * (self.partial(2, a, b) - a * p1)).max(0.0)
287    }
288}
289
290#[cfg(test)]
291mod tests {
292    use super::*;
293
294    #[test]
295    fn identities() {
296        for shape in [0.5, 0.8, 1.0, 1.5, 2.0, 4.0] {
297            let d = Loglogistic::new(shape, 100.0).unwrap();
298            for u in [1.0, 100.0, 20_000.0] {
299                // Far in the tail, 1 - F(u) has too few digits to invert.
300                if d.cdf(u) < 0.99 {
301                    let x = d.quantile(d.cdf(u)).unwrap();
302                    assert!((x / u - 1.0).abs() < 1e-12);
303                }
304                if shape > 1.0 {
305                    let total = d.lev(u) + d.stop_loss(u);
306                    assert!((total / d.mean() - 1.0).abs() < 1e-12, "{shape} {u}");
307                }
308            }
309            let stacked = d.layer(500.0, 0.0) + d.layer(500.0, 500.0);
310            assert!((stacked / d.layer(1000.0, 0.0) - 1.0).abs() < 1e-12);
311            assert!(d.layer_second_moment(200.0, 300.0) > 0.0);
312        }
313        let d = Loglogistic::new(4.0, 3.0).unwrap();
314        let m2 = d.layer_second_moment(f64::INFINITY, 0.0);
315        assert!((m2 / (d.variance() + d.mean().powi(2)) - 1.0).abs() < 1e-12);
316        assert!(Loglogistic::new(0.0, 1.0).is_err());
317    }
318
319    /// Layer moments by the tail series and by limited moments agree.
320    #[test]
321    fn layer_second_moment_branches() {
322        let d = Loglogistic::new(3.0, 50.0).unwrap();
323        let (a, b) = (40.0, 90.0);
324        let direct = d.lev(b) - d.lev(a);
325        assert!((d.layer(b - a, a) / direct - 1.0).abs() < 1e-12);
326        let m = |j: i32, u: f64| d.limited_moment(j, u);
327        let direct = m(2, b) - m(2, a) - 2.0 * a * (m(1, b) - m(1, a));
328        let tails = d.layer_second_moment(b - a, a);
329        assert!((tails / direct - 1.0).abs() < 1e-11);
330    }
331
332    #[test]
333    fn beta_lower_matches_integration() {
334        // B(a, b; x) by the midpoint rule on a fine grid, for b of both signs.
335        for (a, b, x) in [
336            (1.5, 0.5, 0.7),
337            (2.0, -1.0, 0.9),
338            (2.25, -0.25, 0.4),
339            (3.0, -1.0, 0.2),
340        ] {
341            let n = 200_000;
342            let h = x / f64::from(n);
343            let num: f64 = (0..n)
344                .map(|i| {
345                    let t = (f64::from(i) + 0.5) * h;
346                    t.powf(a - 1.0) * (1.0 - t).powf(b - 1.0) * h
347                })
348                .sum();
349            let got = beta_lower(a, b, x);
350            assert!(
351                (got / num - 1.0).abs() < 1e-8,
352                "{a} {b} {x}: {got} vs {num}"
353            );
354        }
355    }
356}
prospicio_prob/loglogistic.rs

prospicio_prob/
loglogistic.rs
