statrs/distribution/
erlang.rs

1use crate::distribution::{Continuous, ContinuousCDF, Gamma, GammaError};
2use crate::statistics::*;
3
4/// Implements the [Erlang](https://en.wikipedia.org/wiki/Erlang_distribution)
5/// distribution
6/// which is a special case of the
7/// [Gamma](https://en.wikipedia.org/wiki/Gamma_distribution)
8/// distribution
9///
10/// # Examples
11///
12/// ```
13/// use statrs::distribution::{Erlang, Continuous};
14/// use statrs::statistics::Distribution;
15/// use statrs::prec;
16///
17/// let n = Erlang::new(3, 1.0).unwrap();
18/// assert_eq!(n.mean().unwrap(), 3.0);
19/// assert!(prec::almost_eq(n.pdf(2.0), 0.270670566473225383788, 1e-15));
20/// ```
21#[derive(Copy, Clone, PartialEq, Debug)]
22pub struct Erlang {
23    g: Gamma,
24}
25
26impl Erlang {
27    /// Constructs a new erlang distribution with a shape (k)
28    /// of `shape` and a rate (λ) of `rate`
29    ///
30    /// # Errors
31    ///
32    /// Returns an error if `shape` or `rate` are `NaN`.
33    /// Also returns an error if `shape == 0` or `rate <= 0.0`
34    ///
35    /// # Examples
36    ///
37    /// ```
38    /// use statrs::distribution::Erlang;
39    ///
40    /// let mut result = Erlang::new(3, 1.0);
41    /// assert!(result.is_ok());
42    ///
43    /// result = Erlang::new(0, 0.0);
44    /// assert!(result.is_err());
45    /// ```
46    pub fn new(shape: u64, rate: f64) -> Result<Erlang, GammaError> {
47        Gamma::new(shape as f64, rate).map(|g| Erlang { g })
48    }
49
50    /// Returns the shape (k) of the erlang distribution
51    ///
52    /// # Examples
53    ///
54    /// ```
55    /// use statrs::distribution::Erlang;
56    ///
57    /// let n = Erlang::new(3, 1.0).unwrap();
58    /// assert_eq!(n.shape(), 3);
59    /// ```
60    pub fn shape(&self) -> u64 {
61        self.g.shape() as u64
62    }
63
64    /// Returns the rate (λ) of the erlang distribution
65    ///
66    /// # Examples
67    ///
68    /// ```
69    /// use statrs::distribution::Erlang;
70    ///
71    /// let n = Erlang::new(3, 1.0).unwrap();
72    /// assert_eq!(n.rate(), 1.0);
73    /// ```
74    pub fn rate(&self) -> f64 {
75        self.g.rate()
76    }
77}
78
79impl std::fmt::Display for Erlang {
80    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
81        write!(f, "E({}, {})", self.rate(), self.shape())
82    }
83}
84
85#[cfg(feature = "rand")]
86#[cfg_attr(docsrs, doc(cfg(feature = "rand")))]
87impl ::rand::distributions::Distribution<f64> for Erlang {
88    fn sample<R: ::rand::Rng + ?Sized>(&self, rng: &mut R) -> f64 {
89        ::rand::distributions::Distribution::sample(&self.g, rng)
90    }
91}
92
93impl ContinuousCDF<f64, f64> for Erlang {
94    /// Calculates the cumulative distribution function for the erlang
95    /// distribution
96    /// at `x`
97    ///
98    /// # Formula
99    ///
100    /// ```text
101    /// γ(k, λx)  (k - 1)!
102    /// ```
103    ///
104    /// where `k` is the shape, `λ` is the rate, and `γ` is the lower
105    /// incomplete gamma function
106    fn cdf(&self, x: f64) -> f64 {
107        self.g.cdf(x)
108    }
109
110    /// Calculates the cumulative distribution function for the erlang
111    /// distribution
112    /// at `x`
113    ///
114    /// # Formula
115    ///
116    /// ```text
117    /// γ(k, λx)  (k - 1)!
118    /// ```
119    ///
120    /// where `k` is the shape, `λ` is the rate, and `γ` is the upper
121    /// incomplete gamma function
122    fn sf(&self, x: f64) -> f64 {
123        self.g.sf(x)
124    }
125
126    /// Calculates the inverse cumulative distribution function for the erlang
127    /// distribution at `x`
128    ///
129    /// # Formula
130    ///
131    /// ```text
132    /// γ^{-1}(k, (k - 1)! x) / λ
133    /// ```
134    ///
135    /// where `k` is the shape, `λ` is the rate, and `γ` is the upper
136    /// incomplete gamma function
137    fn inverse_cdf(&self, p: f64) -> f64 {
138        self.g.inverse_cdf(p)
139    }
140}
141
142impl Min<f64> for Erlang {
143    /// Returns the minimum value in the domain of the
144    /// erlang distribution representable by a double precision
145    /// float
146    ///
147    /// # Formula
148    ///
149    /// ```text
150    /// 0
151    /// ```
152    fn min(&self) -> f64 {
153        self.g.min()
154    }
155}
156
157impl Max<f64> for Erlang {
158    /// Returns the maximum value in the domain of the
159    /// erlang distribution representable by a double precision
160    /// float
161    ///
162    /// # Formula
163    ///
164    /// ```text
165    /// f64::INFINITY
166    /// ```
167    fn max(&self) -> f64 {
168        self.g.max()
169    }
170}
171
172impl Distribution<f64> for Erlang {
173    /// Returns the mean of the erlang distribution
174    ///
175    /// # Remarks
176    ///
177    /// Returns `shape` if `rate == f64::INFINITY`. This behavior
178    /// is borrowed from the Math.NET implementation
179    ///
180    /// # Formula
181    ///
182    /// ```text
183    /// k / λ
184    /// ```
185    ///
186    /// where `k` is the shape and `λ` is the rate
187    fn mean(&self) -> Option<f64> {
188        self.g.mean()
189    }
190
191    /// Returns the variance of the erlang distribution
192    ///
193    /// # Formula
194    ///
195    /// ```text
196    /// k / λ^2
197    /// ```
198    ///
199    /// where `α` is the shape and `λ` is the rate
200    fn variance(&self) -> Option<f64> {
201        self.g.variance()
202    }
203
204    /// Returns the entropy of the erlang distribution
205    ///
206    /// # Formula
207    ///
208    /// ```text
209    /// k - ln(λ) + ln(Γ(k)) + (1 - k) * ψ(k)
210    /// ```
211    ///
212    /// where `k` is the shape, `λ` is the rate, `Γ` is the gamma function,
213    /// and `ψ` is the digamma function
214    fn entropy(&self) -> Option<f64> {
215        self.g.entropy()
216    }
217
218    /// Returns the skewness of the erlang distribution
219    ///
220    /// # Formula
221    ///
222    /// ```text
223    /// 2 / sqrt(k)
224    /// ```
225    ///
226    /// where `k` is the shape
227    fn skewness(&self) -> Option<f64> {
228        self.g.skewness()
229    }
230}
231
232impl Mode<Option<f64>> for Erlang {
233    /// Returns the mode for the erlang distribution
234    ///
235    /// # Remarks
236    ///
237    /// Returns `shape` if `rate ==f64::INFINITY`. This behavior
238    /// is borrowed from the Math.NET implementation
239    ///
240    /// # Formula
241    ///
242    /// ```text
243    /// (k - 1) / λ
244    /// ```
245    ///
246    /// where `k` is the shape and `λ` is the rate
247    fn mode(&self) -> Option<f64> {
248        self.g.mode()
249    }
250}
251
252impl Continuous<f64, f64> for Erlang {
253    /// Calculates the probability density function for the erlang distribution
254    /// at `x`
255    ///
256    /// # Remarks
257    ///
258    /// Returns `NAN` if any of `shape` or `rate` are `INF`
259    /// or if `x` is `INF`
260    ///
261    /// # Formula
262    ///
263    /// ```text
264    /// (λ^k / Γ(k)) * x^(k - 1) * e^(-λ * x)
265    /// ```
266    ///
267    /// where `k` is the shape, `λ` is the rate, and `Γ` is the gamma function
268    fn pdf(&self, x: f64) -> f64 {
269        self.g.pdf(x)
270    }
271
272    /// Calculates the log probability density function for the erlang
273    /// distribution
274    /// at `x`
275    ///
276    /// # Remarks
277    ///
278    /// Returns `NAN` if any of `shape` or `rate` are `INF`
279    /// or if `x` is `INF`
280    ///
281    /// # Formula
282    ///
283    /// ```text
284    /// ln((λ^k / Γ(k)) * x^(k - 1) * e ^(-λ * x))
285    /// ```
286    ///
287    /// where `k` is the shape, `λ` is the rate, and `Γ` is the gamma function
288    fn ln_pdf(&self, x: f64) -> f64 {
289        self.g.ln_pdf(x)
290    }
291}
292
293#[rustfmt::skip]
294#[cfg(test)]
295mod tests {
296    use super::*;
297    use crate::distribution::internal::*;
298    use crate::testing_boiler;
299
300    testing_boiler!(shape: u64, rate: f64; Erlang; GammaError);
301
302    #[test]
303    fn test_create() {
304        create_ok(1, 0.1);
305        create_ok(1, 1.0);
306        create_ok(10, 10.0);
307        create_ok(10, 1.0);
308        create_ok(10, f64::INFINITY);
309    }
310
311    #[test]
312    fn test_bad_create() {
313        let invalid = [
314            (0, 1.0, GammaError::ShapeInvalid),
315            (1, 0.0, GammaError::RateInvalid),
316            (1, f64::NAN, GammaError::RateInvalid),
317            (1, -1.0, GammaError::RateInvalid),
318        ];
319
320        for (s, r, err) in invalid {
321            test_create_err(s, r, err);
322        }
323    }
324
325    #[test]
326    fn test_continuous() {
327        test::check_continuous_distribution(&create_ok(1, 2.5), 0.0, 20.0);
328        test::check_continuous_distribution(&create_ok(2, 1.5), 0.0, 20.0);
329        test::check_continuous_distribution(&create_ok(3, 0.5), 0.0, 20.0);
330    }
331}