Образец из многомерного нормального / гауссовского распределения в C ++ - PullRequest
Новая
       75

Образец из многомерного нормального / гауссовского распределения в C ++

30 голосов
/ 26 мая 2011

Я искал удобный способ выборки из многомерного нормального распределения.Кто-нибудь знает готовый фрагмент кода для этого?Для матриц / векторов я бы предпочел использовать Boost или Eigen или другую феноменальную библиотеку, с которой я не знаком, но я мог бы использовать GSL взажать.Я также хотел бы, чтобы метод принимал неотрицательные -определенные ковариационные матрицы вместо того, чтобы требовать положительно определенной (например, как с разложением Холецкого).Это существует в MATLAB, NumPy и других, но мне было трудно найти готовое решение C / C ++.

Если мне придется реализовать это самостоятельно, я буду ворчать, но это нормально.Если я сделаю это, Википедия заставит это звучать так, как я должен

  1. генерировать n 0-среднее, единичная дисперсия, независимые нормальные сэмплы (повышение сделаетthis)
  2. найти собственное разложение ковариационной матрицы
  3. масштабировать каждый из n отсчетов по корню квадратному из соответствующего собственного значения
  4. поверните вектор отсчетов, предварительно умножив масштабированный вектор на матрицу ортонормированных собственных векторов, найденных по разложению

Я бы хотел, чтобы это работало быстро.Есть ли у кого-то интуиция, когда стоит проверить, положительна ли ковариационная матрица, и если да, то использовать вместо нее Cholesky?

Ответы [ 3 ]

23 голосов
/ 27 декабря 2012

Поскольку этот вопрос получил много просмотров, я подумал, что напишу код для окончательного ответа, который я нашел, частично, путем публикации на форумах 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 << "Mean\n" << mean << std::endl;
  std::cout << "Covar\n" << covar << std::endl;
  std::cout << "Samples\n" << samples << std::endl;
}
9 голосов
/ 25 октября 2016

Вот класс для генерации многомерных нормальных случайных величин в Eigen, который использует генерацию случайных чисел C ++ 11 и избегает вещей Eigen::internal, используя 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); });
    }
};

Может использоваться как 

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;
0 голосов
/ 01 июня 2011

Как насчет выполнения SVD, а затем проверки, является ли матрица PD? Обратите внимание, что для этого не требуется вычислять факторизацию Чольски. Хотя я думаю, что SVD медленнее, чем Cholskey, но они оба должны быть кубическими по количеству флопов.

...