statrs/distribution/
exponential.rs

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

statrs/distribution/
exponential.rs
