Random number reproducibility from Poisson distribution

General Usage
, ,
1

I want to sample a series of random numbers from Poisson distributions, such as

using Distributions, Random

Random.seed!(1)
rand(Poisson(10))  # 9
rand(Poisson(11))  # 12
rand(Poisson(12))  # 9
rand(Poisson(13))  # 17

I would then like to be able to change some values in this sequence, but have unchanged lines return the same result as before:

Random.seed!(1)
rand(Poisson(10))  # 9
rand(Poisson(20))  # 22
rand(Poisson(12))  # 9
rand(Poisson(13))  # 17

However, if any values >=6 switch to <6, or vice versa, then this disrupts the seed for subsequent values:

Random.seed!(1)

rand(Poisson(10))  # 9
rand(Poisson(5))  # 2
rand(Poisson(12))  # 16
rand(Poisson(13))  # 13

I have traced the source of this difference to line 138 of Distributions.jl/src/univariate/discrete/poisson.jl at master · JuliaStats/Distributions.jl · GitHub, where different samplers are used on different sides of poissonsampler_threshold = 6.

Is there a way to have the unchanged lines return the same value regardless of whether other lines switch around this threshold?

2

Your problem is not that the algorithm changes at some threshold. Your problem is instead that the number of random bits extracted from the RNG is non-uniform and dependent on both RNG state and parameters.

This also holds for the sampling algorithm beyond the threshold.

It is indeed hard to avoid that issue: If you consume at most k bits from the random stream, then all probabilities are integer multiples of 2^-k.

And random bits are cheap – none of the sampling algorithms are especially stingy in their random consumption.

The appropriate generic way is to fork the RNG state, call whatever sampling algorithm with the forked RNG, and then discard the fork.

3 Likes
3

Thank you. I’ll have a go at that.