我想以数值方式评估线性生死过程的转移概率
其中
是二项式系数和
对于大多数参数组合，我能够以可接受的数值误差（使用对数和 Kahan-Neumaier 求和算法）对其进行评估。
当加数在符号上交替出现并且数值误差占总和时会出现问题（在这种情况下，条件数趋于无穷大）。这发生在
例如，我在评估p(1000, 2158, 72.78045, 0.02, 0.01). 它应该是 0 但我得到的值非常大log(p) ≈ 99.05811，这对于概率来说是不可能的。
我尝试以许多不同的方式重构总和，并使用各种“精确”求和算法，例如Zhu-Hayes。我总是得到大致相同的错误值，这让我认为问题不在于我对数字求和的方式，而是每个加数的内部表示。
由于二项式系数，值很容易溢出。我尝试进行线性变换，以使每个（绝对）元素保持在最低正态数和 1 之间的总和中。它没有帮助，我认为这是因为许多类似量级的代数运算。
我现在处于死胡同，不知道如何进行。我可以使用任意精度的算术库，但是对于我的马尔可夫链蒙特卡罗应用程序来说，计算成本太高了。
当我们无法在 IEEE-754 双精度中以足够好的精度存储部分总和时，是否有适当的方法或技巧来评估这些总和？
这是一个基本的工作示例，我仅按最大值重新调整值并使用 Kahan 求和算法求和。显然，大多数值最终都是 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))