Before you use a random number generation library, it’s a good idea to do a few tests. The tests might uncover a flaw in the library, but more likely they’ll uncover a flaw in your understanding of the library. Parameterization expectations in particular are a common source of errors.
The following code does a crude test of the C++ TR1 random number generation classes. It computes the mean of 10,000 samples from each of the supported distributions and does a normal approximation to see whether the sample mean is roughly what one would expect. Don’t be surprised if values fall outside the expected interval. But if you run the code several times with different RNG seeds, you shouldn’t expect the same sample to go outside its predicted range very often.
This code is hardly a rigorous test of the distributions. But it provides sample code for working with the distributions and tests your understanding of the parameterizations.
#include <random>
#include <iostream>
#include <string>
template <typename TDistribution>
double sample_mean(TDistribution dist, int sample_size = 10000)
{
std::tr1::mt19937 mt; // Mersenne Twister generator
double sum = 0;
for (int i = 0; i < sample_size; ++i)
sum += dist(mt);
return sum / sample_size;
}
template <typename TDistribution>
void test_mean(TDistribution dist, std::string name, double true_mean, double true_variance)
{
double true_stdev = sqrt(true_variance);
std::cout << "Testing " << name << " distribution\n";
int sample_size = 10000;
double mean = sample_mean(dist, sample_size);
std::cout << "Computed mean: " << mean << "\n";
double lower = true_mean - true_stdev/sqrt((double)sample_size);
double upper = true_mean + true_stdev/sqrt((double)sample_size);
std::cout << "Expected a value between " << lower << " and " << upper << "\n\n";
}
int main()
{
int n;
double p, lambda, shape, mu, sigma;
n = 5; p = 0.3;
std::tr1::binomial_distribution<int, double> binomial(n, p);
test_mean(binomial, "binomial", n*p, n*p*(1-p));
lambda = 4.0;
std::tr1::exponential_distribution<double> exponential(lambda);
test_mean(exponential, "exponential", 1.0/lambda, 1.0/(lambda*lambda));
shape = 3.0;
std::tr1::gamma_distribution<double> gamma(shape);
test_mean(gamma, "gamma", shape, shape);
p = 0.1;
std::tr1::geometric_distribution<int, double> geometric(p);
test_mean(geometric, "geometric", 1.0/p, (1.0 - p)/(p*p));
mu = 3.0; sigma = 4.0;
std::tr1::normal_distribution<double> normal(mu, sigma);
test_mean(normal, "normal", mu, sigma*sigma);
lambda = 7.0;
std::tr1::poisson_distribution<double> poisson(7.0);
test_mean(poisson, "poisson", lambda, lambda);
p = 0.6;
std::tr1::bernoulli_distribution bernoulli(p);
test_mean(bernoulli, "bernoulli", p, p*(1-p));
return (0);
}