how to ask boost::normal_distribution to generate a large vector of random variables without for loop

Question

I have a large vector, I want to add noise with normal distribution to it. what I am doing now is trivial for loop:

for (int i=0 ; i<size ; i++){    //size is of order 1000000
    boost::mt19937 gen;
    gen.seed(i);
    boost::normal_distribution<>  nd(0.0 , 1.0);
    boost::variate_generator<boost::mt19937&, boost::normal_distribution<> >
                                                         randNormal (gen,nd);

    noise = randNormal();  
    nsyData[i] = data[i] + sd*noise; 
    }

Is there an efficient way to do this? like what MATLAB does?

Answer 1

Here's my take on it:

#include <boost/random/normal_distribution.hpp>
#include <boost/random.hpp>

int main()
{
    boost::mt19937 gen(42); // seed it once
    boost::normal_distribution<double> nd(0.0, 1.0);
    boost::variate_generator<boost::mt19937&, boost::normal_distribution<double> > randNormal(gen, nd);

    std::vector<double> data(100000, 0.0), nsyData;
    nsyData.reserve(data.size());

    double sd = 415*randNormal();
    std::transform(data.begin(), data.end(), std::back_inserter(nsyData), 
            [sd,&randNormal](double data) { return data + sd*randNormal(); });

}

Note that you were seeding the mersenne twister on each iteration of the loop. I'm afraid this totally killed any quality guarantees of the generated random numbers. Seed your generator once. (Use a different seed, e.g. from a random_device if you need it to be non-deterministic, obviously).

See this Live On Coliru

Update After some back and forth in the comments, here's a c++03 version that should not actually be too bad while still being comprehensible:

#include <boost/random/normal_distribution.hpp>
#include <boost/random.hpp>
#include <boost/bind.hpp>

struct Xfrm { 
    typedef double result_type;

    template <typename Gen>
    double operator()(double sd, Gen& randNormal, double data) const {
        return data + sd*randNormal();
    }
};

int main()
{
    boost::mt19937 gen(42); // seed it once
    boost::normal_distribution<double> nd(0.0, 1.0);
    boost::variate_generator<boost::mt19937&, boost::normal_distribution<double> > randNormal(gen, nd);

    std::vector<double> data(100000, 0.0), nsyData;
    nsyData.reserve(data.size());

    double sd = 415*randNormal();

    std::transform(data.begin(), data.end(), std::back_inserter(nsyData), 
            boost::bind(Xfrm(), sd, boost::ref(randNormal), ::_1));
}