来自C ++中多元正态/高斯分布的样本
||
我一直在寻找从多元正态分布采样的便捷方法。有谁知道一个现成的代码片段可以做到这一点？对于矩阵/向量，我更喜欢使用Boost或Eigen或其他我不熟悉的非凡的库，但是我可以适时使用GSL。如果该方法接受非负定协方差矩阵而不是需要正定（例如，如Cholesky分解），我也很喜欢它。这在MATLAB，NumPy和其他工具中都存在，但是我很难找到现成的C / C ++解决方案。
如果我必须自己实现它，我会抱怨，但这很好。如果我这样做了，维基百科让我听起来应该
生成n个0均值，单位方差的独立正常样本（boost将这样做）
求出协方差矩阵的特征分解
按相应特征值的平方根缩放n个样本中的每个样本
通过将缩放后的向量预先乘以分解发现的正交本征向量矩阵来旋转样本向量
我希望这个工作很快。是否有人应该何时去检查协方差矩阵是否为正，是否有直觉？如果是，则改用Cholesky？
没有找到相关结果
已邀请:
3 个回复
缔恃钨
#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; }烫珊
避免了东西：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;视蕉梁拌客