線形の生死過程の遷移確率を数値的に評価したい
どこ
で 二項係数と
ほとんどのパラメーターの組み合わせについて、(対数と Kahan-Neumaier 加算アルゴリズムを使用して) 許容可能な数値誤差で評価できます。
加数の符号が交互になり、数値誤差が合計を支配する場合に問題が発生します (この場合、条件数は無限大になる傾向があります)。これは次の場合に発生します。
たとえば、評価に問題がありますp(1000, 2158, 72.78045, 0.02, 0.01)。0 のはずですが、非常に大きな値log(p) ≈ 99.05811が得られました。これは、確率としては不可能です。
さまざまな方法で合計をリファクタリングし、Zhu-Hayesなどのさまざまな「正確な」合計アルゴリズムを使用してみました。私はいつもほぼ同じ間違った値を取得しているため、問題は数値を合計する方法ではなく、各加数の内部表現にあると思います。
二項係数のため、値は簡単にオーバーフローします。各 (絶対) 要素を最小正規数と 1 の間の合計に保つために、線形変換を試みました。
私は今行き止まりにあり、先に進む方法がわかりません。任意精度の算術ライブラリを使用できますが、マルコフ連鎖モンテカルロ アプリケーションには計算コストが高すぎます。
IEEE-754 double に十分な精度で部分和を格納できない場合に、そのような和を評価する適切な方法またはトリックはありますか?
これは、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))