在 C++ 中评估多元正态/高斯密度

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 与其协方差矩阵的“平方根”的倒数预乘起来更“稳定”。从统计上讲，这为您提供了一个均值为零的随机向量，并将单位矩阵作为其协方差矩阵。我的问题是：

1. 这真的是最好的方法吗？
2. 哪个是“最好的”平方根技术，以及
3. 我该怎么做？ （最好使用 Eigen）

Answer 1

广告 1：“视情况而定”。例如，如果您的协方差矩阵具有特殊结构，可以轻松计算其逆矩阵，或者如果维度非常小，则 可以更快、更稳定地实际计算逆矩阵。

广告 2：通常，Cholesky 分解可以完成这项工作。如果您的协方差确实是正定的（即，不接近半定矩阵），请分解 covMat = L*L^T 并计算 squaredNorm(L\(x-mu))（其中 x=A\b 表示“求解 A*x=b for x”）。当然，如果你的协方差是固定的，你应该只计算一次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()。

广告 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();
}