在 C++ 中评估多元 Normal/Gaussian 密度

Question

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

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);
}

这只是抄写 formula。但是我得到了很多 0、nans 和 infs。我怀疑它来自非常接近于零的 covMat.determinant() 部分。

听说用x-meanVec预乘其协方差矩阵的"square root"的倒数更"stable"。从统计学上讲，这为您提供了一个随机向量，该向量的均值为零并且具有单位矩阵作为其协方差矩阵。我的问题是：

这真的是最好的方法吗？
这是 "best" 平方根技术，
我该怎么做？（最好使用 Eigen）

Answer 1

广告 1："Depends"。例如，如果你的协方差矩阵有一个特殊的结构，可以很容易地计算它的逆，或者如果维度非常小，它可以更快更稳定地实际计算逆。

广告 2：通常，Cholesky 分解可以完成这项工作。如果您的协方差确实是正定的（即不接近半定矩阵），请分解 covMat = L*L^T 并计算 squaredNorm(L\(x-mu))（其中 x=A\b 表示 "Solve A*x=b for x"）。当然，如果你的协方差是固定的，你应该只计算一次 L （也可能反转它）。您还应该使用 L 来计算 sqrt(covMat.determinant())，因为计算行列式否则需要再次分解 covMat。小改进：代替 pow(inv_sqrt_2pi, covMat.rows()) 计算 logSqrt2Pi=log(sqrt(2*pi)) 然后 return exp(-0.5*quadform - covMat.rows()*logSqrt2Pi) / L.determinant().

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

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++ 中评估多元 Normal/Gaussian 密度

evaluate multivariate Normal/Gaussian Density in c++

c++

eigen

normal-distribution