prospicio_prob/
weibull.rs

1//! The Weibull distribution.
2
3use prospicio_core::{Error, Result};
4use prospicio_math::special::{gamma_inc, ln_gamma};
5
6use crate::distribution::{Distribution, check_probability};
7use crate::severity::Severity;
8
9/// Weibull distribution with shape `k` and scale `λ`:
10/// `P(X > x) = exp(-(x/λ)^k)`, as SciPy's `weibull_min(c=k, scale=λ)` and
11/// R's `dweibull(shape, scale)`.
12///
13/// Below shape 1 the tail is heavier than exponential (but every moment
14/// exists); above 1 it is lighter. Layer moments use the regularized
15/// incomplete gamma: with `z = (u/λ)^k`,
16/// `E[X^j; X > u] = λ^j Γ(1 + j/k) Q(1 + j/k, z)`.
17///
18/// ```
19/// use prospicio_prob::{Distribution, Severity, Weibull};
20///
21/// // Shape 1 is the exponential.
22/// let e = Weibull::new(1.0, 2.0).unwrap();
23/// assert!((e.mean() - 2.0).abs() < 1e-14);
24/// assert!((e.stop_loss(3.0) - 2.0 * (-1.5f64).exp()).abs() < 1e-14);
25/// ```
26#[derive(Debug, Clone, Copy, PartialEq)]
27pub struct Weibull {
28    shape: f64,
29    scale: f64,
30}
31
32impl Weibull {
33    /// Weibull with shape `k > 0` and scale `λ > 0`.
34    pub fn new(shape: f64, scale: f64) -> Result<Self> {
35        for (name, v) in [("shape", shape), ("scale", scale)] {
36            if !(v.is_finite() && v > 0.0) {
37                return Err(Error::InvalidParameter {
38                    name,
39                    value: v,
40                    reason: "must be finite and positive",
41                });
42            }
43        }
44        Ok(Self { shape, scale })
45    }
46
47    /// Shape `k`.
48    pub fn shape(&self) -> f64 {
49        self.shape
50    }
51
52    /// Scale `λ`.
53    pub fn scale(&self) -> f64 {
54        self.scale
55    }
56
57    /// `E[X^j]`.
58    fn raw_moment(&self, j: f64) -> f64 {
59        self.scale.powf(j) * ln_gamma(1.0 + j / self.shape).exp()
60    }
61
62    /// `E[X^j; X > u] - u^j S(u)` for `j` in 1 and 2: the tail part of
63    /// `E[X^j] - E[min(X, u)^j]`.
64    fn tail_moment(&self, j: i32, u: f64) -> f64 {
65        if u <= 0.0 {
66            return self.raw_moment(f64::from(j)) - u.powi(j);
67        }
68        if u == f64::INFINITY {
69            return 0.0;
70        }
71        let z = (u / self.scale).powf(self.shape);
72        let jf = f64::from(j);
73        let q = gamma_inc(1.0 + jf / self.shape, z).1;
74        (self.raw_moment(jf) * q - u.powi(j) * (-z).exp()).max(0.0)
75    }
76}
77
78impl Distribution for Weibull {
79    fn mean(&self) -> f64 {
80        self.raw_moment(1.0)
81    }
82
83    fn variance(&self) -> f64 {
84        let m = self.mean();
85        self.raw_moment(2.0) - m * m
86    }
87
88    fn cdf(&self, x: f64) -> f64 {
89        if x <= 0.0 {
90            return 0.0;
91        }
92        -(-(x / self.scale).powf(self.shape)).exp_m1()
93    }
94
95    fn survival(&self, x: f64) -> f64 {
96        if x <= 0.0 {
97            return 1.0;
98        }
99        (-(x / self.scale).powf(self.shape)).exp()
100    }
101
102    /// `λ (-ln(1 - p))^(1/k)`, with `ln(1 - p)` taken as `ln_1p(-p)`.
103    fn quantile(&self, p: f64) -> Result<f64> {
104        check_probability(p)?;
105        Ok(self.scale * (-(-p).ln_1p()).powf(1.0 / self.shape))
106    }
107}
108
109impl Severity for Weibull {
110    /// `E[min(X, u)] = λ Γ(1 + 1/k) P(1 + 1/k, z) + u e^(-z)`.
111    fn lev(&self, limit: f64) -> f64 {
112        if limit <= 0.0 {
113            return limit;
114        }
115        if limit == f64::INFINITY {
116            return self.mean();
117        }
118        let z = (limit / self.scale).powf(self.shape);
119        self.mean() * gamma_inc(1.0 + 1.0 / self.shape, z).0 + limit * (-z).exp()
120    }
121
122    /// From the tail, so it does not cancel against the mean.
123    fn stop_loss(&self, retention: f64) -> f64 {
124        if retention <= 0.0 {
125            return self.mean() - retention;
126        }
127        self.tail_moment(1, retention)
128    }
129
130    /// LEV differences below the mean, stop-loss differences above it.
131    fn layer(&self, limit: f64, attachment: f64) -> f64 {
132        let a = attachment.max(0.0);
133        if a <= self.mean() {
134            self.lev(a + limit) - self.lev(a)
135        } else {
136            self.stop_loss(a) - self.stop_loss(a + limit)
137        }
138    }
139
140    /// `t_2(a) - t_2(b) - 2a (t_1(a) - t_1(b))` with `b = a + limit`.
141    fn layer_second_moment(&self, limit: f64, attachment: f64) -> f64 {
142        let a = attachment.max(0.0);
143        let b = a + limit;
144        let t = |j: i32, u: f64| self.tail_moment(j, u);
145        (t(2, a) - t(2, b) - 2.0 * a * (t(1, a) - t(1, b))).max(0.0)
146    }
147}
148
149#[cfg(test)]
150mod tests {
151    use super::*;
152
153    #[test]
154    fn identities() {
155        let d = Weibull::new(0.7, 1000.0).unwrap();
156        for u in [10.0, 1000.0, 20_000.0] {
157            assert!((d.lev(u) + d.stop_loss(u) - d.mean()).abs() < 1e-10 * d.mean());
158            let x = d.quantile(d.cdf(u)).unwrap();
159            assert!((x / u - 1.0).abs() < 1e-12);
160        }
161        let m2 = d.layer_second_moment(f64::INFINITY, 0.0);
162        assert!((m2 / (d.variance() + d.mean().powi(2)) - 1.0).abs() < 1e-12);
163        let stacked = d.layer(500.0, 0.0) + d.layer(500.0, 500.0);
164        assert!((stacked - d.layer(1000.0, 0.0)).abs() < 1e-9);
165        assert!(Weibull::new(0.0, 1.0).is_err());
166    }
167}
prospicio_prob/weibull.rs

prospicio_prob/
weibull.rs
