Bagaimana cara menambahkan istilah eksponensial besar secara andal tanpa kesalahan luapan?

Masalah yang sangat umum di Markov Chain Monte Carlo melibatkan probabilitas komputasi yang merupakan jumlah dari istilah eksponensial besar,

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

dimana komponen dapat berkisar dari sangat kecil sampai yang sangat besar. Pendekatan saya adalah memperhitungkan faktor eksponensial terbesar sehingga:$a$$K := \max_{i}(a_{i})$

Pendekatan ini masuk akal jika semua elemen $a$ besar, tetapi bukan ide yang bagus jika tidak. Tentu saja, elemen yang lebih kecil tidak berkontribusi pada jumlah floating-point, tapi saya tidak yakin bagaimana cara mengatasinya dengan andal. Dalam kode R, pendekatan saya terlihat seperti:

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

Tampaknya masalah yang cukup umum bahwa harus ada solusi standar, tetapi saya tidak yakin apa itu. Terima kasih atas sarannya.

Ada solusi langsung dengan hanya dua melewati data:

pertama

yang memberitahu Anda bahwa, jika ada istilah, maka $n$

Karena Anda mungkin tidak memiliki mendekati hampir , Anda tidak perlu khawatir meluap dalam perhitungan dalam presisi ganda .$n$$10^{20}$

Dengan demikian, hitung dan kemudian solusi Anda adalah .$\tau$$e^K \tau$

Untuk menjaga ketepatan saat Anda menambahkan ganda bersama-sama Anda perlu menggunakan Penjumlahan Kahan , ini adalah perangkat lunak yang setara dengan membawa register.

Ini bagus untuk sebagian besar nilai, tetapi jika Anda mendapatkan overflow maka Anda mencapai batas presisi ganda IEEE 754 yang sekitar . Pada titik ini Anda memerlukan representasi baru. Anda dapat mendeteksi kelebihan pada waktu tambahan dengan dan juga mendeteksi eksponen ke besar untuk dievaluasi oleh . Pada titik ini Anda dapat memodifikasi interpretasi ganda dengan menggeser eksponen dan melacak pergeseran ini.$e^{709.783}$doubleMax - sumSoFar < valueToAddexponent > 709.783

Ini untuk sebagian besar adalah mirip dengan pendekatan Anda mengimbangi eksponen, tetapi versi ini disimpan di basis 2 dan tidak memerlukan pencarian awal untuk menemukan eksponen terbesar. Karenanya .$\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)

Pendekatan Anda solid.

Anda tidak perlu tahu persis , cukup baik untuk menghindari overflow. Jadi, Anda mungkin dapat memperkirakan analitis sebelum Anda melakukan pengambilan sampel MCMC.$K$$K$

Ada paket R yang memasok implementasi yang cepat dan efisien dari "trik log-sum-exp"

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

Fungsi logSumExp menerima vektor numerik lX dan output log (jumlah (exp (lX))) sambil menghindari masalah underflow dan overflow menggunakan metode yang Anda jelaskan.