<p>I finally solved the problem and wrote a paper describing the solution in detail: <a href="https://arxiv.org/abs/1909.10765" rel="nofollow noreferrer">https://arxiv.org/abs/1909.10765</a></p>

<p>Briefly, divide and multiply each addend by the first term in the sum to obtain</p>

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

<p>where <em>ω(a, b, t, λ, μ)</em> is the first term in the series and <em>2F1</em> is the Gaussian hypergeometric function.
Hypergeometric function <em>2F1(-a, -b; -(a + b - k); -z)</em> (<em>a</em> and <em>b</em> positive integers, <em>k &lt;= 1</em>, <em>z</em> real number) can be computed with the following three-term recurrence relation (TTRR):</p>

<p><em>u(a, b, k) y(b + 1) + v(a, b, k, z) y(b) + w(b, k, z) y(b - 1) = 0</em></p>

<p>where</p>

<p><em>u(a, b, k) = (a + b + 1 − k) (a + b − k)</em></p>

<p><em>v(a, b, k, z) = − (a + b − k) (a + b + 1 − k + (a − b) z)</em></p>

<p><em>w(b, k, z) = − b (b − k) z</em></p>

<p>If <em>b > a</em> swap the two variables (that is <em>a' = max(a, b)</em> and <em>b' = min(a, b)</em>).</p>

<p>Compute the recurrence in a forward manner starting with values <em>y(0) = 2F1(-a, 0; -(a - k); -z) = 1</em> and <em>y(1) = 2F1(-a, -1; -(a + 1 - k); -z) = 1 + (a z) / (a + 1 - k)</em>.</p>

<p>I implemented the previous algorithm in the <em>Julia</em> package <a href="https://github.com/albertopessia/SimpleBirthDeathProcess.jl" rel="nofollow noreferrer">SimpleBirthDeathProcess</a>.</p>
How to evaluate an alternating series when addends contain rounding errors?
<p>I want to numerically evaluate the transition probability of a linear Birth and Death process</p>

<p><img src="https://i.stack.imgur.com/WA7Rh.png" height="60"></p>

<p>where <sub><img src="https://i.stack.imgur.com/QklZw.png" height="32"></sub> is the binomial coefficient and</p>

<p><img src="https://i.stack.imgur.com/FG177.png" height="90"></p>

<p>I am able to evaluate it with an acceptable numerical error (using logarithms and the Kahan-Neumaier summation algorithm) for most parameter combinations.</p>

<p>Problems arise when addends alternate in sign and numerical error dominates the sum (condition number tends to infinity in this case). This happens when</p>

<p><img src="https://i.stack.imgur.com/yiJpC.png" height="32"></p>

<p>For example, I have problems evaluating <code>p(1000, 2158, 72.78045, 0.02, 0.01)</code>. It should be 0 but I get the very big value <code>log(p) ≈ 99.05811</code>, which is impossible for a probability.</p>

<p>I tried refactoring the sum in many different ways and using various "precise" summation algorithms such as <a href="https://dl.acm.org/citation.cfm?id=1824815" rel="nofollow noreferrer">Zhu-Hayes</a>. 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.</p>

<p>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.</p>

<p>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.</p>

<p>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?</p>

<p>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.</p>



<pre class="lang-julia prettyprint-override"><code># 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 λ &gt; μ &gt; 0
function log_trans_prob(a, b, t, λ, μ)
  n = a + b
  k = min(a, b)

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

  if b &gt; 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 &lt; 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) &gt;= abs(val)) ? ((res - tot) + val) : ((val - tot) + res)
      res = tot
    end

    res += err

    if res &lt; 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 =&gt; 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))
</code></pre>