C++ TR1: how to use the normal_distribution?

Question

I'm trying to use the C++ STD TechnicalReport1 extensions to generate numbers following a normal distribution, but this code (adapted from this article):

mt19937 eng;
eng.seed(SEED);

normal_distribution<double> dist;
// XXX if I use the one below it exits the for loop
// uniform_int<int> dist(1, 52);

for (unsigned int i = 0; i < 1000; ++i) {
  cout << "Generating " << i << "-th value" << endl;
  cout << dist(eng) << endl;
}

only prints 1 "Generating..." log message, then never exits the for loop! If I use the distribution I commented out instead, it terminates, so I'm wondering what I'm doing wrong. Any idea?

Thanks a lot!

Undertone

Answer 1

This definitely would not hang the program. But, not sure if it really meets your needs.

 #include <random>
 #include <iostream>

 using namespace std;

 typedef std::tr1::ranlux64_base_01 Myeng; 

 typedef std::tr1::normal_distribution<double> Mydist; 

 int main() 
 { 
      Myeng eng; 
      eng.seed(1000);
      Mydist dist(1,10); 

      dist.reset(); // discard any cached values 
      for (int i = 0; i < 10; i++)
      {
           std::cout << "a random value == " << (int)dist(eng) << std::endl; 
      }

 return (0); 
 }

Schnorr

thanks man, it works like a charm, but I'm wondering why with this engine it works, and not with the other.. – Undertone

Obviously the only difference is your using the mt19937 number generator whereas Jagannath uses the std::tr1::ranlux64_base_01. Logically, I guess the bug may be in your implementation of the mt19937 object ( algo which I had never heard about before you did, thx for this :-) ) that is not part of std library. – Unrelenting

Is it possible to vectorize such a for loop when drawing random numbers? I recall that you cannot vectorize a loop which has a function call. – Frenchman

Answer 2

I have had the same issue with the code originally posted and investigated the GNU implementation of

first some observations: with g++-4.4 and using the code hangs, with g++-4.5 and using -std=c++0x (i.e. not TR1 but the real thing) above code works

IMHO, there was a change between TR1 and c++0x with regard to adaptors between random number generation and consumption of random numbers -- mt19937 produces integers, normal_distribution consumes doubles

the c++0x uses adaption automatically, the g++ TR1 code does not

in order to get your code working with g++-4.4 and TR1, do the following

std::tr1::mt19937 prng(seed);
std::tr1::normal_distribution<double> normal;
std::tr1::variate_generator<std::tr1::mt19937, std::tr1::normal_distribution<double> > randn(prng,normal);
double r = randn();

Omeara

Answer 3

If your TR1 random number generation implementation is buggy, you can avoid TR1 by writing your own normal generator as follows.

Generate two uniform (0, 1) random samples u and v using any random generator you trust. Then let r = sqrt( -2 log(u) ) and return x = r sin(2 pi v). (This is called the Box-Mueller method.)

If you need normal samples samples with mean mu and standard deviation sigma, return sigma*x + mu instead of just x.

Irrelative

Just tried this. It runs super fast -- tested with 1M samples and it presented almost perfect stats for samples within 1-sigma, 2-sigma, etc. – Haslet

Answer 4

While this appears to be a bug, a quick confirmation would be to pass the default 0.0, 1.0 parameters. normal_distribution<double>::normal_distribution() should equal normal_distribution<double>::normal_distribution(0.0, 1.0)

Colporteur

it does not work either, it still remains stuck performing the first computation.. – Undertone