mercredi 10 mars 2021

(Why) is the std::binomial_distribution biased for large probabilities p and slow for small n?

I want to generate binomially distributed random numbers in c++. Speed is a major concern. Not knowing a lot about random number generators, I use the standard libraries' tools. My code looks like something below:

#include <random>

static std::random_device random_dev;
static std::mt19937 random_generator{random_dev()};
std::binomial_distribution<int> binomial_generator;

void RandomInit(int s) {
    //I create the generator object here to save time. Does this make sense?
    binomial_generator = std::binomial_distribution<int>(1, 0.5);
    random_generator.seed(s);
}

int binomrand(int n, double p) {
    binomial_generator.param(std::binomial_distribution<int>::param_type(n, p));
    return binomial_generator(random_generator);
}

To test my implementation, I have built a cython wrapper and then executed and timed the function from within python. For reference I have also implemented a "stupid" binomial distribution, which jest returns the sum of Bernoulli trials.

int binomrand2(int n, double p) {
    int result = 0;
    for (int i = 0; i<n; i++) {
        if (_Random() < p) //_Random is a thoroughly tested custom random number generator on U[0,1)
            result++;
    }

    return result;
}

Timing showed that the latter implementation is about 50% faster than the former if n < 25. Furthermore, for p = 0.95, the former yielded significantly biased results (the mean over 1000000 trials for n = 40 was 38.23037; standard deviation is 0.0014; the result was reproducable with different seeds).

Is this a (known) issue with the standard library's functions or is my implementation wrong? What could I do to achieve my goal of obtaining accurate results with high efficiency?

I am open to suggestions outside the realm of the standard library, but I may not be able to use external software libraries.




Aucun commentaire:

Enregistrer un commentaire