# General approach to infinite sums

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.

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.

• 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? Oct 9 '21 at 0:11
• Can you provide some range for $\nu$ and $\lambda$, is it over all values or can you restrict yourself to some ranges?
– PC1
Oct 9 '21 at 0:38
• 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. Oct 9 '21 at 4:28
• 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. Oct 11 '21 at 7:56
• Thanks for the comments. I've updated the question to reflect some of the issues raised. Oct 11 '21 at 18:21

$$e^z = I_0(z) + 2 \sum_{k=1}^{\infty} I_k(z)$$
$$\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}$$.