statrs/distribution/
poisson.rs

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

statrs/distribution/
poisson.rs
