statrs/distribution/
poisson.rs

1use crate::distribution::{Discrete, DiscreteCDF};
2use crate::function::{factorial, gamma};
3use crate::statistics::*;
4use std::f64;
5
6/// Implements the [Poisson](https://en.wikipedia.org/wiki/Poisson_distribution)
7/// distribution
8///
9/// # Examples
10///
11/// ```
12/// use statrs::distribution::{Poisson, Discrete};
13/// use statrs::statistics::Distribution;
14/// use statrs::prec;
15///
16/// let n = Poisson::new(1.0).unwrap();
17/// assert_eq!(n.mean().unwrap(), 1.0);
18/// assert!(prec::almost_eq(n.pmf(1), 0.367879441171442, 1e-15));
19/// ```
20#[derive(Copy, Clone, PartialEq, Debug)]
21pub struct Poisson {
22    lambda: f64,
23}
24
25/// Represents the errors that can occur when creating a [`Poisson`].
26#[derive(Copy, Clone, PartialEq, Eq, Debug, Hash)]
27#[non_exhaustive]
28pub enum PoissonError {
29    /// The lambda is NaN, zero or less than zero.
30    LambdaInvalid,
31}
32
33impl std::fmt::Display for PoissonError {
34    #[cfg_attr(coverage_nightly, coverage(off))]
35    fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
36        match self {
37            PoissonError::LambdaInvalid => write!(f, "Lambda is NaN, zero or less than zero"),
38        }
39    }
40}
41
42impl std::error::Error for PoissonError {}
43
44impl Poisson {
45    /// Constructs a new poisson distribution with a rate (λ)
46    /// of `lambda`
47    ///
48    /// # Errors
49    ///
50    /// Returns an error if `lambda` is `NaN` or `lambda <= 0.0`
51    ///
52    /// # Examples
53    ///
54    /// ```
55    /// use statrs::distribution::Poisson;
56    ///
57    /// let mut result = Poisson::new(1.0);
58    /// assert!(result.is_ok());
59    ///
60    /// result = Poisson::new(0.0);
61    /// assert!(result.is_err());
62    /// ```
63    pub fn new(lambda: f64) -> Result<Poisson, PoissonError> {
64        if lambda.is_nan() || lambda <= 0.0 {
65            Err(PoissonError::LambdaInvalid)
66        } else {
67            Ok(Poisson { lambda })
68        }
69    }
70
71    /// Returns the rate (λ) of the poisson distribution
72    ///
73    /// # Examples
74    ///
75    /// ```
76    /// use statrs::distribution::Poisson;
77    ///
78    /// let n = Poisson::new(1.0).unwrap();
79    /// assert_eq!(n.lambda(), 1.0);
80    /// ```
81    pub fn lambda(&self) -> f64 {
82        self.lambda
83    }
84}
85
86impl std::fmt::Display for Poisson {
87    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
88        write!(f, "Pois({})", self.lambda)
89    }
90}
91
92#[cfg(feature = "rand")]
93#[cfg_attr(docsrs, doc(cfg(feature = "rand")))]
94impl ::rand::distributions::Distribution<u64> for Poisson {
95    /// Generates one sample from the Poisson distribution either by
96    /// Knuth's method if lambda < 30.0 or Rejection method PA by
97    /// A. C. Atkinson from the Journal of the Royal Statistical Society
98    /// Series C (Applied Statistics) Vol. 28 No. 1. (1979) pp. 29 - 35
99    /// otherwise
100    fn sample<R: ::rand::Rng + ?Sized>(&self, rng: &mut R) -> u64 {
101        sample_unchecked(rng, self.lambda) as u64
102    }
103}
104
105#[cfg(feature = "rand")]
106#[cfg_attr(docsrs, doc(cfg(feature = "rand")))]
107impl ::rand::distributions::Distribution<f64> for Poisson {
108    /// Generates one sample from the Poisson distribution either by
109    /// Knuth's method if lambda < 30.0 or Rejection method PA by
110    /// A. C. Atkinson from the Journal of the Royal Statistical Society
111    /// Series C (Applied Statistics) Vol. 28 No. 1. (1979) pp. 29 - 35
112    /// otherwise
113    fn sample<R: ::rand::Rng + ?Sized>(&self, rng: &mut R) -> f64 {
114        sample_unchecked(rng, self.lambda)
115    }
116}
117
118impl DiscreteCDF<u64, f64> for Poisson {
119    /// Calculates the cumulative distribution function for the poisson
120    /// distribution at `x`
121    ///
122    /// # Formula
123    ///
124    /// ```text
125    /// P(x + 1, λ)
126    /// ```
127    ///
128    /// where `λ` is the rate and `P` is the upper regularized gamma function
129    fn cdf(&self, x: u64) -> f64 {
130        gamma::gamma_ur(x as f64 + 1.0, self.lambda)
131    }
132
133    /// Calculates the survival function for the poisson
134    /// distribution at `x`
135    ///
136    /// # Formula
137    ///
138    /// ```text
139    /// P(x + 1, λ)
140    /// ```
141    ///
142    /// where `λ` is the rate and `P` is the lower regularized gamma function
143    fn sf(&self, x: u64) -> f64 {
144        gamma::gamma_lr(x as f64 + 1.0, self.lambda)
145    }
146}
147
148impl Min<u64> for Poisson {
149    /// Returns the minimum value in the domain of the poisson distribution
150    /// representable by a 64-bit integer
151    ///
152    /// # Formula
153    ///
154    /// ```text
155    /// 0
156    /// ```
157    fn min(&self) -> u64 {
158        0
159    }
160}
161
162impl Max<u64> for Poisson {
163    /// Returns the maximum value in the domain of the poisson distribution
164    /// representable by a 64-bit integer
165    ///
166    /// # Formula
167    ///
168    /// ```text
169    /// 2^63 - 1
170    /// ```
171    fn max(&self) -> u64 {
172        u64::MAX
173    }
174}
175
176impl Distribution<f64> for Poisson {
177    /// Returns the mean of the poisson distribution
178    ///
179    /// # Formula
180    ///
181    /// ```text
182    /// λ
183    /// ```
184    ///
185    /// where `λ` is the rate
186    fn mean(&self) -> Option<f64> {
187        Some(self.lambda)
188    }
189
190    /// Returns the variance of the poisson distribution
191    ///
192    /// # Formula
193    ///
194    /// ```text
195    /// λ
196    /// ```
197    ///
198    /// where `λ` is the rate
199    fn variance(&self) -> Option<f64> {
200        Some(self.lambda)
201    }
202
203    /// Returns the entropy of the poisson distribution
204    ///
205    /// # Formula
206    ///
207    /// ```text
208    /// (1 / 2) * ln(2πeλ) - 1 / (12λ) - 1 / (24λ^2) - 19 / (360λ^3)
209    /// ```
210    ///
211    /// where `λ` is the rate
212    fn entropy(&self) -> Option<f64> {
213        Some(
214            0.5 * (2.0 * f64::consts::PI * f64::consts::E * self.lambda).ln()
215                - 1.0 / (12.0 * self.lambda)
216                - 1.0 / (24.0 * self.lambda * self.lambda)
217                - 19.0 / (360.0 * self.lambda * self.lambda * self.lambda),
218        )
219    }
220
221    /// Returns the skewness of the poisson distribution
222    ///
223    /// # Formula
224    ///
225    /// ```text
226    /// λ^(-1/2)
227    /// ```
228    ///
229    /// where `λ` is the rate
230    fn skewness(&self) -> Option<f64> {
231        Some(1.0 / self.lambda.sqrt())
232    }
233}
234
235impl Median<f64> for Poisson {
236    /// Returns the median of the poisson distribution
237    ///
238    /// # Formula
239    ///
240    /// ```text
241    /// floor(λ + 1 / 3 - 0.02 / λ)
242    /// ```
243    ///
244    /// where `λ` is the rate
245    fn median(&self) -> f64 {
246        (self.lambda + 1.0 / 3.0 - 0.02 / self.lambda).floor()
247    }
248}
249
250impl Mode<Option<u64>> for Poisson {
251    /// Returns the mode of the poisson distribution
252    ///
253    /// # Formula
254    ///
255    /// ```text
256    /// floor(λ)
257    /// ```
258    ///
259    /// where `λ` is the rate
260    fn mode(&self) -> Option<u64> {
261        Some(self.lambda.floor() as u64)
262    }
263}
264
265impl Discrete<u64, f64> for Poisson {
266    /// Calculates the probability mass function for the poisson distribution at
267    /// `x`
268    ///
269    /// # Formula
270    ///
271    /// ```text
272    /// (λ^x * e^(-λ)) / x!
273    /// ```
274    ///
275    /// where `λ` is the rate
276    fn pmf(&self, x: u64) -> f64 {
277        (-self.lambda + x as f64 * self.lambda.ln() - factorial::ln_factorial(x)).exp()
278    }
279
280    /// Calculates the log probability mass function for the poisson
281    /// distribution at
282    /// `x`
283    ///
284    /// # Formula
285    ///
286    /// ```text
287    /// ln((λ^x * e^(-λ)) / x!)
288    /// ```
289    ///
290    /// where `λ` is the rate
291    fn ln_pmf(&self, x: u64) -> f64 {
292        -self.lambda + x as f64 * self.lambda.ln() - factorial::ln_factorial(x)
293    }
294}
295
296/// Generates one sample from the Poisson distribution either by
297/// Knuth's method if lambda < 30.0 or Rejection method PA by
298/// A. C. Atkinson from the Journal of the Royal Statistical Society
299/// Series C (Applied Statistics) Vol. 28 No. 1. (1979) pp. 29 - 35
300/// otherwise
301#[cfg(feature = "rand")]
302#[cfg_attr(docsrs, doc(cfg(feature = "rand")))]
303pub fn sample_unchecked<R: ::rand::Rng + ?Sized>(rng: &mut R, lambda: f64) -> f64 {
304    if lambda < 30.0 {
305        let limit = (-lambda).exp();
306        let mut count = 0.0;
307        let mut product: f64 = rng.gen();
308        while product >= limit {
309            count += 1.0;
310            product *= rng.gen::<f64>();
311        }
312        count
313    } else {
314        let c = 0.767 - 3.36 / lambda;
315        let beta = f64::consts::PI / (3.0 * lambda).sqrt();
316        let alpha = beta * lambda;
317        let k = c.ln() - lambda - beta.ln();
318
319        loop {
320            let u: f64 = rng.gen();
321            let x = (alpha - ((1.0 - u) / u).ln()) / beta;
322            let n = (x + 0.5).floor();
323            if n < 0.0 {
324                continue;
325            }
326
327            let v: f64 = rng.gen();
328            let y = alpha - beta * x;
329            let temp = 1.0 + y.exp();
330            let lhs = y + (v / (temp * temp)).ln();
331            let rhs = k + n * lambda.ln() - factorial::ln_factorial(n as u64);
332            if lhs <= rhs {
333                return n;
334            }
335        }
336    }
337}
338
339#[rustfmt::skip]
340#[cfg(test)]
341mod tests {
342    use super::*;
343    use crate::distribution::internal::*;
344    use crate::testing_boiler;
345
346    testing_boiler!(lambda: f64; Poisson; PoissonError);
347
348    #[test]
349    fn test_create() {
350        create_ok(1.5);
351        create_ok(5.4);
352        create_ok(10.8);
353    }
354
355    #[test]
356    fn test_bad_create() {
357        create_err(f64::NAN);
358        create_err(-1.5);
359        create_err(0.0);
360    }
361
362    #[test]
363    fn test_mean() {
364        let mean = |x: Poisson| x.mean().unwrap();
365        test_exact(1.5, 1.5, mean);
366        test_exact(5.4, 5.4, mean);
367        test_exact(10.8, 10.8, mean);
368    }
369
370    #[test]
371    fn test_variance() {
372        let variance = |x: Poisson| x.variance().unwrap();
373        test_exact(1.5, 1.5, variance);
374        test_exact(5.4, 5.4, variance);
375        test_exact(10.8, 10.8, variance);
376    }
377
378    #[test]
379    fn test_entropy() {
380        let entropy = |x: Poisson| x.entropy().unwrap();
381        test_absolute(1.5, 1.531959153102376331946, 1e-15, entropy);
382        test_absolute(5.4, 2.244941839577643504608, 1e-15, entropy);
383        test_exact(10.8, 2.600596429676975222694, entropy);
384    }
385
386    #[test]
387    fn test_skewness() {
388        let skewness = |x: Poisson| x.skewness().unwrap();
389        test_absolute(1.5, 0.8164965809277260327324, 1e-15, skewness);
390        test_absolute(5.4, 0.4303314829119352094644, 1e-16, skewness);
391        test_absolute(10.8, 0.3042903097250922852539, 1e-16, skewness);
392    }
393
394    #[test]
395    fn test_median() {
396        let median = |x: Poisson| x.median();
397        test_exact(1.5, 1.0, median);
398        test_exact(5.4, 5.0, median);
399        test_exact(10.8, 11.0, median);
400    }
401
402    #[test]
403    fn test_mode() {
404        let mode = |x: Poisson| x.mode().unwrap();
405        test_exact(1.5, 1, mode);
406        test_exact(5.4, 5, mode);
407        test_exact(10.8, 10, mode);
408    }
409
410    #[test]
411    fn test_min_max() {
412        let min = |x: Poisson| x.min();
413        let max = |x: Poisson| x.max();
414        test_exact(1.5, 0, min);
415        test_exact(5.4, 0, min);
416        test_exact(10.8, 0, min);
417        test_exact(1.5, u64::MAX, max);
418        test_exact(5.4, u64::MAX, max);
419        test_exact(10.8, u64::MAX, max);
420    }
421
422    #[test]
423    fn test_pmf() {
424        let pmf = |arg: u64| move |x: Poisson| x.pmf(arg);
425        test_absolute(1.5, 0.334695240222645000000000000000, 1e-15, pmf(1));
426        test_absolute(1.5, 0.000003545747740570180000000000, 1e-20, pmf(10));
427        test_absolute(1.5, 0.000000000000000304971208961018, 1e-30, pmf(20));
428        test_absolute(5.4, 0.024389537090108400000000000000, 1e-17, pmf(1));
429        test_absolute(5.4, 0.026241240591792300000000000000, 1e-16, pmf(10));
430        test_absolute(5.4, 0.000000825202200316548000000000, 1e-20, pmf(20));
431        test_absolute(10.8, 0.000220314636840657000000000000, 1e-18, pmf(1));
432        test_absolute(10.8, 0.121365183659420000000000000000, 1e-15, pmf(10));
433        test_absolute(10.8, 0.003908139778574110000000000000, 1e-16, pmf(20));
434    }
435
436    #[test]
437    fn test_ln_pmf() {
438        let ln_pmf = |arg: u64| move |x: Poisson| x.ln_pmf(arg);
439        test_absolute(1.5, -1.09453489189183485135413967177, 1e-15, ln_pmf(1));
440        test_absolute(1.5, -12.5497614919938728510400000000, 1e-14, ln_pmf(10));
441        test_absolute(1.5, -35.7263142985901000000000000000, 1e-13, ln_pmf(20));
442        test_exact(5.4, -3.71360104642977159156055355910, ln_pmf(1));
443        test_absolute(5.4, -3.64042303737322774736223038530, 1e-15, ln_pmf(10));
444        test_absolute(5.4, -14.0076373893489089949388000000, 1e-14, ln_pmf(20));
445        test_absolute(10.8, -8.42045386586982559781714423000, 1e-14, ln_pmf(1));
446        test_absolute(10.8, -2.10895123177378079525424989992, 1e-14, ln_pmf(10));
447        test_absolute(10.8, -5.54469377815000936289610059500, 1e-14, ln_pmf(20));
448    }
449
450    #[test]
451    fn test_cdf() {
452        let cdf = |arg: u64| move |x: Poisson| x.cdf(arg);
453        test_absolute(1.5, 0.5578254003710750000000, 1e-15, cdf(1));
454        test_absolute(1.5, 0.9999994482467640000000, 1e-15, cdf(10));
455        test_exact(1.5, 1.0, cdf(20));
456        test_absolute(5.4, 0.0289061180327211000000, 1e-16, cdf(1));
457        test_absolute(5.4, 0.9774863006897650000000, 1e-15, cdf(10));
458        test_absolute(5.4, 0.9999997199928290000000, 1e-15, cdf(20));
459        test_absolute(10.8, 0.0002407141402518290000, 1e-16, cdf(1));
460        test_absolute(10.8, 0.4839692359955690000000, 1e-15, cdf(10));
461        test_absolute(10.8, 0.9961800769608090000000, 1e-15, cdf(20));
462    }
463
464    #[test]
465    fn test_sf() {
466        let sf = |arg: u64| move |x: Poisson| x.sf(arg);
467        test_absolute(1.5, 0.44217459962892536, 1e-15, sf(1));
468        test_absolute(1.5, 0.0000005517532358246565, 1e-15, sf(10));
469        test_absolute(1.5, 2.3372210700347092e-17, 1e-15, sf(20));
470        test_absolute(5.4, 0.971093881967279, 1e-16, sf(1));
471        test_absolute(5.4, 0.022513699310235582, 1e-15, sf(10));
472        test_absolute(5.4, 0.0000002800071708975261, 1e-15, sf(20));
473        test_absolute(10.8, 0.9997592858597482, 1e-16, sf(1));
474        test_absolute(10.8, 0.5160307640044303, 1e-15, sf(10));
475        test_absolute(10.8, 0.003819923039191422, 1e-15, sf(20));
476    }
477
478    #[test]
479    fn test_discrete() {
480        test::check_discrete_distribution(&create_ok(0.3), 10);
481        test::check_discrete_distribution(&create_ok(4.5), 30);
482    }
483}
statrs/distribution/poisson.rs

statrs/distribution/
poisson.rs
