No momento, tenho a seguinte função para avaliar a densidade gaussiana:

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

Isso está apenas transcrevendo oFórmula. No entanto, recebo muitos 0, nans e infs. Eu suspeito que está vindo docovMat.determinant() porção sendo muito próxima de zero.

Ouvi dizer que é mais "estável" pré-multiplicarx-meanVec pelo inverso de uma "raiz quadrada" de sua matriz de covariância. Estatisticamente, isso fornece um vetor aleatório com média zero e tem a matriz de identidade como sua matriz de covariância. Minhas perguntas são:

Este é realmente o melhor caminho a percorrer?Qual é a "melhor" técnica de raiz quadrada eComo eu faço isso? (De preferência usando Eigen)