diff options
author | Daniel Mueller <deso@posteo.net> | 2019-01-02 21:14:10 -0800 |
---|---|---|
committer | Daniel Mueller <deso@posteo.net> | 2019-01-02 21:14:10 -0800 |
commit | ecf3474223ca3d16a10f12dc2272e3b0ed72c1bb (patch) | |
tree | 03134a683791176b49ef5c92e8d6acd24c3b5a9b /rand/src/distributions/poisson.rs | |
parent | 686f61b75055ecb02baf9d9449525ae447a3bed1 (diff) | |
download | nitrocli-ecf3474223ca3d16a10f12dc2272e3b0ed72c1bb.tar.gz nitrocli-ecf3474223ca3d16a10f12dc2272e3b0ed72c1bb.tar.bz2 |
Update nitrokey crate to 0.2.3
This change updates the nitrokey crate to version 0.2.3. This version
bumps the rand crate used to 0.6.1, which in turn requires an additional
set of dependencies.
Import subrepo nitrokey/:nitrokey at b3e2adc5bb1300441ca74cc7672617c042f3ea31
Import subrepo rand/:rand at 73613ff903512e9503e41cc8ba9eae76269dc598
Import subrepo rustc_version/:rustc_version at 0294f2ba2018bf7be672abd53db351ce5055fa02
Import subrepo semver-parser/:semver-parser at 750da9b11a04125231b1fb293866ca036845acee
Import subrepo semver/:semver at 5eb6db94fa03f4d5c64a625a56188f496be47598
Diffstat (limited to 'rand/src/distributions/poisson.rs')
-rw-r--r-- | rand/src/distributions/poisson.rs | 157 |
1 files changed, 157 insertions, 0 deletions
diff --git a/rand/src/distributions/poisson.rs b/rand/src/distributions/poisson.rs new file mode 100644 index 0000000..1244caa --- /dev/null +++ b/rand/src/distributions/poisson.rs @@ -0,0 +1,157 @@ +// Copyright 2018 Developers of the Rand project. +// Copyright 2016-2017 The Rust Project Developers. +// +// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or +// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license +// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your +// option. This file may not be copied, modified, or distributed +// except according to those terms. + +//! The Poisson distribution. + +use Rng; +use distributions::{Distribution, Cauchy}; +use distributions::utils::log_gamma; + +/// The Poisson distribution `Poisson(lambda)`. +/// +/// This distribution has a density function: +/// `f(k) = lambda^k * exp(-lambda) / k!` for `k >= 0`. +/// +/// # Example +/// +/// ``` +/// use rand::distributions::{Poisson, Distribution}; +/// +/// let poi = Poisson::new(2.0); +/// let v = poi.sample(&mut rand::thread_rng()); +/// println!("{} is from a Poisson(2) distribution", v); +/// ``` +#[derive(Clone, Copy, Debug)] +pub struct Poisson { + lambda: f64, + // precalculated values + exp_lambda: f64, + log_lambda: f64, + sqrt_2lambda: f64, + magic_val: f64, +} + +impl Poisson { + /// Construct a new `Poisson` with the given shape parameter + /// `lambda`. Panics if `lambda <= 0`. + pub fn new(lambda: f64) -> Poisson { + assert!(lambda > 0.0, "Poisson::new called with lambda <= 0"); + let log_lambda = lambda.ln(); + Poisson { + lambda, + exp_lambda: (-lambda).exp(), + log_lambda, + sqrt_2lambda: (2.0 * lambda).sqrt(), + magic_val: lambda * log_lambda - log_gamma(1.0 + lambda), + } + } +} + +impl Distribution<u64> for Poisson { + fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> u64 { + // using the algorithm from Numerical Recipes in C + + // for low expected values use the Knuth method + if self.lambda < 12.0 { + let mut result = 0; + let mut p = 1.0; + while p > self.exp_lambda { + p *= rng.gen::<f64>(); + result += 1; + } + result - 1 + } + // high expected values - rejection method + else { + let mut int_result: u64; + + // we use the Cauchy distribution as the comparison distribution + // f(x) ~ 1/(1+x^2) + let cauchy = Cauchy::new(0.0, 1.0); + + loop { + let mut result; + let mut comp_dev; + + loop { + // draw from the Cauchy distribution + comp_dev = rng.sample(cauchy); + // shift the peak of the comparison ditribution + result = self.sqrt_2lambda * comp_dev + self.lambda; + // repeat the drawing until we are in the range of possible values + if result >= 0.0 { + break; + } + } + // now the result is a random variable greater than 0 with Cauchy distribution + // the result should be an integer value + result = result.floor(); + int_result = result as u64; + + // this is the ratio of the Poisson distribution to the comparison distribution + // the magic value scales the distribution function to a range of approximately 0-1 + // since it is not exact, we multiply the ratio by 0.9 to avoid ratios greater than 1 + // this doesn't change the resulting distribution, only increases the rate of failed drawings + let check = 0.9 * (1.0 + comp_dev * comp_dev) + * (result * self.log_lambda - log_gamma(1.0 + result) - self.magic_val).exp(); + + // check with uniform random value - if below the threshold, we are within the target distribution + if rng.gen::<f64>() <= check { + break; + } + } + int_result + } + } +} + +#[cfg(test)] +mod test { + use distributions::Distribution; + use super::Poisson; + + #[test] + fn test_poisson_10() { + let poisson = Poisson::new(10.0); + let mut rng = ::test::rng(123); + let mut sum = 0; + for _ in 0..1000 { + sum += poisson.sample(&mut rng); + } + let avg = (sum as f64) / 1000.0; + println!("Poisson average: {}", avg); + assert!((avg - 10.0).abs() < 0.5); // not 100% certain, but probable enough + } + + #[test] + fn test_poisson_15() { + // Take the 'high expected values' path + let poisson = Poisson::new(15.0); + let mut rng = ::test::rng(123); + let mut sum = 0; + for _ in 0..1000 { + sum += poisson.sample(&mut rng); + } + let avg = (sum as f64) / 1000.0; + println!("Poisson average: {}", avg); + assert!((avg - 15.0).abs() < 0.5); // not 100% certain, but probable enough + } + + #[test] + #[should_panic] + fn test_poisson_invalid_lambda_zero() { + Poisson::new(0.0); + } + + #[test] + #[should_panic] + fn test_poisson_invalid_lambda_neg() { + Poisson::new(-10.0); + } +} |