General approach to infinite sums

Question

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\}$$

so

$$ 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.

2. 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$

3. 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.

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}$.