rand_distr/
weibull.rs

1// Copyright 2018 Developers of the Rand project.
2//
3// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
4// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
5// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
6// option. This file may not be copied, modified, or distributed
7// except according to those terms.
8
9//! The Weibull distribution `Weibull(λ, k)`
10
11use crate::{Distribution, OpenClosed01};
12use core::fmt;
13use num_traits::Float;
14use rand::Rng;
15
16/// The [Weibull distribution](https://en.wikipedia.org/wiki/Weibull_distribution) `Weibull(λ, k)`.
17///
18/// This is a family of continuous probability distributions with
19/// scale parameter `λ` (`lambda`) and shape parameter `k`. It is used
20/// to model reliability data, life data, and accelerated life testing data.
21///
22/// # Density function
23///
24/// `f(x; λ, k) = (k / λ) * (x / λ)^(k - 1) * exp(-(x / λ)^k)` for `x >= 0`.
25///
26/// # Plot
27///
28/// The following plot shows the Weibull distribution with various values of `λ` and `k`.
29///
30/// ![Weibull distribution](https://raw.githubusercontent.com/rust-random/charts/main/charts/weibull.svg)
31///
32/// # Example
33/// ```
34/// use rand::prelude::*;
35/// use rand_distr::Weibull;
36///
37/// let val: f64 = rand::rng().sample(Weibull::new(1., 10.).unwrap());
38/// println!("{}", val);
39/// ```
40///
41/// # Numerics
42///
43/// For small `k` like `< 0.005`, even with `f64` a significant number of samples will be so small that they underflow to `0.0`
44/// or so big they overflow to `inf`. This is a limitation of the floating point representation and not specific to this implementation.
45#[derive(Clone, Copy, Debug, PartialEq)]
46#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
47pub struct Weibull<F>
48where
49    F: Float,
50    OpenClosed01: Distribution<F>,
51{
52    inv_shape: F,
53    scale: F,
54}
55
56/// Error type returned from [`Weibull::new`].
57#[derive(Clone, Copy, Debug, PartialEq, Eq)]
58pub enum Error {
59    /// `scale <= 0` or `nan`.
60    ScaleTooSmall,
61    /// `shape <= 0` or `nan`.
62    ShapeTooSmall,
63}
64
65impl fmt::Display for Error {
66    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
67        f.write_str(match self {
68            Error::ScaleTooSmall => "scale is not positive in Weibull distribution",
69            Error::ShapeTooSmall => "shape is not positive in Weibull distribution",
70        })
71    }
72}
73
74#[cfg(feature = "std")]
75impl std::error::Error for Error {}
76
77impl<F> Weibull<F>
78where
79    F: Float,
80    OpenClosed01: Distribution<F>,
81{
82    /// Construct a new `Weibull` distribution with given `scale` and `shape`.
83    pub fn new(scale: F, shape: F) -> Result<Weibull<F>, Error> {
84        if !(scale > F::zero()) {
85            return Err(Error::ScaleTooSmall);
86        }
87        if !(shape > F::zero()) {
88            return Err(Error::ShapeTooSmall);
89        }
90        Ok(Weibull {
91            inv_shape: F::from(1.).unwrap() / shape,
92            scale,
93        })
94    }
95}
96
97impl<F> Distribution<F> for Weibull<F>
98where
99    F: Float,
100    OpenClosed01: Distribution<F>,
101{
102    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> F {
103        let x: F = rng.sample(OpenClosed01);
104        self.scale * (-x.ln()).powf(self.inv_shape)
105    }
106}
107
108#[cfg(test)]
109mod tests {
110    use super::*;
111
112    #[test]
113    #[should_panic]
114    fn invalid() {
115        Weibull::new(0., 0.).unwrap();
116    }
117
118    #[test]
119    fn sample() {
120        let scale = 1.0;
121        let shape = 2.0;
122        let d = Weibull::new(scale, shape).unwrap();
123        let mut rng = crate::test::rng(1);
124        for _ in 0..1000 {
125            let r = d.sample(&mut rng);
126            assert!(r >= 0.);
127        }
128    }
129
130    #[test]
131    fn value_stability() {
132        fn test_samples<F: Float + fmt::Debug, D: Distribution<F>>(
133            distr: D,
134            zero: F,
135            expected: &[F],
136        ) {
137            let mut rng = crate::test::rng(213);
138            let mut buf = [zero; 4];
139            for x in &mut buf {
140                *x = rng.sample(&distr);
141            }
142            assert_eq!(buf, expected);
143        }
144
145        test_samples(
146            Weibull::new(1.0, 1.0).unwrap(),
147            0f32,
148            &[0.041495778, 0.7531094, 1.4189332, 0.38386202],
149        );
150        test_samples(
151            Weibull::new(2.0, 0.5).unwrap(),
152            0f64,
153            &[
154                1.1343478702739669,
155                0.29470010050655226,
156                0.7556151370284702,
157                7.877212340241561,
158            ],
159        );
160    }
161
162    #[test]
163    fn weibull_distributions_can_be_compared() {
164        assert_eq!(Weibull::new(1.0, 2.0), Weibull::new(1.0, 2.0));
165    }
166}