stats_claw/distributions/positive/
weibull.rs

1//! Weibull distribution numerics, for the
2//! [`WeibullDistribution`].
3//!
4//! Equivalent to `scipy.stats.weibull_min(c = shape_parameter, scale =
5//! scale_parameter)`. The density, CDF, and quantile are closed forms; the
6//! moments use the gamma function via `exp(ln_gamma(..))`, and sampling is the
7//! inverse-CDF transform.
8//!
9//! # Examples
10//!
11//! ```
12//! use stats_claw::distributions::Cdf;
13//! use stats_claw::distributions::WeibullDistribution;
14//!
15//! // With shape=1 the Weibull is Exponential(1/scale).
16//! let d = WeibullDistribution { shape_parameter: 1.0, scale_parameter: 2.0, ..Default::default() };
17//! // cdf(2) = 1 - exp(-1) ≈ 0.6321.
18//! assert!((d.cdf(2.0) - (1.0 - (-1.0f64).exp())).abs() < 1e-12);
19//! ```
20
21use super::super::{Cdf, Moments, Pdf, Quantile, Sample};
22use crate::distributions::WeibullDistribution;
23use crate::rng::SplitMix64;
24use crate::special::ln_gamma;
25
26impl WeibullDistribution {
27    /// `(x / scale)` raised to the shape, the recurring `(x/λ)^k` term.
28    fn scaled_pow(&self, x: f64) -> f64 {
29        (x / self.scale_parameter).powf(self.shape_parameter)
30    }
31}
32
33impl Pdf for WeibullDistribution {
34    fn pdf(&self, x: f64) -> f64 {
35        if x < 0.0 {
36            return 0.0;
37        }
38        let (k, lambda) = (self.shape_parameter, self.scale_parameter);
39        let z = x / lambda;
40        (k / lambda) * z.powf(k - 1.0) * (-self.scaled_pow(x)).exp()
41    }
42}
43
44impl Cdf for WeibullDistribution {
45    fn cdf(&self, x: f64) -> f64 {
46        if x < 0.0 {
47            0.0
48        } else {
49            (-self.scaled_pow(x)).exp().mul_add(-1.0, 1.0)
50        }
51    }
52}
53
54impl Quantile for WeibullDistribution {
55    fn quantile(&self, p: f64) -> f64 {
56        self.scale_parameter * (-(1.0 - p).ln()).powf(1.0 / self.shape_parameter)
57    }
58}
59
60impl Moments for WeibullDistribution {
61    fn mean(&self) -> Option<f64> {
62        Some(self.scale_parameter * gamma(1.0 + 1.0 / self.shape_parameter))
63    }
64    fn variance(&self) -> Option<f64> {
65        let g1 = gamma(1.0 + 1.0 / self.shape_parameter);
66        let g2 = gamma(1.0 + 2.0 / self.shape_parameter);
67        Some(self.scale_parameter * self.scale_parameter * g2.mul_add(1.0, -(g1 * g1)))
68    }
69}
70
71impl Sample for WeibullDistribution {
72    fn sample(&self, rng: &mut SplitMix64) -> f64 {
73        self.quantile(rng.next_f64())
74    }
75}
76
77/// Gamma function `Γ(x)` for `x > 0`, via `exp(ln_gamma(x))`.
78fn gamma(x: f64) -> f64 {
79    ln_gamma(x).exp()
80}
81
82#[cfg(test)]
83mod tests {
84    use super::*;
85
86    /// With shape 1 the Weibull reduces to an exponential with rate `1/scale`.
87    #[test]
88    fn shape_one_is_exponential() {
89        let d = WeibullDistribution {
90            shape_parameter: 1.0,
91            scale_parameter: 2.0,
92            ..Default::default()
93        };
94        // exponential cdf at x=2 with rate 0.5 is 1 - e^{-1}
95        assert!((d.cdf(2.0) - (1.0 - (-1.0f64).exp())).abs() < 1e-12);
96    }
97
98    /// The mean is `scale · Γ(1 + 1/shape)`.
99    #[test]
100    fn mean_matches_gamma_form() {
101        let d = WeibullDistribution {
102            shape_parameter: 2.0,
103            scale_parameter: 1.0,
104            ..Default::default()
105        };
106        // Γ(1.5) = √π/2 ≈ 0.886227
107        assert!(matches!(d.mean(), Some(m) if (m - 0.886_226_925).abs() < 1e-6));
108    }
109}
stats_claw/distributions/positive/weibull.rs

stats_claw/distributions/positive/
weibull.rs
