rand_distr/
exponential.rs

1// Copyright 2018 Developers of the Rand project.
2// Copyright 2013 The Rust Project Developers.
3//
4// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
5// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
6// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
7// option. This file may not be copied, modified, or distributed
8// except according to those terms.
9
10//! The exponential distribution `Exp(λ)`.
11
12use crate::utils::ziggurat;
13use crate::{Distribution, ziggurat_tables};
14use core::fmt;
15use num_traits::Float;
16use rand::{Rng, RngExt};
17
18/// The standard exponential distribution `Exp(1)`.
19///
20/// This is equivalent to `Exp::new(1.0)` or sampling with
21/// `-rng.gen::<f64>().ln()`, but faster.
22///
23/// See [`Exp`](crate::Exp) for the general exponential distribution.
24///
25/// # Plot
26///
27/// The following plot illustrates the exponential distribution with `λ = 1`.
28///
29/// ![Exponential distribution](https://raw.githubusercontent.com/rust-random/charts/main/charts/exponential_exp1.svg)
30///
31/// # Example
32///
33/// ```
34/// use rand::prelude::*;
35/// use rand_distr::Exp1;
36///
37/// let val: f64 = rand::rng().sample(Exp1);
38/// println!("{}", val);
39/// ```
40///
41/// # Notes
42///
43/// Implemented via the ZIGNOR variant[^1] of the Ziggurat method. The exact
44/// description in the paper was adjusted to use tables for the exponential
45/// distribution rather than normal.
46///
47/// [^1]: Jurgen A. Doornik (2005). [*An Improved Ziggurat Method to
48///       Generate Normal Random Samples*](
49///       https://www.doornik.com/research/ziggurat.pdf).
50///       Nuffield College, Oxford
51#[derive(Clone, Copy, Debug)]
52#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
53pub struct Exp1;
54
55impl Distribution<f32> for Exp1 {
56    #[inline]
57    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f32 {
58        // TODO: use optimal 32-bit implementation
59        let x: f64 = self.sample(rng);
60        x as f32
61    }
62}
63
64// This could be done via `-rng.gen::<f64>().ln()` but that is slower.
65impl Distribution<f64> for Exp1 {
66    #[inline]
67    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
68        #[inline]
69        fn pdf(x: f64) -> f64 {
70            (-x).exp()
71        }
72        #[inline]
73        fn zero_case<R: Rng + ?Sized>(rng: &mut R, _u: f64) -> f64 {
74            ziggurat_tables::ZIG_EXP_R - rng.random::<f64>().ln()
75        }
76
77        ziggurat(
78            rng,
79            false,
80            &ziggurat_tables::ZIG_EXP_X,
81            &ziggurat_tables::ZIG_EXP_F,
82            pdf,
83            zero_case,
84        )
85    }
86}
87
88/// The [exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution) `Exp(λ)`.
89///
90/// The exponential distribution is a continuous probability distribution
91/// with rate parameter `λ` (`lambda`). It describes the time between events
92/// in a [`Poisson`](crate::Poisson) process, i.e. a process in which
93/// events occur continuously and independently at a constant average rate.
94///
95/// See [`Exp1`](crate::Exp1) for an optimised implementation for `λ = 1`.
96///
97/// # Density function
98///
99/// `f(x) = λ * exp(-λ * x)` for `x > 0`, when `λ > 0`.
100///
101/// For `λ = 0`, all samples yield infinity (because a Poisson process
102/// with rate 0 has no events).
103///
104/// # Plot
105///
106/// The following plot illustrates the exponential distribution with
107/// various values of `λ`.
108/// The `λ` parameter controls the rate of decay as `x` approaches infinity,
109/// and the mean of the distribution is `1/λ`.
110///
111/// ![Exponential distribution](https://raw.githubusercontent.com/rust-random/charts/main/charts/exponential.svg)
112///
113/// # Example
114///
115/// ```
116/// use rand_distr::{Exp, Distribution};
117///
118/// let exp = Exp::new(2.0).unwrap();
119/// let v = exp.sample(&mut rand::rng());
120/// println!("{} is from a Exp(2) distribution", v);
121/// ```
122#[derive(Clone, Copy, Debug, PartialEq)]
123#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
124pub struct Exp<F>
125where
126    F: Float,
127    Exp1: Distribution<F>,
128{
129    /// `lambda` stored as `1/lambda`, since this is what we scale by.
130    lambda_inverse: F,
131}
132
133/// Error type returned from [`Exp::new`].
134#[derive(Clone, Copy, Debug, PartialEq, Eq)]
135pub enum Error {
136    /// `lambda < 0` or is `-0.0` is `nan`.
137    LambdaTooSmall,
138}
139
140impl fmt::Display for Error {
141    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
142        f.write_str(match self {
143            Error::LambdaTooSmall => {
144                "lambda is negative (including -0.0) or NaN in exponential distribution"
145            }
146        })
147    }
148}
149
150#[cfg(feature = "std")]
151impl std::error::Error for Error {}
152
153impl<F: Float> Exp<F>
154where
155    F: Float,
156    Exp1: Distribution<F>,
157{
158    /// Construct a new `Exp` with the given shape parameter
159    /// `lambda`.
160    ///
161    /// # Remarks
162    ///
163    /// For custom types `N` implementing the [`Float`] trait,
164    /// the case `lambda = 0` is handled as follows: each sample corresponds
165    /// to a sample from an `Exp1` multiplied by `1 / 0`. Primitive types
166    /// yield infinity, since `1 / 0 = infinity`.
167    ///
168    /// The case `lambda = N::neg_zero()` returns an error, because `-0.0` is typically
169    /// produced through underflow of computations with a negative ideal result, and
170    /// because `1 / -0.0 = -infinity`.
171    #[inline]
172    pub fn new(lambda: F) -> Result<Exp<F>, Error> {
173        if lambda.is_sign_negative() || lambda.is_nan() {
174            return Err(Error::LambdaTooSmall);
175        }
176        Ok(Exp {
177            lambda_inverse: F::one() / lambda,
178        })
179    }
180}
181
182impl<F> Distribution<F> for Exp<F>
183where
184    F: Float,
185    Exp1: Distribution<F>,
186{
187    fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> F {
188        rng.sample(Exp1) * self.lambda_inverse
189    }
190}
191
192#[cfg(test)]
193mod test {
194    use super::*;
195
196    #[test]
197    fn test_exp() {
198        let exp = Exp::new(10.0).unwrap();
199        let mut rng = crate::test::rng(221);
200        for _ in 0..1000 {
201            assert!(exp.sample(&mut rng) >= 0.0);
202        }
203    }
204    #[test]
205    fn test_zero() {
206        let d = Exp::new(0.0).unwrap();
207        assert_eq!(d.sample(&mut crate::test::rng(21)), f64::infinity());
208    }
209    #[test]
210    #[should_panic]
211    fn test_exp_invalid_lambda_neg() {
212        Exp::new(-10.0).unwrap();
213    }
214
215    #[test]
216    #[should_panic]
217    fn test_exp_invalid_lambda_nan() {
218        Exp::new(f64::nan()).unwrap();
219    }
220
221    #[test]
222    fn exponential_distributions_can_be_compared() {
223        assert_eq!(Exp::new(1.0), Exp::new(1.0));
224    }
225}
rand_distr/exponential.rs

rand_distr/
exponential.rs
