statrs/distribution/
weibull.rs

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

statrs/distribution/
weibull.rs
