Jak niezawodnie dodawać duże wyrażenia wykładnicze bez błędów przepełnienia?

Bardzo częstym problemem w Markov Chain Monte Carlo jest obliczanie prawdopodobieństw, które są sumą dużych terminów wykładniczych,

$e^{a_1} + e^{a_2} + ...$

gdzie składniki może wahać się od bardzo małych do bardzo dużych. Moje podejście na uwzględnieniu największego wykładniczego terminu , aby:$a$$K := \max_{i}(a_{i})$

e ' ≡ e się 1 + wiadomości e się 2 + . .

Takie podejście jest uzasadnione, jeżeli wszystkie elementy są duże, ale nie taki dobry pomysł, jeśli nie są. Oczywiście, mniejsze elementy i tak nie przyczyniają się do sumy zmiennoprzecinkowej, ale nie jestem pewien, jak sobie z nimi poradzić. W kodzie R moje podejście wygląda następująco:$a$

if ( max(abs(a)) > max(a) )
  K <-  min(a)
else
  K <- max(a)
ans <- log(sum(exp(a-K))) + K

Wydaje się dość powszechnym problemem, że powinno istnieć standardowe rozwiązanie, ale nie jestem pewien, co to jest. Dziękuję za wszelkie sugestie.

Istnieje proste rozwiązanie z dwoma przejściami danych:

Najpierw oblicz

co mówi ci, że jeśli nie ma terminów, to Σ i e i ≤ n e K$n$

Ponieważ prawdopodobnie nie masz tak blisko jak nawet , nie powinieneś się martwić przepełnieniem obliczeń z podwójną precyzją .10 20 τ : = ∑ i e a i - K ≤ $n$$10^{20}$

Zatem oblicz / a wtedy twoim rozwiązaniem jest .e K$\tau$$e^K \tau$

Aby zachować precyzję podczas sumowania podwójnego, musisz użyć Kahan Summation , jest to oprogramowanie równoważne z rejestrem przenoszenia.

Jest to w porządku dla większości wartości, ale jeśli dostajesz przepełnienie, to osiągasz limit podwójnej precyzji IEEE 754, co byłoby około . W tym momencie potrzebujesz nowej reprezentacji. Możesz wykryć przepełnienie w czasie dodawania, a także wykryć duże wykładniki do oceny . W tym momencie możesz zmodyfikować interpretację podwójnego, przesuwając wykładnik i śledząc tę ​​zmianę.$e^{709.783}$doubleMax - sumSoFar < valueToAddexponent > 709.783

To w przeważającej części jest podobne do twojego podejścia wykładniczego, ale ta wersja jest przechowywana w bazie 2 i nie wymaga wstępnego wyszukiwania w celu znalezienia największego wykładnika. Stąd .$\mathit{value}\times2^\mathit{shift}$

#!/usr/bin/env python
from math import exp, log, ceil

doubleMAX = (1.0 + (1.0 - (2 ** -52))) * (2 ** (2 ** 10 - 1))

def KahanSumExp(expvalues):
  expvalues.sort() # gives precision improvement in certain cases 
  shift = 0 
  esum = 0.0 
  carry = 0.0 
  for exponent in expvalues:
    if exponent - shift * log(2) > 709.783:
      n = ceil((exponent - shift * log(2) - 709.783)/log(2))
      shift += n
      carry /= 2*n
      esum /= 2*n
    elif exponent - shift * log(2) < -708.396:
      n = floor((exponent - shift * log(2) - -708.396)/log(2))
      shift += n
      carry *= 2*n
      esum *= 2*n
    exponent -= shift * log(2)
    value = exp(exponent) - carry 
    if doubleMAX - esum < value:
      shift += 1
      esum /= 2
      value /= 2
    tmp = esum + value 
    carry = (tmp - esum) - value 
    esum = tmp
  return esum, shift

values = [10, 37, 34, 0.1, 0.0004, 34, 37.1, 37.2, 36.9, 709, 710, 711]
value, shift = KahanSumExp(values)
print "{0} x 2^{1}".format(value, shift)

Twoje podejście jest solidne.

Nie musisz dokładnie znać , wystarczy, aby uniknąć przepełnienia. Być może będziesz w stanie oszacować analitycznie przed wykonaniem próbkowania MCMC.$K$$K$

Istnieje pakiet R, który zapewnia szybką i wydajną implementację „sztuczki log-sum-exp”

http://www.inside-r.org/packages/cran/matrixStats/docs/logSumExp

Funkcja logSumExp akceptuje wektor numeryczny lX i log danych wyjściowych (suma (exp (lX))), unikając problemów z przepełnieniem i przepełnieniem przy użyciu opisanej metody.