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

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乘以其协方差矩阵的“平方根”的逆数更为“稳定”。从统计上讲,这将为您提供均值为零的随机 vector ,并且具有单位矩阵作为其协方差矩阵。我的问题是:
  • 这真的是最好的方法吗?
  • 是“最佳”平方根技术,而
  • 我该怎么做? (最好使用Eigen)
  • 最佳答案

    广告1:“取决于”。例如,如果您的协方差矩阵具有可以轻松计算其逆数的特殊结构,或者如果维数很小,则实际计算逆数可能会更快,更稳定。

    广告2:通常,Cholesky分解可以完成这项工作。如果您的协方差是真正的正定的(即不接近半定矩阵),则分解covMat = L*L^T并计算squaredNorm(L\(x-mu))(其中x=A\b的意思是“为A*x=b解决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();
    }
    

    关于c++ - 在C++中评估多元正态/高斯密度,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/41538095/

    10-09 07:35