C++中多变量正态/高斯分布的样本

Sample from multivariate normal/Gaussian distribution in C++

本文关键字：高斯分布 样本 变量 C++      更新时间：2023-10-16

我一直在寻找一种从多元正态分布中采样的方便方法。有人知道有现成的代码片段可以做到这一点吗？对于矩阵/向量，我更喜欢使用Boost或Eigen或另一个我不熟悉的现象库，但我可以在必要时使用GSL。我也希望该方法接受非负-定协方差矩阵，而不是要求正定（例如，与Cholesky分解一样）。这在MATLAB、NumPy和其他软件中都存在，但我很难找到现成的C/C++解决方案。

如果我必须自己实施，我会抱怨，但没关系。如果我这样做，维基百科会让我觉得我应该

  1. 生成n0-均值、单位方差、独立正态样本（boost将执行此操作）
  2. 求协方差矩阵的特征分解
  3. 通过相应特征值的平方根对n样本中的每个样本进行缩放
  4. 通过将缩放后的向量与分解后的正交特征向量矩阵相乘，旋转样本向量

我希望这能很快奏效。有人有直觉知道什么时候值得检查协方差矩阵是否为正，如果是，请使用Cholesky？

由于这个问题已经获得了很多浏览量，我想我应该发布我找到的最终答案的代码，部分原因是通过发布到Eigen论坛。该代码将Boost用于单变量正态，将Eigen用于矩阵处理。它感觉相当非正统，因为它涉及到使用"内部"命名空间，但它是有效的。如果有人提出一种方法，我愿意改进它。

#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 << "Meann" << mean << std::endl;
  std::cout << "Covarn" << covar << std::endl;
  std::cout << "Samplesn" << samples << std::endl;
}

这里有一个在Eigen中生成多变量正态随机变量的类，它使用C++11随机数生成，并通过使用Eigen::MatrixBase::unaryExpr():避免了Eigen::internal的填充

struct normal_random_variable
{
    normal_random_variable(Eigen::MatrixXd const& covar)
        : normal_random_variable(Eigen::VectorXd::Zero(covar.rows()), covar)
    {}
    normal_random_variable(Eigen::VectorXd const& mean, Eigen::MatrixXd const& covar)
        : mean(mean)
    {
        Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> eigenSolver(covar);
        transform = eigenSolver.eigenvectors() * eigenSolver.eigenvalues().cwiseSqrt().asDiagonal();
    }
    Eigen::VectorXd mean;
    Eigen::MatrixXd transform;
    Eigen::VectorXd operator()() const
    {
        static std::mt19937 gen{ std::random_device{}() };
        static std::normal_distribution<> dist;
        return mean + transform * Eigen::VectorXd{ mean.size() }.unaryExpr([&](auto x) { return dist(gen); });
    }
};

它可以用作

int size = 2;
Eigen::MatrixXd covar(size,size);
covar << 1, .5,
        .5, 1;
normal_random_variable sample { covar };
std::cout << sample() << std::endl;
std::cout << sample() << std::endl;

对于现成的解决方案，armadillo C++库支持使用函数mvnrnd（）从多变量高斯分布（甚至从正半定协方差矩阵）进行采样。

做SVD然后检查矩阵是否为PD怎么样？请注意，这并不需要计算Cholskey因子分解。尽管，我认为SVD比Cholskey慢，但它们在失败次数上都必须是立方的。

相关文章：