는 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() 부분이 0에 매우 가깝다고 생각됩니다.

나는 공분산 행렬의 "제곱근"의 역수를 곱하여 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을 x"으로 해석합니다.)을 계산합니다. 물론 공분산이 고정 된 경우 L을 한 번만 계산해야합니다 (반전시킬 수도 있습니다). 행렬식을 계산하면 covMat을 다시 분해해야하므로 sqrt(covMat.determinant())도 계산하려면 L을 사용해야합니다. 사소한 개선 : pow(inv_sqrt_2pi, covMat.rows()) 대신 logSqrt2Pi=log(sqrt(2*pi))을 계산 한 다음 exp(-0.5*quadform - covMat.rows()*logSqrt2Pi)/L.determinant()을 반환하십시오.

광고 3 :이 나중에 아이겐 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(); 
} 

을