I want to numerically evaluate the transition probability of a linear Birth and Death process
where 
is the binomial coefficient and
I am able to evaluate it with an acceptable numerical error (using logarithms and the Kahan-Neumaier summation algorithm) for most parameter combinations.
Problems arise when addends alternate in sign and numerical error dominates the sum (condition number tends to infinity in this case). This happens when
For example, I have problems evaluating p(1000, 2158, 72.78045, 0.02, 0.01). It should be 0 but I get the very big value log(p) ≈ 99.05811, which is impossible for a probability.
I tried refactoring the sum in many different ways and using various "precise" summation algorithms such as Zhu-Hayes. I always get approximately the same wrong value, making me think that the problem is not the way I sum the numbers but the inner representation of each addend.
Because of binomial coefficients, values easily overflow. I tried with a linear transformation in order to keep each (absolute) element in the sum between the lowest normal number and 1. It didn't help and I think it is because of many algebraic operations of similar magnitudes.
I am now at a dead end and don't know how to proceed. I could use arbitrary precision arithmentic libraries but the computational cost is too high for my Markov Chain Monte Carlo application.
Is there a proper way or tricks to evaluate such sums when we cannot store partial sums at a good-enough precision in a IEEE-754 double?
Here is a basic working example where I only rescale the values by the maximum and sum with Kahan summation algorithm. Obviously, most values end up being subnormals with a 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))
