c++ - C++中多元正态/高斯分布的样本

我一直在寻找从多元正态分布采样的便捷方法。有谁知道一个现成的代码片段可以做到这一点？对于矩阵/vector ，我更喜欢使用Boost或Eigen或其他我不熟悉的非凡的库，但可以在适当的时候使用GSL。如果该方法接受非负定协方差矩阵而不是要求正定(例如，如Cholesky分解)，我也很喜欢它。这在MATLAB，NumPy等中都存在，但是我很难找到现成的C/C++解决方案。

如果我必须自己实现它，我会提示，但这很好。如果我这样做，就应该像我一样Wikipedia makes it sound

生成 n 0均值，单位方差，独立的正常样本(增强将执行此操作)

找到协方差矩阵

的特征分解

通过对应特征值

的平方根来缩放每个和样本

通过将缩放后的 vector 与分解

找到的正交本征 vector 矩阵预先乘以旋转样本 vector

我希望这个工作很快。是否有人有直觉，何时应该检查协方差矩阵是否为正，如果是，则改用Cholesky？

最佳答案

既然这个问题引起了很多见解，所以我认为我应该为posting to the Eigen forums找到的最终答案贴上代码。该代码使用Boost表示单变量正态，使用Eigen表示矩阵。感觉很不合常规，因为它涉及使用“内部” namespace ，但是它可以工作。如果有人提出建议，我愿意改进它。

#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++ - C++中多元正态/高斯分布的样本，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/6142576/