c++ - tabla - ejercicios resueltos de distribucion normal de probabilidades

Muestra de la distribución multivariable normal/gaussiana en C++

Question

He estado buscando una manera conveniente de muestrear desde una distribución normal multivariada. ¿Alguien sabe de un fragmento de código disponible para hacer eso? Para matrices / vectores, prefiero 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 aceptara matrices de covarianza no negativas -definitas en lugar de requerir un resultado positivo-definido (por ejemplo, como con la descomposición de Cholesky). Esto existe en MATLAB, NumPy y otros, pero me ha costado mucho encontrar una solución de C / C ++ ya preparada.

Si tengo que implementarlo yo mismo, me quejaré, pero eso está bien. Si hago eso, Wikipedia hace que suene como debería

1. generar n 0-media, unidad-varianza, muestras normales independientes (boost hará esto)
2. encontrar la descomposición propia de la matriz de covarianza
3. escalar cada una de las n muestras por la raíz cuadrada del valor propio correspondiente
4. 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 de cuándo valdría la pena comprobar si la matriz de covarianza es positiva y, de ser así, usar Cholesky en su lugar?

Answer 1

¿Qué hay de hacer una SVD y luego verificar si la matriz es PD? Tenga en cuenta que esto no requiere que calcule la factorización de Cholskey. Aunque, creo que la SVD es más lenta que Cholskey, pero ambas deben ser cúbicas en número de fracasos.

Answer 2

Aquí hay una clase para generar variables aleatorias normales multivariadas en Eigen que utiliza la generación de números aleatorios de C ++ 11 y evita las cosas Eigen::internal de Eigen::MatrixBase::unaryExpr() utilizando 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); });
    }
};

Puede ser utilizado 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;

Answer 3

Dado que esta pregunta ha acumulado muchas vistas, pensé que publicaría el código para la respuesta final que encontré, en parte, publicando en los foros de Eigen . El código utiliza Boost para el normal univariado y Eigen para el manejo de matrices. 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;
}