How do I plug a boost::random QRNG into a C++ standard distribution?

80 Views Asked by At

I'm trying to use a boost::random::sobol QRNG to sample one of the standard C++ distributions listed here https://en.cppreference.com/w/cpp/numeric/random, but the results are not correct.

Suppose I wish to sample a binomial distribution. According to https://en.cppreference.com/w/cpp/numeric/random/binomial_distribution, I can do this:

#include <iostream>
#include <map>
#include <random>

int main()
{
  std::random_device rd;
  std::mt19937 gen(rd());  // <-- PRNG

  // perform 4 trials, each succeeds 1 in 2 times
  std::binomial_distribution<> d(4, 0.5);

  std::map<int, int> hist;
  for (int n = 0; n != 10000; ++n)
    ++hist[d(gen)];

  for (auto const& [x, y] : hist)
    std::cout << x << ' ' << std::string(y / 100, '*') << '\n';
}

and that works well, producing output something like this:

0 ******
1 *************************
2 *************************************
3 ************************
4 *****

However, I want to use a quasi-random number generator (in particular, a boost::random::sobol). If I've understood the documentation and the hints in Running std::normal_distribution with user-defined random generator (unfortunately, that page includes a dead link to what looks like a code example :/) correctly (and clearly I haven't), I can adapt my code above, and do this:

#include <iostream>
#include <map>
#include <random>
#include <boost/random/sobol.hpp>

int main()
{
  boost::random::sobol gen(1);  // QRNG: 1D sobol sequence
  gen.seed(237);  // just set the seed to something

  // perform 4 trials, each succeeds 1 in 2 times
  std::binomial_distribution<> d(4, 0.5);

  std::map<int, int> hist;
  for (int n = 0; n != 10000; ++n)
    ++hist[d(gen)];  // <-- plug the QRNG in here

  for (auto const& [x, y] : hist)
    std::cout << x << ' ' << std::string(y / 100, '*') << '\n';
}

This code runs, which (I think) means that the boost::random::sobol satisfies the requirements of the uniform_random_bit_generator concept, but the results this time are:

2 ****************************************************************************************************

which is clearly wrong (I've tried a bunch of standard distributions and they all look similarly awful). So I've obviously misunderstood something. Any advice very much appreciated.

If it makes any difference, I'm working in Windows with the latest version of C++ (I think it's the experimental C++23) in Visual Studio, and use boost 1.80.

1

There are 1 best solutions below

0
T.C. On

This is going to be dependent on the library implementation of the distribution, and in general you cannot expect this to work - whether with boost or std. By definition, a sequence of numbers obtained from a QRNG is not random or uniform - it is required to have low discrepancy.

For example, if I follow the definition of the binomial distribution literally and perform 4 trials using 4 numbers drawn from the distribution like so:

  const auto cutoff = gen.max() / 2;

  std::map<int, int> hist;
  for (int n = 0; n != 10000; ++n) {
    int passes = 0;
    for (int i = 0; i < 4; ++i) {
        passes += (gen() <= cutoff);
    }
    ++hist[passes];
  }

Then we get the output you report - and it's easy to see why: the "low-discrepancy" part means that half of the four numbers drawn are likely in the first half of the range (so the trial succeeds) and the remaining half of them are in the second half. (This, incidentally, is how MSVC implements binomial_distribution for small numbers of trials.)

Now, this works if the distribution consumes one number from the engine per number generated (typical for inverse transform sampling); this is the case for boost::binomial_distribution if the number of expected successes is small, but is generally not the case for arbitrary distributions.