оценить многомерную нормальную/Gaussian плотности в 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: Обычно разложение Холески выполняет эту работу. Если ваша ковариантность действительно положительно определена (т. Е. Не закрывается с полуопределенной матрицей), разложите covMat = L*L^T и вычислите squaredNorm(L\(x-mu)) (где x=A\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().

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