]>
Commit | Line | Data |
---|---|---|
0731742a XL |
1 | // Copyright 2018 Developers of the Rand project. |
2 | // Copyright 2016-2017 The Rust Project Developers. | |
b7449926 XL |
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 Poisson distribution. | |
416331ca | 11 | #![allow(deprecated)] |
b7449926 | 12 | |
416331ca | 13 | use crate::distributions::utils::log_gamma; |
dfeec247 XL |
14 | use crate::distributions::{Cauchy, Distribution}; |
15 | use crate::Rng; | |
b7449926 XL |
16 | |
17 | /// The Poisson distribution `Poisson(lambda)`. | |
18 | /// | |
19 | /// This distribution has a density function: | |
20 | /// `f(k) = lambda^k * exp(-lambda) / k!` for `k >= 0`. | |
dfeec247 | 21 | #[deprecated(since = "0.7.0", note = "moved to rand_distr crate")] |
b7449926 XL |
22 | #[derive(Clone, Copy, Debug)] |
23 | pub struct Poisson { | |
24 | lambda: f64, | |
25 | // precalculated values | |
26 | exp_lambda: f64, | |
27 | log_lambda: f64, | |
28 | sqrt_2lambda: f64, | |
29 | magic_val: f64, | |
30 | } | |
31 | ||
32 | impl Poisson { | |
33 | /// Construct a new `Poisson` with the given shape parameter | |
34 | /// `lambda`. Panics if `lambda <= 0`. | |
35 | pub fn new(lambda: f64) -> Poisson { | |
36 | assert!(lambda > 0.0, "Poisson::new called with lambda <= 0"); | |
37 | let log_lambda = lambda.ln(); | |
38 | Poisson { | |
39 | lambda, | |
40 | exp_lambda: (-lambda).exp(), | |
41 | log_lambda, | |
42 | sqrt_2lambda: (2.0 * lambda).sqrt(), | |
43 | magic_val: lambda * log_lambda - log_gamma(1.0 + lambda), | |
44 | } | |
45 | } | |
46 | } | |
47 | ||
48 | impl Distribution<u64> for Poisson { | |
49 | fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> u64 { | |
50 | // using the algorithm from Numerical Recipes in C | |
51 | ||
52 | // for low expected values use the Knuth method | |
53 | if self.lambda < 12.0 { | |
54 | let mut result = 0; | |
55 | let mut p = 1.0; | |
56 | while p > self.exp_lambda { | |
57 | p *= rng.gen::<f64>(); | |
58 | result += 1; | |
59 | } | |
60 | result - 1 | |
61 | } | |
62 | // high expected values - rejection method | |
63 | else { | |
64 | let mut int_result: u64; | |
65 | ||
66 | // we use the Cauchy distribution as the comparison distribution | |
67 | // f(x) ~ 1/(1+x^2) | |
68 | let cauchy = Cauchy::new(0.0, 1.0); | |
69 | ||
70 | loop { | |
71 | let mut result; | |
72 | let mut comp_dev; | |
73 | ||
74 | loop { | |
75 | // draw from the Cauchy distribution | |
76 | comp_dev = rng.sample(cauchy); | |
77 | // shift the peak of the comparison ditribution | |
78 | result = self.sqrt_2lambda * comp_dev + self.lambda; | |
79 | // repeat the drawing until we are in the range of possible values | |
80 | if result >= 0.0 { | |
81 | break; | |
82 | } | |
83 | } | |
84 | // now the result is a random variable greater than 0 with Cauchy distribution | |
85 | // the result should be an integer value | |
86 | result = result.floor(); | |
87 | int_result = result as u64; | |
88 | ||
89 | // this is the ratio of the Poisson distribution to the comparison distribution | |
90 | // the magic value scales the distribution function to a range of approximately 0-1 | |
91 | // since it is not exact, we multiply the ratio by 0.9 to avoid ratios greater than 1 | |
92 | // this doesn't change the resulting distribution, only increases the rate of failed drawings | |
dfeec247 XL |
93 | let check = 0.9 |
94 | * (1.0 + comp_dev * comp_dev) | |
b7449926 XL |
95 | * (result * self.log_lambda - log_gamma(1.0 + result) - self.magic_val).exp(); |
96 | ||
97 | // check with uniform random value - if below the threshold, we are within the target distribution | |
98 | if rng.gen::<f64>() <= check { | |
99 | break; | |
100 | } | |
101 | } | |
102 | int_result | |
103 | } | |
104 | } | |
105 | } | |
106 | ||
107 | #[cfg(test)] | |
108 | mod test { | |
b7449926 | 109 | use super::Poisson; |
dfeec247 | 110 | use crate::distributions::Distribution; |
b7449926 XL |
111 | |
112 | #[test] | |
dfeec247 | 113 | #[cfg_attr(miri, ignore)] // Miri is too slow |
b7449926 XL |
114 | fn test_poisson_10() { |
115 | let poisson = Poisson::new(10.0); | |
416331ca | 116 | let mut rng = crate::test::rng(123); |
b7449926 XL |
117 | let mut sum = 0; |
118 | for _ in 0..1000 { | |
119 | sum += poisson.sample(&mut rng); | |
120 | } | |
121 | let avg = (sum as f64) / 1000.0; | |
122 | println!("Poisson average: {}", avg); | |
123 | assert!((avg - 10.0).abs() < 0.5); // not 100% certain, but probable enough | |
124 | } | |
125 | ||
126 | #[test] | |
127 | fn test_poisson_15() { | |
128 | // Take the 'high expected values' path | |
129 | let poisson = Poisson::new(15.0); | |
416331ca | 130 | let mut rng = crate::test::rng(123); |
b7449926 XL |
131 | let mut sum = 0; |
132 | for _ in 0..1000 { | |
133 | sum += poisson.sample(&mut rng); | |
134 | } | |
135 | let avg = (sum as f64) / 1000.0; | |
136 | println!("Poisson average: {}", avg); | |
137 | assert!((avg - 15.0).abs() < 0.5); // not 100% certain, but probable enough | |
138 | } | |
139 | ||
140 | #[test] | |
141 | #[should_panic] | |
142 | fn test_poisson_invalid_lambda_zero() { | |
143 | Poisson::new(0.0); | |
144 | } | |
145 | ||
146 | #[test] | |
147 | #[should_panic] | |
148 | fn test_poisson_invalid_lambda_neg() { | |
149 | Poisson::new(-10.0); | |
150 | } | |
151 | } |