In order to draw random number from a Poisson distribution in C++, it is generally advised to use
RNG_type rng;
std::poisson_distribution<size_t> d(1e-6);
auto r = d(rng);
At each call of the std::poisson_distribution object, an entire sequence of random bits is consumed (e.g. 32 bits with std::mt19937, 64 bits for std::mt19937_64). It strikes me that with such low mean (mean = 1e-6), the vast majority of times, only a few bits are enough to determine that the value to return is 0. The other bits could then be cached for later use.
Assuming that a sequence of bits set to true is associated to a high returned value from the Poisson distribution, when using a mean of 1e-6, any sequence not starting with 19 trues necessarily returns a zero! Indeed,
1 - 1/2^19 < P(0, 1e-6) < 1 - 1/2^20
, where P(n, r) denotes the probability of drawing n from a Poisson distribution with mean r. An algorithm that does not waste bits would use one bit half of the time, two bits a quarter of the times, three bits an eighth of the times, ....
Is there an algorithm out there that can improve performance by consuming as few bits as possible when drawing Poisson numbers? Is there another way to improve performance compared to std::poisson_distribution when we consider a low mean?
In response to @Jarod42's comment who said
Wonder if using fewer bits don't break equiprobability...
I don't think it would break equiprobability. In a vague attempt to test it, I consider the same question with a simple bernoulli distribution. I am sampling true with a probability 1/2^4 and sampling false with a probability 1 - 1/2^4. The function drawWithoutWastingBits stops as soon as it sees a true in the cache and the function drawWastingBits consumes 4 bits whatever these bits are.
#include <iostream>
#include <vector>
#include <string>
#include <algorithm>
#include <random>
bool drawWithoutWastingBits(std::vector<bool>& cache, size_t& cache_index)
{
/*
Get a true with probability 1/2^4 (=1/16=0.0625) and a false otherwise
*/
size_t nbTrues = 0;
while (cache[cache_index])
{
++nbTrues;
++cache_index;
if (nbTrues == 4)
{
return true;
}
}
++cache_index;
return false;
}
bool drawWastingBits(std::vector<bool>& cache, size_t& cache_index)
{
/*
Get a true with probability 1/2^4 (=1/16=0.0625) and a false otherwise
*/
bool isAnyTrue = false;
for (size_t i = 0 ; i < 4; ++i)
{
if (cache[cache_index])
{
isAnyTrue = true;
}
++cache_index;
}
return !isAnyTrue;
}
int main()
{
/*
Just cache a lot of bits in advance in `cache`. The same sequence of bits will be used by both function.
I am just caching way enough bits to make sure they don't run out of bits below
I made sure to have the same number of zeros and ones so that any deviation is caused by the methodology and not by the RNG
*/
// Produce cache
std::vector<bool> cache;
size_t nbBitsToCache = 1e7;
cache.reserve(nbBitsToCache);
for (size_t i = 0 ; i < nbBitsToCache/2 ; ++i)
{
cache.push_back(false);
cache.push_back(true);
}
// Shuffle cache
{
std::mt19937 mt(std::random_device{}());
std::shuffle(cache.begin(), cache.end(), mt);
}
// Draw without wasting bits
{
size_t nbDraws = 1e6;
size_t cache_index = 0;
std::pair<size_t, size_t> outcomes = {0,0};
for (size_t r = 0 ; r < nbDraws ; ++r)
{
drawWithoutWastingBits(cache, cache_index) ? ++outcomes.first : ++outcomes.second;
assert(cache_index <= cache.size());
}
assert(outcomes.first + outcomes.second == nbDraws);
std::cout << "Draw Without Wasting Bits: prob true = " << (double)outcomes.first / nbDraws << "\n";
}
// Draw wasting bits
{
size_t nbDraws = 1e6;
size_t cache_index = 0;
std::pair<size_t, size_t> outcomes = {0,0};
for (size_t r = 0 ; r < nbDraws ; ++r)
{
drawWastingBits(cache, cache_index) ? ++outcomes.first : ++outcomes.second;
assert(cache_index <= cache.size());
}
assert(outcomes.first + outcomes.second == nbDraws);
std::cout << "Draw Wit Wasting Bits: prob true = " << (double)outcomes.first / nbDraws << "\n";
}
}
Possible output
Draw Without Wasting Bits: prob true = 0.062832
Draw Wit Wasting Bits: prob true = 0.062363
