Skip to content

Commit

Permalink
Poisson u64 sampling (#1498)
Browse files Browse the repository at this point in the history
This addresses #1497 by adding
`Distribution<u64>`
It also solves #1312 by not
allowing `lambda` bigger than `1.844e19` (this also makes them always
fit into `u64`)
  • Loading branch information
benjamin-lieser authored Oct 1, 2024
1 parent f263820 commit e2092e9
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 15 deletions.
2 changes: 2 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ You may also find the [Upgrade Guide](https://rust-random.github.io/book/update.
- Add `UniformUsize` and use to make `Uniform` for `usize` portable (#1487)
- Remove support for generating `isize` and `usize` values with `Standard`, `Uniform` and `Fill` and usage as a `WeightedAliasIndex` weight (#1487)
- Require `Clone` and `AsRef` bound for `SeedableRng::Seed`. (#1491)
- Implement `Distribution<u64>` for `Poisson<f64>` (#1498)
- Limit the maximal acceptable lambda for `Poisson` to solve (#1312) (#1498)

## [0.9.0-alpha.1] - 2024-03-18
- Add the `Slice::num_choices` method to the Slice distribution (#1402)
Expand Down
44 changes: 29 additions & 15 deletions rand_distr/src/poisson.rs
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,6 @@ use rand::Rng;
/// This distribution has density function:
/// `f(k) = λ^k * exp(-λ) / k!` for `k >= 0`.
///
/// # Known issues
///
/// See documentation of [`Poisson::new`].
///
/// # Plot
///
/// The following plot shows the Poisson distribution with various values of `λ`.
Expand All @@ -40,7 +36,7 @@ use rand::Rng;
/// use rand_distr::{Poisson, Distribution};
///
/// let poi = Poisson::new(2.0).unwrap();
/// let v = poi.sample(&mut rand::thread_rng());
/// let v: f64 = poi.sample(&mut rand::thread_rng());
/// println!("{} is from a Poisson(2) distribution", v);
/// ```
#[derive(Clone, Copy, Debug, PartialEq)]
Expand All @@ -52,20 +48,23 @@ where

/// Error type returned from [`Poisson::new`].
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
// Marked non_exhaustive to allow a new error code in the solution to #1312.
#[non_exhaustive]
pub enum Error {
/// `lambda <= 0`
ShapeTooSmall,
/// `lambda = ∞` or `lambda = nan`
NonFinite,
/// `lambda` is too large, see [Poisson::MAX_LAMBDA]
ShapeTooLarge,
}

impl fmt::Display for Error {
fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
f.write_str(match self {
Error::ShapeTooSmall => "lambda is not positive in Poisson distribution",
Error::NonFinite => "lambda is infinite or nan in Poisson distribution",
Error::ShapeTooLarge => {
"lambda is too large in Poisson distribution, see Poisson::MAX_LAMBDA"
}
})
}
}
Expand Down Expand Up @@ -125,14 +124,7 @@ where
/// Construct a new `Poisson` with the given shape parameter
/// `lambda`.
///
/// # Known issues
///
/// Although this method should return an [`Error`] on invalid parameters,
/// some (extreme) values of `lambda` are known to return a [`Poisson`]
/// object which hangs when [sampled](Distribution::sample).
/// Large (less extreme) values of `lambda` may result in successful
/// sampling but with reduced precision.
/// See [#1312](https://github.com/rust-random/rand/issues/1312).
/// The maximum allowed lambda is [MAX_LAMBDA](Self::MAX_LAMBDA).
pub fn new(lambda: F) -> Result<Poisson<F>, Error> {
if !lambda.is_finite() {
return Err(Error::NonFinite);
Expand All @@ -145,11 +137,25 @@ where
let method = if lambda < F::from(12.0).unwrap() {
Method::Knuth(KnuthMethod::new(lambda))
} else {
if lambda > F::from(Self::MAX_LAMBDA).unwrap() {
return Err(Error::ShapeTooLarge);
}
Method::Rejection(RejectionMethod::new(lambda))
};

Ok(Poisson(method))
}

/// The maximum supported value of `lambda`
///
/// This value was selected such that
/// `MAX_LAMBDA + 1e6 * sqrt(MAX_LAMBDA) < 2^64 - 1`,
/// thus ensuring that the probability of sampling a value larger than
/// `u64::MAX` is less than 1e-1000.
///
/// Applying this limit also solves
/// [#1312](https://github.com/rust-random/rand/issues/1312).
pub const MAX_LAMBDA: f64 = 1.844e19;
}

impl<F> Distribution<F> for KnuthMethod<F>
Expand Down Expand Up @@ -232,6 +238,14 @@ where
}
}

impl Distribution<u64> for Poisson<f64> {
#[inline]
fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> u64 {
// `as` from float to int saturates
<Poisson<f64> as Distribution<f64>>::sample(self, rng) as u64
}
}

#[cfg(test)]
mod test {
use super::*;
Expand Down

0 comments on commit e2092e9

Please sign in to comment.