Как оценить чередующийся ряд, если в аддентах есть ошибки округления?

Question

Я хочу, чтобы численно оценить вероятность перехода линейного процесса рождения и смерти

где является биномиальным коэффициентом и

Я могу оценить его с приемлемой числовой ошибкой (используя логарифмы и алгоритм суммирования Кахана-Неймайера) для большинства комбинаций параметров.

Проблемы возникают, когда слагаемые чередуются по знаку и числовая ошибка доминирует в сумме (в этом случае число условия стремится к бесконечности). Это происходит когда

Например, у меня проблемы с оценкой p(1000, 2158, 72.78045, 0.02, 0.01), Должно быть 0, но я получаю очень большое значение log(p) ≈ 99.05811, что невозможно по вероятности.

Я попытался провести рефакторинг суммы разными способами и с использованием различных "точных" алгоритмов суммирования, таких как Чжу-Хейс. Я всегда получаю примерно одно и то же неправильное значение, заставляя меня думать, что проблема не в способе суммирования чисел, а в внутреннем представлении каждого добавления.

Из-за биномиальных коэффициентов значения легко переполняются. Я попытался с помощью линейного преобразования, чтобы сохранить каждый (абсолютный) элемент в сумме между наименьшим нормальным числом и 1. Это не помогло, и я думаю, что это из-за многих алгебраических операций схожих величин.

Я сейчас в тупике и не знаю, как поступить. Я мог бы использовать арифметические библиотеки произвольной точности, но вычислительные затраты слишком высоки для моего приложения "Марковская цепь Монте-Карло".

Есть ли правильный способ или приемы для оценки таких сумм, когда мы не можем хранить частичные суммы с достаточно хорошей точностью в двойном IEEE-754?

Вот основной рабочий пример, где я только изменяю значения по максимуму и суммирую с помощью алгоритма суммирования Кахана. Очевидно, что большинство значений в конечном итоге являются субнормальными с Float64.

# this is the logarithm of the absolute value of element h
@inline function log_addend(a, b, h, lα, lβ, lγ)
  log(a) + lgamma(a + b - h) - lgamma(h + 1) - lgamma(a - h + 1) -
  lgamma(b - h + 1) + (a - h) * lα + (b - h) * lβ + h * lγ
end

# this is the logarithm of the ratio between (absolute) elements i and j
@inline function log_ratio(a, b, i, j, q)
  lgamma(j + 1) + lgamma(a - j + 1) + lgamma(b - j + 1) + lgamma(a + b - i) -
  lgamma(i + 1) - lgamma(a - i + 1) - lgamma(b - i + 1) - lgamma(a + b - j) +
  (j - i) * q
end

# function designed to handle the case of an alternating series with λ > μ > 0
function log_trans_prob(a, b, t, λ, μ)
  n = a + b
  k = min(a, b)

  ω = μ / λ
  η = exp((μ - λ) * t)

  if b > zero(b)
    lβ = log1p(-η) - log1p(-ω * η)
    lα = log(μ) + lβ - log(λ)
    lγ = log(ω - η) - log1p(-ω * η)
    q = lα + lβ - lγ

    # find the index of the maximum addend in the sum
    # use a numerically stable method for solving quadratic equations
    x = exp(q)
    y = 2 * x / (1 + x) - n
    z = ((b - x) * a - x * b) / (1 + x)

    sup = if y < zero(y)
      ceil(typeof(a), 2 * z / (-y + sqrt(y^2 - 4 * z)))
    else
      ceil(typeof(a), (-y - sqrt(y^2 - 4 * z)) / 2)
    end

    # Kahan summation algorithm
    val = zero(t)
    tot = zero(t)
    err = zero(t)
    res = zero(t)
    for h in 0:k
      # the problem happens here when we call the `exp` function
      # My guess is that log_ratio is either very big or very small and its
      # `exp` cannot be properly represented by Float64
      val = (-one(t))^h * exp(log_ratio(a, b, h, sup, q))
      tot = res + val
      # Neumaier modification
      err += (abs(res) >= abs(val)) ? ((res - tot) + val) : ((val - tot) + res)
      res = tot
    end

    res += err

    if res < zero(res)
      # sum cannot be negative (they are probabilities), it might be because of
      # rounding errors
      res = zero(res)
    end

    log_addend(a, b, sup, lα, lβ, lγ) + log(res)
  else
    a * (log(μ) + log1p(-η) - log(λ) - log1p(-ω * η))
  end
end

# ≈ 99.05810564477483 => impossible
log_trans_prob(1000, 2158, 72.78045, 0.02, 0.01)

# increasing precision helps but it is too slow for applications
log_trans_prob(BigInt(1000), BigInt(2158), BigFloat(72.78045), BigFloat(0.02),
               BigFloat(0.01))

Answer 1

Я наконец решил проблему и написал статью, подробно описывающую решение: https://arxiv.org/abs/1909.10765

Вкратце, разделите и умножьте каждое слагаемое на первое слагаемое в сумме, чтобы получить

p(a, b, t, λ, μ) = ω(a, b, t, λ, μ) 2F1(-a, -b; -(a + b - 1); -z(t, λ, μ))

где ω(a, b, t, λ, μ) - первое слагаемое в ряду, а 2F1 - гипергеометрическая функция Гаусса. Гипергеометрическая функция 2F1(-a, -b; -(a + b - k); -z) (положительные целые числаa и b, k <= 1, z действительное число) может быть вычислена с помощью следующего трехчленного рекуррентного соотношения (TTRR):

u(a, b, k) y(b + 1) + v(a, b, k, z) y(b) + w(b, k, z) y(b - 1) = 0

где

u (a, b, k) = (a + b + 1 - k) (a + b - k)

v (a, b, k, z) = - (a + b - k) (a + b + 1 - k + (a - b) z)

w (b, k, z) = - b (b - k) z

Если b > a поменять местами две переменные (то есть a' = max(a, b) и b' = min (a, b)).

Вычислите повторение прямым способом, начиная со значений y(0) = 2F1(-a, 0; -(a - k); -z) = 1 и y(1) = 2F1(-a, -1; -(а + 1 - к); -z) = 1 + (аз) / (а + 1 - к).

Я реализовал предыдущий алгоритм в пакете Julia SimpleBirthDeathProcess.