Muestra de distribución multivariante normal/gaussiana en C++
He estado buscando una forma conveniente de muestrear de una distribución normal multivariante. ¿Alguien sabe de un fragmento de código fácilmente disponible para hacer eso? Para matrices/vectores, preferiría usar Boost o Eigen u otra biblioteca fenomenal con la que no estoy familiarizado, pero podría usar GSL en un apuro. También me gustaría si el método aceptado no negativo-matrices de covarianza definidas en lugar de requerir positivo-definido (por ejemplo, como con el Cholesky descomposición). Esto existe en MATLAB, NumPy y otros, pero he tenido dificultades para encontrar una solución C/C++ lista para usar.
Si tengo que implementarlo yo mismo, me quejaré, pero está bien. Si hago eso, Wikipedia hace que suene como debería
- generar n 0-media, unidad de varianza, muestras normales independientes (boost hará esto)
- encontrar la descomposición propia de la matriz de covarianza
- escala cada una de las muestras n por la raíz cuadrada del valor propio correspondiente
- gire el vector de muestras multiplicando previamente el vector escalado por la matriz de vectores propios ortonormales encontrados por la descomposición
Me gustaría que esto funcione rápidamente. ¿Alguien tiene una intuición para cuándo valdría la pena verificar si la matriz de covarianza es positiva, y si es así, use Cholesky en su lugar?
3 answers
Dado que esta pregunta ha obtenido muchas opiniones, pensé en enviar un código para la respuesta final que encontré, en parte, mediante publicar en los foros de Eigen. El código utiliza Boost para la normal univariada y Eigen para el manejo de la matriz. Se siente bastante poco ortodoxo, ya que implica el uso del espacio de nombres" interno", pero funciona. Estoy abierto a mejorarlo si alguien sugiere una manera.
#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;
}
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2012-12-26 21:22:01
Aquí hay una clase para generar variables aleatorias normales multivariantes en Eigen que usa C++11 generación de números aleatorios y evita las cosas Eigen::internal usando Eigen::MatrixBase::unaryExpr():
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); });
}
};
Se puede utilizar como
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;
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2017-01-18 14:25:45
¿Qué hay de hacer un SVD y luego comprobar si la matriz es PD? Tenga en cuenta que esto no requiere que calcule la factorización Cholskey. Aunque, creo que SVD es más lento que Cholskey, pero ambos deben ser cúbicos en número de flops.
Warning: date(): Invalid date.timezone value 'Europe/Kyiv', we selected the timezone 'UTC' for now. in /var/www/agent_stack/data/www/ajaxhispano.com/template/agent.layouts/content.php on line 61
2011-06-01 16:48:45