来自C ++中多变量正态/高斯分布的样本 [英] Sample from multivariate normal/Gaussian distribution in C++

查看:709
本文介绍了来自C ++中多变量正态/高斯分布的样本的处理方法,对大家解决问题具有一定的参考价值,需要的朋友们下面随着小编来一起学习吧!

问题描述

我一直在寻找一种方便的方法来从多元正态分布中抽样。有人知道一个容易获得的代码片段吗?对于矩阵/向量,我更喜欢使用 Boost Eigen 或另一个我不熟悉的现象库,但我可以使用 GSL 。如果方法接受非负的无限协方差矩阵,而不是需要正定(例如,与Cholesky分解一样),我也喜欢它。这存在于MATLAB,NumPy等,但我很难找到一个现成的C / C ++解决方案。

I've been hunting for a convenient way to sample from a multivariate normal distribution. Does anyone know of a readily available code snippet to do that? For matrices/vectors, I'd prefer to use Boost or Eigen or another phenomenal library I'm not familiar with, but I could use GSL in a pinch. I'd also like it if the method accepted nonnegative-definite covariance matrices rather than requiring positive-definite (e.g., as with the Cholesky decomposition). This exists in MATLAB, NumPy, and others, but I've had a hard time finding a ready-made C/C++ solution.

如果我自己实现它,我会抱怨,但是很好。如果我这样做,维基百科使之听起来像我应该

If I have to implement it myself, I'll grumble but that's fine. If I do that, Wikipedia makes it sound like I should


  1. 生成 n 0均值,单位方差,独立正态样本(boost将执行此操作)

  2. - 协方差矩阵的分解

  3. 通过对应的特征值的平方根来缩放每个 n 样本

  4. 通过将缩放的矢量乘以分解中​​找到的正交特征向量矩阵,旋转样本向量

  1. generate n 0-mean, unit-variance, independent normal samples (boost will do this)
  2. find the eigen-decomposition of the covariance matrix
  3. scale each of the n samples by the square-root of the corresponding eigenvalue
  4. rotate the vector of samples by pre-multiplying the scaled vector by the matrix of orthonormal eigenvectors found by the decomposition

工作快。有人有一个直觉,什么时候值得检查协方差矩阵是否为正,如果是,使用Cholesky代替?

I would like this to work quickly. Does someone have an intuition for when it would be worthwhile to check to see if the covariance matrix is positive, and if so, use Cholesky instead?

推荐答案

由于这个问题已经获得了很多意见,我想我会发布代码的最终答案,我发现,在一定程度上,发布到Eigen论坛。代码使用Boost作为单变量法线和Eigen用于矩阵处理。它感觉相当不正常,因为它涉及使用内部命名空间,但它的工作原理。

Since this question has garnered a lot of views, I thought I'd post code for the final answer that I found, in part, by posting to the Eigen forums. The code uses Boost for the univariate normal and Eigen for matrix handling. It feels rather unorthodox, since it involves using the "internal" namespace, but it works. I'm open to improving it if someone suggests a way.

#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;
}

这篇关于来自C ++中多变量正态/高斯分布的样本的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!

查看全文
登录 关闭
扫码关注1秒登录
发送“验证码”获取 | 15天全站免登陆