Как надежно добавить большие экспоненциальные термины без ошибок переполнения?

Очень распространенная проблема в цепочке Маркова Монте-Карло включает вычисление вероятностей, которые являются суммой больших экспоненциальных членов,

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

$a$$K := \max_{i}(a_{i})$

е ' ≡ е в 1 + Ē в 2 + . ,

Такой подход является разумным , если все элементы велики, но не такая хорошая идея , если они не являются. Конечно, меньшие элементы в любом случае не влияют на сумму с плавающей запятой, но я не уверен, как с ними надежно справиться. В коде R мой подход выглядит так:$a$

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

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

Существует простое решение только с двумя проходами через данные:

Сначала вычислите

который говорит вам, что, если есть терминов, то Σ я е я ≤ п е К$n$

Так как вы, вероятно, не имеете столь же большого, как даже , вам не нужно беспокоиться о переполнении в вычислениях с двойной точностью ,$n$$10^{20}$

Таким образом, вычислите и тогда ваше решение будет .$\tau$$e^K \tau$

Чтобы сохранить точность при добавлении двойных чисел, вам необходимо использовать Kahan Summation , это программное обеспечение, эквивалентное наличию регистра переноса.

$e^{709.783}$doubleMax - sumSoFar < valueToAddexponent > 709.783

$\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)

Ваш подход твердый.

$K$$K$

Существует пакет R, который обеспечивает быструю и эффективную реализацию "трюка log-sum-exp"

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

Функция logSumExp принимает числовой вектор lX и выводит log (sum (exp (lX))), избегая при этом проблем переполнения и переполнения, используя описанный вами метод.