Skip to content

Commit

Permalink
rand_distr: Add support for f64 to Poisson
Browse files Browse the repository at this point in the history
Internally, we generate `f64`, so it makes sense to let the user access
that result directly, instead of forcing a conversion from `u64`.
  • Loading branch information
vks committed May 15, 2019
1 parent 31aad14 commit c03e2c8
Showing 1 changed file with 21 additions and 15 deletions.
36 changes: 21 additions & 15 deletions rand_distr/src/poisson.rs
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ use crate::utils::log_gamma;
/// use rand_distr::{Poisson, Distribution};
///
/// let poi = Poisson::new(2.0).unwrap();
/// let v = poi.sample(&mut rand::thread_rng());
/// let v: u64 = poi.sample(&mut rand::thread_rng());
/// println!("{} is from a Poisson(2) distribution", v);
/// ```
#[derive(Clone, Copy, Debug)]
Expand Down Expand Up @@ -62,30 +62,28 @@ impl Poisson {
}
}

impl Distribution<u64> for Poisson {
fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> u64 {
impl Distribution<f64> for Poisson {
fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> f64 {
// 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 result = 0.;
let mut p = 1.0;
while p > self.exp_lambda {
p *= rng.gen::<f64>();
result += 1;
result += 1.;
}
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).unwrap();
let mut result;

loop {
let mut result;
let mut comp_dev;

loop {
Expand All @@ -101,9 +99,6 @@ impl Distribution<u64> for Poisson {
// now the result is a random variable greater than 0 with Cauchy distribution
// the result should be an integer value
result = result.floor();
assert!(result >= 0.);
assert!(result <= ::core::u64::MAX as f64);
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
Expand All @@ -117,11 +112,20 @@ impl Distribution<u64> for Poisson {
break;
}
}
int_result
result
}
}
}

impl Distribution<u64> for Poisson {
fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> u64 {
let result: f64 = self.sample(rng);
assert!(result >= 0.);
assert!(result <= ::core::u64::MAX as f64);
result as u64
}
}

#[cfg(test)]
mod test {
use crate::Distribution;
Expand All @@ -133,7 +137,8 @@ mod test {
let mut rng = crate::test::rng(123);
let mut sum = 0;
for _ in 0..1000 {
sum += poisson.sample(&mut rng);
let s: u64 = poisson.sample(&mut rng);
sum += s;
}
let avg = (sum as f64) / 1000.0;
println!("Poisson average: {}", avg);
Expand All @@ -147,7 +152,8 @@ mod test {
let mut rng = crate::test::rng(123);
let mut sum = 0;
for _ in 0..1000 {
sum += poisson.sample(&mut rng);
let s: u64 = poisson.sample(&mut rng);
sum += s;
}
let avg = (sum as f64) / 1000.0;
println!("Poisson average: {}", avg);
Expand Down

0 comments on commit c03e2c8

Please sign in to comment.