prospicio_prob/
lognormal.rs

1//! The lognormal distribution.
2
3use prospicio_core::{Error, Result};
4use prospicio_math::special::{norm_cdf, norm_quantile};
5
6use crate::distribution::{Distribution, check_probability};
7use crate::severity::Severity;
8
9/// Lognormal distribution: `ln X ~ Normal(meanlog, sdlog^2)`.
10///
11/// Parameterized as in SciPy (`s = sdlog`, `scale = exp(meanlog)`), R
12/// (`meanlog`, `sdlog`) and actuar.
13///
14/// # Example
15///
16/// ```
17/// use prospicio_prob::{Distribution, Lognormal};
18///
19/// let d = Lognormal::from_mean_cv(1000.0, 0.5).unwrap();
20/// assert!((d.mean() - 1000.0).abs() < 1e-9);
21/// assert!((d.std_dev() - 500.0).abs() < 1e-9);
22/// ```
23#[derive(Debug, Clone, Copy, PartialEq)]
24pub struct Lognormal {
25    meanlog: f64,
26    sdlog: f64,
27}
28
29impl Lognormal {
30    /// Lognormal with the given log-scale mean and standard deviation.
31    pub fn new(meanlog: f64, sdlog: f64) -> Result<Self> {
32        if !meanlog.is_finite() {
33            return Err(Error::InvalidParameter {
34                name: "meanlog",
35                value: meanlog,
36                reason: "must be finite",
37            });
38        }
39        if !sdlog.is_finite() || sdlog <= 0.0 {
40            return Err(Error::InvalidParameter {
41                name: "sdlog",
42                value: sdlog,
43                reason: "must be finite and positive",
44            });
45        }
46        Ok(Self { meanlog, sdlog })
47    }
48
49    /// Lognormal with the given mean and coefficient of variation, the way
50    /// severity assumptions are usually stated.
51    pub fn from_mean_cv(mean: f64, cv: f64) -> Result<Self> {
52        if !mean.is_finite() || mean <= 0.0 {
53            return Err(Error::InvalidParameter {
54                name: "mean",
55                value: mean,
56                reason: "must be finite and positive",
57            });
58        }
59        if !cv.is_finite() || cv <= 0.0 {
60            return Err(Error::InvalidParameter {
61                name: "cv",
62                value: cv,
63                reason: "must be finite and positive",
64            });
65        }
66        let sigma2 = cv.mul_add(cv, 1.0).ln();
67        Self::new(mean.ln() - 0.5 * sigma2, sigma2.sqrt())
68    }
69
70    /// Log-scale mean.
71    pub fn meanlog(&self) -> f64 {
72        self.meanlog
73    }
74
75    /// Log-scale standard deviation.
76    pub fn sdlog(&self) -> f64 {
77        self.sdlog
78    }
79}
80
81impl Distribution for Lognormal {
82    fn mean(&self) -> f64 {
83        (self.meanlog + 0.5 * self.sdlog * self.sdlog).exp()
84    }
85
86    fn variance(&self) -> f64 {
87        let s2 = self.sdlog * self.sdlog;
88        s2.exp_m1() * (2.0 * self.meanlog + s2).exp()
89    }
90
91    fn cdf(&self, x: f64) -> f64 {
92        if x <= 0.0 {
93            return 0.0;
94        }
95        norm_cdf((x.ln() - self.meanlog) / self.sdlog)
96    }
97
98    fn survival(&self, x: f64) -> f64 {
99        if x <= 0.0 {
100            return 1.0;
101        }
102        norm_cdf((self.meanlog - x.ln()) / self.sdlog)
103    }
104
105    fn quantile(&self, p: f64) -> Result<f64> {
106        check_probability(p)?;
107        Ok((self.meanlog + self.sdlog * norm_quantile(p)).exp())
108    }
109}
110
111impl Severity for Lognormal {
112    /// `E[min(X, d)] = e^(mu + s^2/2) Phi((ln d - mu - s^2) / s) + d (1 - Phi((ln d - mu) / s))`.
113    fn lev(&self, limit: f64) -> f64 {
114        if limit <= 0.0 {
115            return limit;
116        }
117        if limit == f64::INFINITY {
118            return self.mean();
119        }
120        let (mu, s) = (self.meanlog, self.sdlog);
121        let z = (limit.ln() - mu) / s;
122        self.mean() * norm_cdf(z - s) + limit * norm_cdf(-z)
123    }
124
125    /// `E[(X - d)+] = e^(mu + s^2/2) Phi((mu + s^2 - ln d) / s) - d Phi((mu - ln d) / s)`,
126    /// with both terms small in the tail, so it keeps full relative
127    /// precision where `mean() - lev(d)` would cancel.
128    fn stop_loss(&self, retention: f64) -> f64 {
129        if retention <= 0.0 {
130            return self.mean() - retention;
131        }
132        if retention == f64::INFINITY {
133            return 0.0;
134        }
135        let (mu, s) = (self.meanlog, self.sdlog);
136        let z = (retention.ln() - mu) / s;
137        (self.mean() * norm_cdf(s - z) - retention * norm_cdf(-z)).max(0.0)
138    }
139
140    /// With `b = a + limit` and `Y` the layer loss,
141    /// `E[Y^2] = E[(X - a)^2; a < X <= b] + limit^2 P(X > b)`, from the
142    /// partial moments `E[X^k; a < X <= b] = e^(k mu + k^2 s^2 / 2)
143    /// (Phi(z_b - k s) - Phi(z_a - k s))`. Differences of `Phi` are taken in
144    /// the upper tail so they keep their precision for high layers.
145    fn layer_second_moment(&self, limit: f64, attachment: f64) -> f64 {
146        let (mu, s) = (self.meanlog, self.sdlog);
147        let a = attachment.max(0.0);
148        let b = a + limit;
149        let z = |x: f64| {
150            if x <= 0.0 {
151                f64::NEG_INFINITY
152            } else {
153                (x.ln() - mu) / s
154            }
155        };
156        let (za, zb) = (z(a), z(b));
157        // Phi(hi) - Phi(lo) for hi >= lo, from the side where both are small.
158        let band = |lo: f64, hi: f64| {
159            if lo > 0.0 {
160                norm_cdf(-lo) - norm_cdf(-hi)
161            } else {
162                norm_cdf(hi) - norm_cdf(lo)
163            }
164        };
165        let m0 = band(za, zb);
166        let m1 = self.mean() * band(za - s, zb - s);
167        let m2 = (2.0 * mu + 2.0 * s * s).exp() * band(za - 2.0 * s, zb - 2.0 * s);
168        let inside = m2 - 2.0 * a * m1 + a * a * m0;
169        let above = if b == f64::INFINITY {
170            0.0
171        } else {
172            limit * limit * norm_cdf(-zb)
173        };
174        (inside + above).max(0.0)
175    }
176}
177
178#[cfg(test)]
179mod tests {
180    use super::*;
181    use prospicio_core::StreamRng;
182
183    #[test]
184    fn severity_identities_and_edges() {
185        let d = Lognormal::new(7.0, 0.5).unwrap();
186        for limit in [100.0, 1_000.0, 5_000.0] {
187            let sum = d.lev(limit) + d.stop_loss(limit);
188            assert!((sum - d.mean()).abs() < 1e-12 * d.mean());
189        }
190        assert_eq!(d.lev(0.0), 0.0);
191        assert_eq!(d.lev(-5.0), -5.0);
192        assert_eq!(d.lev(f64::INFINITY), d.mean());
193        assert_eq!(d.stop_loss(f64::INFINITY), 0.0);
194        assert_eq!(d.stop_loss(0.0), d.mean());
195        assert_eq!(d.layer(f64::INFINITY, 0.0), d.mean());
196        // Layer limits stack: 1000 xs 0 + 1000 xs 1000 = 2000 xs 0.
197        let stacked = d.layer(1_000.0, 0.0) + d.layer(1_000.0, 1_000.0);
198        assert!((stacked - d.layer(2_000.0, 0.0)).abs() < 1e-9);
199    }
200
201    #[test]
202    fn stop_loss_keeps_precision_in_the_tail() {
203        // Retention 1135 is about the 1 - 1e-12 quantile of Lognormal(0, 1).
204        let d = Lognormal::new(0.0, 1.0).unwrap();
205        let sl = d.stop_loss(1135.0);
206        // mpmath at 40 digits, closed form and quadrature agree.
207        let exact = 1.796211496502316e-10;
208        assert!((sl / exact - 1.0).abs() < 1e-12, "{sl}");
209        // The naive difference keeps only about 6 digits here (relative
210        // error 9e-7), against 2e-14 for the direct formula.
211        let naive = d.mean() - d.lev(1135.0);
212        assert!((naive / exact - 1.0).abs() > 1e-8);
213    }
214
215    #[test]
216    fn rejects_bad_parameters() {
217        assert!(Lognormal::new(0.0, 0.0).is_err());
218        assert!(Lognormal::new(f64::NAN, 1.0).is_err());
219        assert!(Lognormal::from_mean_cv(-1.0, 0.5).is_err());
220    }
221
222    #[test]
223    fn quantile_edges() {
224        let d = Lognormal::new(0.0, 1.0).unwrap();
225        assert_eq!(d.quantile(0.0), Ok(0.0));
226        assert_eq!(d.quantile(1.0), Ok(f64::INFINITY));
227        assert_eq!(d.quantile(1.1), Err(Error::InvalidProbability(1.1)));
228        assert_eq!(d.cdf(-1.0), 0.0);
229    }
230
231    /// Same draws are asserted in `python/tests` and `R/prospicio/tests`,
232    /// proving all front ends share this kernel.
233    const PINNED_SAMPLE: [f64; 3] = [1.0007760893701914, 1.6293872534754683, 1.0763869265482304];
234
235    #[test]
236    fn sample_is_pinned() {
237        let d = Lognormal::new(0.0, 1.0).unwrap();
238        assert_eq!(d.sample(&mut StreamRng::new(42, 3), 3), PINNED_SAMPLE);
239    }
240
241    #[test]
242    fn sample_is_reproducible_and_centred() {
243        let d = Lognormal::from_mean_cv(100.0, 0.3).unwrap();
244        let a = d.sample(&mut StreamRng::new(1, 0), 200_000);
245        let b = d.sample(&mut StreamRng::new(1, 0), 200_000);
246        assert_eq!(a, b);
247        let mean = a.iter().sum::<f64>() / a.len() as f64;
248        // Standard error of the mean is 30 / sqrt(200_000) ~ 0.067.
249        assert!((mean - 100.0).abs() < 0.35, "mean {mean}");
250    }
251}
prospicio_prob/lognormal.rs

prospicio_prob/
lognormal.rs
