如何解决来自C ++中多元正态/高斯分布的样本
|| 我一直在寻找从多元正态分布采样的便捷方法。有谁知道一个现成的代码片段可以做到这一点？对于矩阵/向量，我更喜欢使用Boost或Eigen或其他我不熟悉的非凡的库，但是我可以适时使用GSL。如果该方法接受非负定协方差矩阵而不是需要正定（例如，如Cholesky分解），我也很喜欢它。这在MATLAB，NumPy和其他工具中都存在，但是我很难找到现成的C / C ++解决方案。 如果我必须自己实现它，我会抱怨，但这很好。如果我这样做了，维基百科让我听起来应该 生成n个0均值，单位方差的独立正常样本（boost将这样做） 求出协方差矩阵的特征分解 按相应特征值的平方根缩放n个样本中的每个样本 通过将缩放后的向量预先乘以分解发现的正交本征向量矩阵来旋转样本向量 我希望这个工作很快。是否有人应该何时去检查协方差矩阵是否为正，是否有直觉？如果是，则改用Cholesky？解决方法
由于这个问题已经引起了很多见解，所以我认为我将为最终找到的答案张贴代码，部分方法是在Eigen论坛上发帖。该代码使用Boost表示单变量正态，使用Eigen表示矩阵。感觉很不合常规，因为它涉及使用\“ internal \”名称空间，但是它可以工作。如果有人建议，我愿意改进它。
#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,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;
}
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,1;
normal_random_variable sample { covar };
std::cout << sample() << std::endl;
std::cout << sample() << std::endl;
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。