stats_claw/distributions/positive/
weibull.rs

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

stats_claw/distributions/positive/
weibull.rs
