statrs/distribution/
weibull.rs

1use crate::distribution::{Continuous, ContinuousCDF};
2use crate::function::gamma;
3use crate::is_zero;
4use crate::statistics::*;
5use crate::{consts, Result, StatsError};
6use rand::Rng;
7use std::f64;
8
9/// Implements the [Weibull](https://en.wikipedia.org/wiki/Weibull_distribution)
10/// distribution
11///
12/// # Examples
13///
14/// ```
15/// use statrs::distribution::{Weibull, Continuous};
16/// use statrs::statistics::Distribution;
17/// use statrs::prec;
18///
19/// let n = Weibull::new(10.0, 1.0).unwrap();
20/// assert!(prec::almost_eq(n.mean().unwrap(),
21/// 0.95135076986687318362924871772654021925505786260884, 1e-15));
22/// assert_eq!(n.pdf(1.0), 3.6787944117144232159552377016146086744581113103177);
23/// ```
24#[derive(Debug, Copy, Clone, PartialEq)]
25pub struct Weibull {
26    shape: f64,
27    scale: f64,
28    scale_pow_shape_inv: f64,
29}
30
31impl Weibull {
32    /// Constructs a new weibull distribution with a shape (k) of `shape`
33    /// and a scale (λ) of `scale`
34    ///
35    /// # Errors
36    ///
37    /// Returns an error if `shape` or `scale` are `NaN`.
38    /// Returns an error if `shape <= 0.0` or `scale <= 0.0`
39    ///
40    /// # Examples
41    ///
42    /// ```
43    /// use statrs::distribution::Weibull;
44    ///
45    /// let mut result = Weibull::new(10.0, 1.0);
46    /// assert!(result.is_ok());
47    ///
48    /// result = Weibull::new(0.0, 0.0);
49    /// assert!(result.is_err());
50    /// ```
51    pub fn new(shape: f64, scale: f64) -> Result<Weibull> {
52        let is_nan = shape.is_nan() || scale.is_nan();
53        match (shape, scale, is_nan) {
54            (_, _, true) => Err(StatsError::BadParams),
55            (_, _, false) if shape <= 0.0 || scale <= 0.0 => Err(StatsError::BadParams),
56            (_, _, false) => Ok(Weibull {
57                shape,
58                scale,
59                scale_pow_shape_inv: scale.powf(-shape),
60            }),
61        }
62    }
63
64    /// Returns the shape of the weibull distribution
65    ///
66    /// # Examples
67    ///
68    /// ```
69    /// use statrs::distribution::Weibull;
70    ///
71    /// let n = Weibull::new(10.0, 1.0).unwrap();
72    /// assert_eq!(n.shape(), 10.0);
73    /// ```
74    pub fn shape(&self) -> f64 {
75        self.shape
76    }
77
78    /// Returns the scale of the weibull distribution
79    ///
80    /// # Examples
81    ///
82    /// ```
83    /// use statrs::distribution::Weibull;
84    ///
85    /// let n = Weibull::new(10.0, 1.0).unwrap();
86    /// assert_eq!(n.scale(), 1.0);
87    /// ```
88    pub fn scale(&self) -> f64 {
89        self.scale
90    }
91}
92
93impl ::rand::distributions::Distribution<f64> for Weibull {
94    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
95        let x: f64 = rng.gen();
96        self.scale * (-x.ln()).powf(1.0 / self.shape)
97    }
98}
99
100impl ContinuousCDF<f64, f64> for Weibull {
101    /// Calculates the cumulative distribution function for the weibull
102    /// distribution at `x`
103    ///
104    /// # Formula
105    ///
106    /// ```ignore
107    /// 1 - e^-((x/λ)^k)
108    /// ```
109    ///
110    /// where `k` is the shape and `λ` is the scale
111    fn cdf(&self, x: f64) -> f64 {
112        if x < 0.0 {
113            0.0
114        } else {
115            -(-x.powf(self.shape) * self.scale_pow_shape_inv).exp_m1()
116        }
117    }
118
119    /// Calculates the survival function for the weibull
120    /// distribution at `x`
121    ///
122    /// # Formula
123    ///
124    /// ```ignore
125    /// e^-((x/λ)^k)
126    /// ```
127    ///
128    /// where `k` is the shape and `λ` is the scale
129    fn sf(&self, x: f64) -> f64 {
130        if x < 0.0 {
131            1.0
132        } else {
133            (-x.powf(self.shape) * self.scale_pow_shape_inv).exp()
134        }
135    }
136}
137
138impl Min<f64> for Weibull {
139    /// Returns the minimum value in the domain of the weibull
140    /// distribution representable by a double precision float
141    ///
142    /// # Formula
143    ///
144    /// ```ignore
145    /// 0
146    /// ```
147    fn min(&self) -> f64 {
148        0.0
149    }
150}
151
152impl Max<f64> for Weibull {
153    /// Returns the maximum value in the domain of the weibull
154    /// distribution representable by a double precision float
155    ///
156    /// # Formula
157    ///
158    /// ```ignore
159    /// INF
160    /// ```
161    fn max(&self) -> f64 {
162        f64::INFINITY
163    }
164}
165
166impl Distribution<f64> for Weibull {
167    /// Returns the mean of the weibull distribution
168    ///
169    /// # Formula
170    ///
171    /// ```ignore
172    /// λΓ(1 + 1 / k)
173    /// ```
174    ///
175    /// where `k` is the shape, `λ` is the scale, and `Γ` is
176    /// the gamma function
177    fn mean(&self) -> Option<f64> {
178        Some(self.scale * gamma::gamma(1.0 + 1.0 / self.shape))
179    }
180    /// Returns the variance of the weibull distribution
181    ///
182    /// # Formula
183    ///
184    /// ```ignore
185    /// λ^2 * (Γ(1 + 2 / k) - Γ(1 + 1 / k)^2)
186    /// ```
187    ///
188    /// where `k` is the shape, `λ` is the scale, and `Γ` is
189    /// the gamma function
190    fn variance(&self) -> Option<f64> {
191        let mean = self.mean()?;
192        Some(self.scale * self.scale * gamma::gamma(1.0 + 2.0 / self.shape) - mean * mean)
193    }
194    /// Returns the entropy of the weibull distribution
195    ///
196    /// # Formula
197    ///
198    /// ```ignore
199    /// γ(1 - 1 / k) + ln(λ / k) + 1
200    /// ```
201    ///
202    /// where `k` is the shape, `λ` is the scale, and `γ` is
203    /// the Euler-Mascheroni constant
204    fn entropy(&self) -> Option<f64> {
205        let entr = consts::EULER_MASCHERONI * (1.0 - 1.0 / self.shape)
206            + (self.scale / self.shape).ln()
207            + 1.0;
208        Some(entr)
209    }
210    /// Returns the skewness of the weibull distribution
211    ///
212    /// # Formula
213    ///
214    /// ```ignore
215    /// (Γ(1 + 3 / k) * λ^3 - 3μσ^2 - μ^3) / σ^3
216    /// ```
217    ///
218    /// where `k` is the shape, `λ` is the scale, and `Γ` is
219    /// the gamma function, `μ` is the mean of the distribution.
220    /// and `σ` the standard deviation of the distribution
221    fn skewness(&self) -> Option<f64> {
222        let mu = self.mean()?;
223        let sigma = self.std_dev()?;
224        let sigma2 = sigma * sigma;
225        let sigma3 = sigma2 * sigma;
226        let skew = (self.scale * self.scale * self.scale * gamma::gamma(1.0 + 3.0 / self.shape)
227            - 3.0 * sigma2 * mu
228            - (mu * mu * mu))
229            / sigma3;
230        Some(skew)
231    }
232}
233
234impl Median<f64> for Weibull {
235    /// Returns the median of the weibull distribution
236    ///
237    /// # Formula
238    ///
239    /// ```ignore
240    /// λ(ln(2))^(1 / k)
241    /// ```
242    ///
243    /// where `k` is the shape and `λ` is the scale
244    fn median(&self) -> f64 {
245        self.scale * f64::consts::LN_2.powf(1.0 / self.shape)
246    }
247}
248
249impl Mode<Option<f64>> for Weibull {
250    /// Returns the median of the weibull distribution
251    ///
252    /// # Formula
253    ///
254    /// ```ignore
255    /// if k == 1 {
256    ///     0
257    /// } else {
258    ///     λ((k - 1) / k)^(1 / k)
259    /// }
260    /// ```
261    ///
262    /// where `k` is the shape and `λ` is the scale
263    fn mode(&self) -> Option<f64> {
264        let mode = if ulps_eq!(self.shape, 1.0) {
265            0.0
266        } else {
267            self.scale * ((self.shape - 1.0) / self.shape).powf(1.0 / self.shape)
268        };
269        Some(mode)
270    }
271}
272
273impl Continuous<f64, f64> for Weibull {
274    /// Calculates the probability density function for the weibull
275    /// distribution at `x`
276    ///
277    /// # Formula
278    ///
279    /// ```ignore
280    /// (k / λ) * (x / λ)^(k - 1) * e^(-(x / λ)^k)
281    /// ```
282    ///
283    /// where `k` is the shape and `λ` is the scale
284    fn pdf(&self, x: f64) -> f64 {
285        if x < 0.0 {
286            0.0
287        } else if is_zero(x) && ulps_eq!(self.shape, 1.0) {
288            1.0 / self.scale
289        } else if x.is_infinite() {
290            0.0
291        } else {
292            self.shape
293                * (x / self.scale).powf(self.shape - 1.0)
294                * (-(x.powf(self.shape)) * self.scale_pow_shape_inv).exp()
295                / self.scale
296        }
297    }
298
299    /// Calculates the log probability density function for the weibull
300    /// distribution at `x`
301    ///
302    /// # Formula
303    ///
304    /// ```ignore
305    /// ln((k / λ) * (x / λ)^(k - 1) * e^(-(x / λ)^k))
306    /// ```
307    ///
308    /// where `k` is the shape and `λ` is the scale
309    fn ln_pdf(&self, x: f64) -> f64 {
310        if x < 0.0 {
311            f64::NEG_INFINITY
312        } else if is_zero(x) && ulps_eq!(self.shape, 1.0) {
313            0.0 - self.scale.ln()
314        } else if x.is_infinite() {
315            f64::NEG_INFINITY
316        } else {
317            self.shape.ln() + (self.shape - 1.0) * (x / self.scale).ln()
318                - x.powf(self.shape) * self.scale_pow_shape_inv
319                - self.scale.ln()
320        }
321    }
322}
323
324#[rustfmt::skip]
325#[cfg(all(test, feature = "nightly"))]
326mod tests {
327    use crate::statistics::*;
328    use crate::distribution::{ContinuousCDF, Continuous, Weibull};
329    use crate::distribution::internal::*;
330    use crate::consts::ACC;
331
332    fn try_create(shape: f64, scale: f64) -> Weibull {
333        let n = Weibull::new(shape, scale);
334        assert!(n.is_ok());
335        n.unwrap()
336    }
337
338    fn create_case(shape: f64, scale: f64) {
339        let n = try_create(shape, scale);
340        assert_eq!(shape, n.shape());
341        assert_eq!(scale, n.scale());
342    }
343
344    fn bad_create_case(shape: f64, scale: f64) {
345        let n = Weibull::new(shape, scale);
346        assert!(n.is_err());
347    }
348
349    fn get_value<F>(shape: f64, scale: f64, eval: F) -> f64
350        where F: Fn(Weibull) -> f64
351    {
352        let n = try_create(shape, scale);
353        eval(n)
354    }
355
356    fn test_case<F>(shape: f64, scale: f64, expected: f64, eval: F)
357        where F: Fn(Weibull) -> f64
358    {
359        let x = get_value(shape, scale, eval);
360        assert_eq!(expected, x);
361    }
362
363    fn test_almost<F>(shape: f64, scale: f64, expected: f64, acc: f64, eval: F)
364        where F: Fn(Weibull) -> f64
365    {
366        let x = get_value(shape, scale, eval);
367        assert_almost_eq!(expected, x, acc);
368    }
369
370    #[test]
371    fn test_create() {
372        create_case(1.0, 0.1);
373        create_case(10.0, 1.0);
374        create_case(11.0, 10.0);
375        create_case(12.0, f64::INFINITY);
376    }
377
378    #[test]
379    fn test_bad_create() {
380        bad_create_case(f64::NAN, 1.0);
381        bad_create_case(1.0, f64::NAN);
382        bad_create_case(f64::NAN, f64::NAN);
383        bad_create_case(1.0, -1.0);
384        bad_create_case(-1.0, 1.0);
385        bad_create_case(-1.0, -1.0);
386        bad_create_case(0.0, 0.0);
387        bad_create_case(0.0, 1.0);
388        bad_create_case(1.0, 0.0);
389    }
390
391    #[test]
392    fn test_mean() {
393        let mean = |x: Weibull| x.mean().unwrap();
394        test_case(1.0, 0.1, 0.1, mean);
395        test_case(1.0, 1.0, 1.0, mean);
396        test_almost(10.0, 10.0, 9.5135076986687318362924871772654021925505786260884, 1e-14, mean);
397        test_almost(10.0, 1.0, 0.95135076986687318362924871772654021925505786260884, 1e-15, mean);
398    }
399
400    #[test]
401    fn test_variance() {
402        let variance = |x: Weibull| x.variance().unwrap();
403        test_almost(1.0, 0.1, 0.01, 1e-16, variance);
404        test_almost(1.0, 1.0, 1.0, 1e-14, variance);
405        test_almost(10.0, 10.0, 1.3100455073468309147154581687505295026863354547057, 1e-12, variance);
406        test_almost(10.0, 1.0, 0.013100455073468309147154581687505295026863354547057, 1e-14, variance);
407    }
408
409    #[test]
410    fn test_entropy() {
411        let entropy = |x: Weibull| x.entropy().unwrap();
412        test_almost(1.0, 0.1, -1.302585092994045684018, 1e-15, entropy);
413        test_case(1.0, 1.0, 1.0, entropy);
414        test_case(10.0, 10.0, 1.519494098411379574546, entropy);
415        test_almost(10.0, 1.0, -0.783090994582666109472, 1e-15, entropy);
416    }
417
418    #[test]
419    fn test_skewnewss() {
420        let skewness = |x: Weibull| x.skewness().unwrap();
421        test_almost(1.0, 0.1, 2.0, 1e-13, skewness);
422        test_almost(1.0, 1.0, 2.0, 1e-13, skewness);
423        test_almost(10.0, 10.0, -0.63763713390314440916597757156663888653981696212127, 1e-11, skewness);
424        test_almost(10.0, 1.0, -0.63763713390314440916597757156663888653981696212127, 1e-11, skewness);
425    }
426
427    #[test]
428    fn test_median() {
429        let median = |x: Weibull| x.median();
430        test_case(1.0, 0.1, 0.069314718055994530941723212145817656807550013436026, median);
431        test_case(1.0, 1.0, 0.69314718055994530941723212145817656807550013436026, median);
432        test_case(10.0, 10.0, 9.6401223546778973665856033763604752124634905617583, median);
433        test_case(10.0, 1.0, 0.96401223546778973665856033763604752124634905617583, median);
434    }
435
436    #[test]
437    fn test_mode() {
438        let mode = |x: Weibull| x.mode().unwrap();
439        test_case(1.0, 0.1, 0.0, mode);
440        test_case(1.0, 1.0, 0.0, mode);
441        test_case(10.0, 10.0, 9.8951925820621439264623017041980483215553841533709, mode);
442        test_case(10.0, 1.0, 0.98951925820621439264623017041980483215553841533709, mode);
443    }
444
445    #[test]
446    fn test_min_max() {
447        let min = |x: Weibull| x.min();
448        let max = |x: Weibull| x.max();
449        test_case(1.0, 1.0, 0.0, min);
450        test_case(1.0, 1.0, f64::INFINITY, max);
451    }
452
453    #[test]
454    fn test_pdf() {
455        let pdf = |arg: f64| move |x: Weibull| x.pdf(arg);
456        test_case(1.0, 0.1, 10.0, pdf(0.0));
457        test_case(1.0, 0.1, 0.00045399929762484851535591515560550610237918088866565, pdf(1.0));
458        test_case(1.0, 0.1, 3.7200759760208359629596958038631183373588922923768e-43, pdf(10.0));
459        test_case(1.0, 1.0, 1.0, pdf(0.0));
460        test_case(1.0, 1.0, 0.36787944117144232159552377016146086744581113103177, pdf(1.0));
461        test_case(1.0, 1.0, 0.000045399929762484851535591515560550610237918088866565, pdf(10.0));
462        test_case(10.0, 10.0, 0.0, pdf(0.0));
463        test_almost(10.0, 10.0, 9.9999999990000000000499999999983333333333750000000e-10, 1e-24, pdf(1.0));
464        test_case(10.0, 10.0, 0.36787944117144232159552377016146086744581113103177, pdf(10.0));
465        test_case(10.0, 1.0, 0.0, pdf(0.0));
466        test_case(10.0, 1.0, 3.6787944117144232159552377016146086744581113103177, pdf(1.0));
467        test_case(10.0, 1.0, 0.0, pdf(10.0));
468    }
469
470    #[test]
471    fn test_ln_pdf() {
472        let ln_pdf = |arg: f64| move |x: Weibull| x.ln_pdf(arg);
473        test_almost(1.0, 0.1, 2.3025850929940456840179914546843642076011014886288, 1e-15, ln_pdf(0.0));
474        test_almost(1.0, 0.1, -7.6974149070059543159820085453156357923988985113712, 1e-15, ln_pdf(1.0));
475        test_case(1.0, 0.1, -97.697414907005954315982008545315635792398898511371, ln_pdf(10.0));
476        test_case(1.0, 1.0, 0.0, ln_pdf(0.0));
477        test_case(1.0, 1.0, -1.0, ln_pdf(1.0));
478        test_case(1.0, 1.0, -10.0, ln_pdf(10.0));
479        test_case(10.0, 10.0, f64::NEG_INFINITY, ln_pdf(0.0));
480        test_almost(10.0, 10.0, -20.723265837046411156161923092159277868409913397659, 1e-14, ln_pdf(1.0));
481        test_case(10.0, 10.0, -1.0, ln_pdf(10.0));
482        test_case(10.0, 1.0, f64::NEG_INFINITY, ln_pdf(0.0));
483        test_almost(10.0, 1.0, 1.3025850929940456840179914546843642076011014886288, 1e-15, ln_pdf(1.0));
484        test_case(10.0, 1.0, -9.999999976974149070059543159820085453156357923988985113712e9, ln_pdf(10.0));
485    }
486
487    #[test]
488    fn test_cdf() {
489        let cdf = |arg: f64| move |x: Weibull| x.cdf(arg);
490        test_case(1.0, 0.1, 0.0, cdf(0.0));
491        test_case(1.0, 0.1, 0.99995460007023751514846440848443944938976208191113, cdf(1.0));
492        test_case(1.0, 0.1, 0.99999999999999999999999999999999999999999996279924, cdf(10.0));
493        test_case(1.0, 1.0, 0.0, cdf(0.0));
494        test_case(1.0, 1.0, 0.63212055882855767840447622983853913255418886896823, cdf(1.0));
495        test_case(1.0, 1.0, 0.99995460007023751514846440848443944938976208191113, cdf(10.0));
496        test_case(10.0, 10.0, 0.0, cdf(0.0));
497        test_almost(10.0, 10.0, 9.9999999995000000000166666666662500000000083333333e-11, 1e-25, cdf(1.0));
498        test_case(10.0, 10.0, 0.63212055882855767840447622983853913255418886896823, cdf(10.0));
499        test_case(10.0, 1.0, 0.0, cdf(0.0));
500        test_case(10.0, 1.0, 0.63212055882855767840447622983853913255418886896823, cdf(1.0));
501        test_case(10.0, 1.0, 1.0, cdf(10.0));
502    }
503
504    #[test]
505    fn test_sf() {
506        let sf = |arg: f64| move |x: Weibull| x.sf(arg);
507        test_case(1.0, 0.1, 1.0, sf(0.0));
508        test_case(1.0, 0.1, 4.5399929762484854e-5, sf(1.0));
509        test_case(1.0, 0.1, 3.720075976020836e-44, sf(10.0));
510        test_case(1.0, 1.0, 1.0, sf(0.0));
511        test_case(1.0, 1.0, 0.36787944117144233, sf(1.0));
512        test_case(1.0, 1.0, 4.5399929762484854e-5, sf(10.0));
513        test_case(10.0, 10.0, 1.0, sf(0.0));
514        test_almost(10.0, 10.0, 0.9999999999, 1e-25, sf(1.0));
515        test_case(10.0, 10.0, 0.36787944117144233, sf(10.0));
516        test_case(10.0, 1.0, 1.0, sf(0.0));
517        test_case(10.0, 1.0, 0.36787944117144233, sf(1.0));
518        test_case(10.0, 1.0, 0.0, sf(10.0));
519    }
520
521    #[test]
522    fn test_continuous() {
523        test::check_continuous_distribution(&try_create(1.0, 0.2), 0.0, 10.0);
524    }
525}
statrs/distribution/weibull.rs

statrs/distribution/
weibull.rs
