Generate random samples from multivariate distribution

[dr86]

Hello.


At the link:

https://stackoverflow.com/questions/6142576/sample-from-multivariate-normal-gaussian-distribution-in-c,

a program to sample values from multivariate normal Gaussian distribution can be found. However, every time when I start program, same values are generated (comparing to previous run of program). I need randomness such that different values are generated every time when I enter program. Here is program:

#include <Eigen/Dense>
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/normal_distribution.hpp>    

/*
  We need a functor that can pretend it's const,
  but to be a good random number generator 
  it needs mutable state.
*/
namespace Eigen {
namespace internal {
template<typename Scalar> 
struct scalar_normal_dist_op 
{
  static boost::mt19937 rng;    // The uniform pseudo-random algorithm
  mutable boost::normal_distribution<Scalar> norm;  // The gaussian combinator

  EIGEN_EMPTY_STRUCT_CTOR(scalar_normal_dist_op)

  template<typename Index>
  inline const Scalar operator() (Index, Index = 0) const { return norm(rng); }
};

template<typename Scalar> boost::mt19937 scalar_normal_dist_op<Scalar>::rng;

template<typename Scalar>
struct functor_traits<scalar_normal_dist_op<Scalar> >
{ enum { Cost = 50 * NumTraits<Scalar>::MulCost, PacketAccess = false, IsRepeatable = false }; };
} // end namespace internal
} // end namespace Eigen

/*
  Draw nn samples from a size-dimensional normal distribution
  with a specified mean and covariance
*/
void main() 
{
  int size = 2; // Dimensionality (rows)
  int nn=5;     // How many samples (columns) to draw
  Eigen::internal::scalar_normal_dist_op<double> randN; // Gaussian functor
  Eigen::internal::scalar_normal_dist_op<double>::rng.seed(1); // Seed the rng

  // Define mean and covariance of the distribution
  Eigen::VectorXd mean(size);       
  Eigen::MatrixXd covar(size,size);

  mean  <<  0,  0;
  covar <<  1, .5,
           .5,  1;

  Eigen::MatrixXd normTransform(size,size);

  Eigen::LLT<Eigen::MatrixXd> cholSolver(covar);

  // We can only use the cholesky decomposition if 
  // the covariance matrix is symmetric, pos-definite.
  // But a covariance matrix might be pos-semi-definite.
  // In that case, we'll go to an EigenSolver
  if (cholSolver.info()==Eigen::Success) {
    // Use cholesky solver
    normTransform = cholSolver.matrixL();
  } else {
    // Use eigen solver
    Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigenSolver(covar);
    normTransform = eigenSolver.eigenvectors() 
                   * eigenSolver.eigenvalues().cwiseSqrt().asDiagonal();
  }

  Eigen::MatrixXd samples = (normTransform 
                           * Eigen::MatrixXd::NullaryExpr(size,nn,randN)).colwise() 
                           + mean;

  std::cout << "Mean\n" << mean << std::endl;
  std::cout << "Covar\n" << covar << std::endl;
  std::cout << "Samples\n" << samples << std::endl;
}

I tried (instead of static boost::mt19937 rng;):

static boost::random::mt19937 rng{static_cast<std::uint32_t>(now)};

and there is error:

error: #error This file requires compiler and library support for the ISO C++ 2011 standard. I will not change compiler, so this option is off.

I also tried to use boost::variate_generator,

struct scalar_normal_dist_op 
{
   static boost::mt19937 rng;    // The uniform pseudo-random algorithm
   mutable boost::normal_distribution<Scalar> norm;  // The gaussian combinator
   boost::variate_generator<boost::mt19937, boost::normal_distribution<> > var;
  
  
  EIGEN_EMPTY_STRUCT_CTOR(scalar_normal_dist_op)

  template<typename Index>
  inline const Scalar operator() (Index, Index = 0) const { return var(rng, norm);}
};

with various combinations of commands, without success. Any advice is appreciated.

Regards,
Darko

[Peter87]

std::uint32_t is C++11. You can probably use boost::uint32_t instead.
Brace initialization {} were very limited before C++11. In this case you need to use regular parentheses instead.

static boost::random::mt19937 rng(static_cast<boost::uint32_t>(now));

[dr86]

Peter,

Thank you for response. When I use command that you recommended, when I compile program, there are still errors (this is the beginning):

dare.cpp:15:38: error: expected identifier before ‘static_cast’
static boost::random::mt19937 rng(static_cast<boost::uint32_t>(now));
^
dare.cpp:15:38: error: expected ‘,’ or ‘...’ before ‘static_cast’
dare.cpp:26:73: error: ‘boost::random::mt19937 Eigen::internal::scalar_normal_dist_op<Scalar>::rng’ is not a static data member of ‘struct Eigen::internal::scalar_normal_dist_op<Scalar>’
template<typename Scalar> boost::mt19937 scalar_normal_dist_op<Scalar>::rng;
^
dare.cpp:26:73: error: template definition of non-template ‘boost::random::mt19937 Eigen::internal::scalar_normal_dist_op<Scalar>::rng’
dare.cpp: In function ‘int main()’:
dare.cpp:43:55: error: request for member ‘seed’ in ‘Eigen::internal::scalar_normal_dist_op<Scalar>::rng<double>’, which is of non-class type ‘boost::random::mt19937(int) {aka boost::random::mersenne_twister_engine<unsigned int, 32ul, 624ul, 397ul, 31ul, 2567483615u, 11ul, 4294967295u, 7ul, 2636928640u, 15ul, 4022730752u, 18ul, 1812433253u>(int)}’
Eigen::internal::scalar_normal_dist_op<double>::rng.seed(1); // Seed the rng

^

[Peter87]

Ah, it's a class template. I obviously wasn't looking at the code closely enough.

Well, in that case the constructor arguments need to be passed where you actually define the static variable, outside the class body.

Something like this:

template<typename Scalar>
boost::mt19937 scalar_normal_dist_op<Scalar>::rng(static_cast<boost::uint32_t>(now));

But now the now variable looks a bit out of place. You might want to pass std::time(nullptr), or whatever you use as seed, directly here.

But hold on a minute. You are actually seeding the random generator already. On line 41. The problem is that you always use 1 as the seed. You might only want to change this line. If you decide to specify the seed when constructing the random engine, as we tried to do earlier, you will have to remove line 41 otherwise it will override the seed that you passed to the constructor.

[dr86]

Peter,

Thank you very much for your help.

Line 41 is modified:

Eigen::internal::scalar_normal_dist_op<double>::rng.seed(time(0)); // Seed the rng


Regards,
Darko