My question is specific to algorithms and models of computation.

I would like to write code to evaluate the following expression quickly and accurately:

$$\log \left( \sum_{i=1}^{\infty}{I_{\nu+i}(2\lambda)} \right)$$

Where $I_\nu(x)$ is the modified bessel function of the first kind.

I have limited experience with numerical summation of infinite series. Is there a general approach/algorithm for this? Kahan summation algorithm? Some other general purpose approach? Any reference material to help me learn are greatly appreciated.

So far my only ideas are to either iteratively sum to N and stop when $\log \left( \sum_{i=1}^{N}{I_{\nu+i}(2\lambda)} \right)$ meets some convergence criteria (though I am not certain how precise it ought to be) or I might try summing $\log \left( \sum_{i=N}^{\infty}{I_{\nu+i}(2\lambda)} \right)$ using a recurrence relation $I_{\nu + 1}(x) = I_{\nu-1}(x) - \frac{2\nu}{x} I_{\nu}(x)$ until the tail value is small enough, and then calculating the partial sum up to $N$.

For those interested, I am working in R.

UPDATE: Based on the comments, I have a few followup points.

  1. Convergence: I am working with the SDF of the absolute value of a Skellam random variable, i.e $P(|W|>\nu)$. In probability theory, the SDF is always between 0 and 1, and it can be shown that the SDF is equal to:

$$P(|W|>\nu) = 2\exp(-2\lambda)\sum_{i=1}^{\infty}{I_{\nu+i}(2\lambda)}, \ \ \lambda > 0, \nu \in \{0,1,2,3, \dots , \infty\}$$


$$ 0 < P(|W|>\nu) \le 1 \\ 0 < 2\exp(-2\lambda)\sum_{i=1}^{\infty}{I_{\nu+i}(2\lambda)} \le 1 \\ 0 < \sum_{i=1}^{\infty}{I_{\nu+i}(2\lambda)} \le \frac{\exp(2\lambda)}{2} $$

demonstrating that the infinite sum does indeed converge.

  1. Range for $\lambda$ and $\nu$: As mentioned in part (1), $\lambda$ can be any real number greater than zero and $\nu$ (the random variable) can take on any value in the non-negative integers. Since we are dealing with probabilities, we will assume $\lambda$ is fixed and attempt to compute the infinite sum for any value of $\nu$. I would PREFER an approach to calculating the value that is general for any fixed value of $\lambda$ and multiple $\nu$, but for argument's sake, anyone can take $\lambda = 2000$ and $\nu = 1700$

  2. Convergence rate: I am not a computer scientist (I am a statistician) so I am not sure that I have a target convergence rate, I only hope to find an accurate approximation of the infinite sum.

  • 1
    $\begingroup$ A usage question: Roughly speaking, are you evaluating the sum for many values of $\lambda$ for a single $\nu$, many values of $\nu$ for a single $\lambda$, or something else? And what sort of accuracy are you looking for? $\endgroup$ Commented Oct 9, 2021 at 0:11
  • 2
    $\begingroup$ Can you provide some range for $\nu$ and $\lambda$, is it over all values or can you restrict yourself to some ranges? $\endgroup$
    – PC1
    Commented Oct 9, 2021 at 0:38
  • 1
    $\begingroup$ First of all, you need to know apriori that the infinite sum is convergent, and that's something that you cannot prove numerically. Next you need to know at what rate it is supposed to be converging, i.e., the residual $r=|\Sigma_1^N - \Sigma_1^{\infty}$| vs N, and that's also something that has to be done analytically. Then you can establish numerically that the residual is on the right asymptote, which would suggest that the series (the sum) is close to converging, and would let you estimate the limit. $\endgroup$ Commented Oct 9, 2021 at 4:28
  • 2
    $\begingroup$ You need to establish convergence theoretically. Having the convergence behavior / order at hand, would be a plus. You can also try to figure this out experimentally. Due to the recurrence relation of $I_\nu$, Clenshaw's Algorithm might work. Of course, you need to truncate the series. From there, you could use convergence acceleration / sequence transformation / extrapolation methods: Wynn's $\varepsilon$, Wynn's $\rho$ method, Richardson extrapolation, etc. $\endgroup$
    – cos_theta
    Commented Oct 11, 2021 at 7:56
  • $\begingroup$ Thanks for the comments. I've updated the question to reflect some of the issues raised. $\endgroup$
    – Lewkrr
    Commented Oct 11, 2021 at 18:21

1 Answer 1


Note the identity for the modified Bessel functions of the first kind,

$ e^z = I_0(z) + 2 \sum_{k=1}^{\infty} I_k(z) $

(Abramowitz and Stegun, Eq. 9.6.34, https://personal.math.ubc.ca/~cbm/aands/page_376.htm)

Using it, we can rewrite the infinite sum as

$ \sum_{i=1}^{\infty} I_{\nu+i}(2\lambda) = \sum_{i=1}^{\infty} I_{i}(2\lambda) - \sum_{i=1}^{\nu} I_{i}(2\lambda) = \frac{1}{2}\left(e^{2\lambda} - I_0(2\lambda)\right) - \sum_{i=1}^{\nu} I_{i}(2\lambda) $

Thus, evaluation of this expression does not require any considerations of convergence, one just needs to calculate the finite sum $\sum_{i=1}^{\nu}$.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.