vendor/rand-0.7.3/src/distributions/poisson.rs

   1 // Copyright 2018 Developers of the Rand project.
   2 // Copyright 2016-2017 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 Poisson distribution.
  11 #![allow(deprecated)]
  12
  13 use crate::distributions::utils::log_gamma;
  14 use crate::distributions::{Cauchy, Distribution};
  15 use crate::Rng;
  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`.
  21 #[deprecated(since = "0.7.0", note = "moved to rand_distr crate")]
  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
  93                 let check = 0.9
  94                     * (1.0 + comp_dev * comp_dev)
  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 {
 109     use super::Poisson;
 110     use crate::distributions::Distribution;
 111
 112     #[test]
 113     #[cfg_attr(miri, ignore)] // Miri is too slow
 114     fn test_poisson_10() {
 115         let poisson = Poisson::new(10.0);
 116         let mut rng = crate::test::rng(123);
 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);
 130         let mut rng = crate::test::rng(123);
 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 }