Численно стабильный способ вычислить нормальное логарифмическое правдоподобие

Question

Я пытаюсь вычислить нормальное логарифмическое правдоподобие, который определяется по формуле:

L = l1+l2+l3+...+ln, 

где

lk = log(1/(sqrt(2*PI)*sigma_k))-0.5*e_k*e_k 

Сигмы вокруг 0.2 и e_k нормально распределены со средним значением 0 и единичной дисперсией, поэтому большинство из них находятся между -2 и 2;

Я попытался следующий код Java (sigma_k упоминалось выше = sigmas.get (к) * Math.sqrt (дт)):

private double new1(List<Double> residuals, List<Double> sigmas, double dt) { 
    double a = 0; 
    for(int i=0; i<sigmas.size(); i++) { 
     a += Math.log(1.0/(Math.sqrt(2*Math.PI*dt)*sigmas.get(i))); 
    } 
    double b = 0; 
    for(int i=0; i<residuals.size(); i++) { 
     b += residuals.get(i)*residuals.get(i); 
    } 
    return a-0.5*b; 
} 

но теоретический максимум ниже максимума меня, делая численная оптимизация, поэтому у меня есть некоторые подозрения, что мой метод субоптимален.

Answer 1

Я не уверен, что это позволит значительно улучшить численную устойчивость, но ваши уравнения можно упростить, используя логарифм законы:

log(a*b) = log(a) + log(b) 
log(1/a) = -log(a) 
log(sqrt(a)) = log(a)/2 

поэтому у вас есть:

lk = -log(2*pi)/2 - log(sigma_k) - 0.5*e_k*e_k 
    = -log(2*pi)/2 - log(dt)/2 - log(sigmas.get(k)) - 0.5*e_k*e_k 
    = -log(2*pi*dt)/2 - log(sigmas.get(k)) - 0.5*e_k*e_k 

Первая постоянна, так в первом цикле вам нужно только сделать a -= log(sigmas.get(k)).

Кроме того, это выглядит подозрительно, что первый цикл равен sigmas.size(), а второй - residuals.size(), а уравнение предлагает, чтобы они имели одинаковую длину.

Answer 2

Примечание: В некоторой вероятности областей/статистика рассчитывается, не принимая журнала, например, в языковых частотах комбинаций.

Следующее упрощает, становясь менее стабильным, но впоследствии преобразует его обратно в сумму журналов или около того.

---

double a = 0; 
for(int i=0; i<sigmas.size(); i++) { 
    a += Math.log(1.0/(Math.sqrt(2*Math.PI*dt)*sigmas.get(i))); 
} 

---

журнал (х) + журнал (у) = лог (х * у)

double a = 1.0; 
for(int i=0; i<sigmas.size(); i++) { 
    a *= 1.0/(Math.sqrt(2*Math.PI*dt)*sigmas.get(i)); 
} 
a = Math.log(a); 

---

(1/х) * (1/y) = 1/(x * y)

double a = 1.0; 
for(int i=0; i<sigmas.size(); i++) { 
    a *= Math.sqrt(2*Math.PI*dt)*sigmas.get(i); 
} 
a = Math.log(1.0/a); 

---

SQRT (х)^п = (х^0,5)^п = х^(п/2)

static import Math.*; 

double a = pow(2*PI*dt, sigmas.size()/2.0); 
for(int i=0; i<sigmas.size(); i++) { 
    a *= sigmas.get(i); 
} 
a = -log(a);