首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >c++中多元正态/高斯密度的评价

c++中多元正态/高斯密度的评价
EN

Stack Overflow用户
提问于 2017-01-08 21:28:13
回答 1查看 2.7K关注 0票数 1

现在,我有以下函数来计算高斯密度:

代码语言:javascript
复制
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()部分,非常接近于零。

我听说,用x-meanVec协方差矩阵的“平方根”的逆乘积它更“稳定”。统计上,这给出了一个均值为零的随机向量,并将恒等矩阵作为它的协方差矩阵。我的问题是:

  1. 这真的是最好的方法吗?
  2. 这是“最好”的平方根技术
  3. 我该怎么做呢?(最好使用特征)
EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2017-01-09 13:17:34

广告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或更高版本中运行:

代码语言:javascript
复制
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();
}
票数 3
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/41538095

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档