在C ++中评估多元正态/高斯密度 [英] evaluate multivariate Normal/Gaussian Density in c++

问题描述

现在，我具有以下函数来评估高斯密度：

Right now I have the following function to evaluate the Gaussian density:

double densities::evalMultivNorm(const Eigen::VectorXd &x, const Eigen::VectorXd &meanVec, const Eigen::MatrixXd &covMat)
{
    double inv_sqrt_2pi = 0.3989422804014327;
    double quadform  = (x - meanVec).transpose() * covMat.inverse() * (x-meanVec);
    double normConst = pow(inv_sqrt_2pi, covMat.rows()) * pow(covMat.determinant(), -.5);
    return normConst * exp(-.5* quadform);
}

这只是转录公式。但是我得到很多0，nans和infs。我怀疑它来自 covMat.determinant（）部分非常接近于零。

This is just transcribing the formula. However I get a lot of 0, nans and infs. I suspect it's coming from the covMat.determinant() portion being very close to zero.

我听说将 x-meanVec 乘以其协方差矩阵的平方根的逆数会更稳定。从统计上讲，这将为您提供均值为零的随机向量，并且具有单位矩阵作为其协方差矩阵。我的问题是：

I have heard that it is more "stable" to premultiply x-meanVec by the inverse of a "square root" of its covariance matrix. Statistically this gives you a random vector that is mean zero and has the identity matrix as its covariance matrix. My questions are:

1. 这真的是最好的方法吗？


2. 最佳的平方根技巧，以及


3. 我该怎么做？ （最好使用Eigen）

推荐答案

广告1：依赖。例如，如果协方差矩阵具有一种特殊的结构，可以轻松地计算其逆数，或者如果维数很小，则 可以更快，更稳定地实际计算逆数。

Ad 1: "Depends". E.g., if your covariance matrix has a special structure which makes it easy to calculate its inverse or if the dimension is very small, it can be faster and more stable to actually calculate the inverse.

广告2：通常，霍尔斯基分解可完成此工作。如果您的协方差是真正的正定的（即不接近半定矩阵），则分解 covMat = L * L ^ T 并计算 squaredNorm（ L\（x-mu））（其中 x = A\b 的意思是解决 A * x = b 表示 x ）。当然，如果您的协方差是固定的，则您应该只计算一次（code> L （也许也可以求反）。您还应该使用 L 来计算 sqrt（covMat.determinant（）），因为计算行列式会另外需要再次分解 covMat 。
较小改进：代替 pow（inv_sqrt_2pi，covMat.rows（））计算 logSqrt2Pi = log（sqrt（2 * pi）） ，然后返回 exp（-0.5 * quadform-covMat.rows（）* logSqrt2Pi）/ L.determinant（）。

Ad 2: Usually, Cholesky decomposition does the job. If your covariance is truly positive definite (i.e., not to close to a semidefinite matrix), decompose covMat = L*L^T and compute squaredNorm(L\(x-mu)) (where x=A\b means "Solve A*x=b for x"). Of course, if your covariance is fixed, you should compute L only once (and perhaps invert it as well). You should use L to compute sqrt(covMat.determinant()) as well, since computing the determinant would otherwise require to decompose covMat again. Minor improvement: instead of pow(inv_sqrt_2pi, covMat.rows()) compute logSqrt2Pi=log(sqrt(2*pi)) and then return exp(-0.5*quadform - covMat.rows()*logSqrt2Pi) / L.determinant().

广告3：应在Eigen 3.2或更高版本中运行：

Ad 3: This should run in Eigen 3.2 or later:

double foo(const Eigen::VectorXd &x, const Eigen::VectorXd &meanVec, const Eigen::MatrixXd &covMat)
{
    // avoid magic numbers in your code. Compilers will be able to compute this at compile time:
    const double logSqrt2Pi = 0.5*std::log(2*M_PI);
    typedef Eigen::LLT<Eigen::MatrixXd> Chol;
    Chol chol(covMat);
    // Handle non positive definite covariance somehow:
    if(chol.info()!=Eigen::Success) throw "decomposition failed!";
    const Chol::Traits::MatrixL& L = chol.matrixL();
    double quadform = (L.solve(x - meanVec)).squaredNorm();
    return std::exp(-x.rows()*logSqrt2Pi - 0.5*quadform) / L.determinant();
}

这篇关于在C ++中评估多元正态/高斯密度的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！