# Approximation of half-integers modified Bessel function of the second kind

I am trying to optimise the calculation of the probability distribution poisson-inverse-gaussian

its calculation involves a half-integers modified Bessel function of the second kind. Here's a formula for the probabilistic expression. Pubblication

For high input positive integer values, the calculation is very slow, as it involves a summation, or a recursion avoiding the Bessel in another formulation.

Is there an approximation for

$$K_{x-1/2}(\alpha); x \gg 0$$

Thanks a lot.

• Maybe this will help? carma.newcastle.edu.au/jon/jmm.pdf – Noam D. Elkies Nov 13 '18 at 3:15
• Do you want a large-$x$ expansion, or a large-$\alpha$ expansion? In any case, see Spherical Bessel functions and their asymptotics. – AccidentalFourierTransform Nov 13 '18 at 3:17
• Thanks. (I will use a probabilistic MCMC algorithm (Stan)) data counts $x$ can be high (e.g., 100.000) and the parameter $\alpha$ is "shape" I suppose is not >> 1 – Stefano Vespucci Nov 13 '18 at 3:33
• What do you mean by "high input positive integer values"? – Amir Sagiv Nov 13 '18 at 3:47
• @AmirSagiv for poisson_inverse_gaussian x = observed counts. So I am analysing gene expression where there could easily be and average of 100K copies for a given gene – Stefano Vespucci Nov 13 '18 at 3:58

Using a simple saddle point approximation one can derive the following: $$[A] \quad K_y(a) \sim \frac{1}{2} \exp{(-r + y\,s)}\sqrt{\frac{\pi}{2r}}\Big( 1 + \text{erf}\big(\sqrt{\frac{r}{2}}\, s\big) \Big) \quad$$ where $$r=\sqrt{a^2+y^2} \text{ and } s=\text{arcsinh}(y/a)$$ I derived it from the known integral relation $$K_y(a) = \int_0^\infty \exp(-a\,\cosh(t))\cosh(y\,t) dt$$ Approximate $$\cosh(yt) \sim 1/2 \exp{(yt)} .$$ The integrand sans the 1/2 is $$\exp(-a\,\cosh(t) + yt).$$ Expand the argument of the exponential around its saddle point $$s = \text{arcsinh}(y/a)$$ and you get $$K_y(a) \sim \frac{1}{2} \exp{(-r + y\,s)} \int_0^\infty \exp{(-\frac{r}{2} (t-s)^2)}\, dt$$ The integral evaluates to an error function, which is shown in [A]. Alternatively, if the argument of the error function is sufficiently large, say, >5, the the expression in the big parentheses can be approximated by 1.

Using Mathematica, for $$y=x-1/2=10000,$$ and $$a=1/2$$ I get about 6 significant figures agreement. For $$a=11/2,$$ I get about 5 sig figs. It is expected that this approximation will deteriorate as $$a$$ becomes comparable to $$x$$ but I haven't studied the limits. Comments seem to indicate that $$x>>a$$ so this formula may be sufficient. One can also get more terms by doing a complete expansion.

• Thanks a lot @skbmoore. This error trend is also true if I truncate the summation in besselK to the interval [y-10, y]. I am currently trying to understand whether the shape parameter alpha can get of the same size as x (e.g., alpha = 100-1000, x = 1000). – Stefano Vespucci Nov 15 '18 at 2:10
• I would up-vote your answer! But I am not allowed :( – Stefano Vespucci Nov 21 '18 at 22:52

As I read Figure 2 in "Algorithm 644: A portable package for Bessel functions of a complex argument and nonnegative order" (D.E.Amos, https://dl.acm.org/citation.cfm?id=214331) you're veering dangerously close into 'overflow' territory.

Some experimentation with scipy.special.kv bears this out. Computing $$K_{1000.5}(x)$$ for $$600\le x < 700$$ gives values ranging from $$K_{1000.5}(600) \approx 4.36527371e+49$$ to $$K_{1000.5}(700) \approx 3.73027792e-30$$. (and plotting the logs produces a nearly-straight line, which suggests that this is in fact in the asymptotic region).

Looking at your original equation, I see an $$x!$$ in a denominator, and a $$(\cdot)^x$$ in a numerator, which aren't exactly going to make the numerics easier. What I would recommend is switching to logspace and using the large-order approximation ( (10.41.2) on DLMF). You'll have to tweak the MCMC update rules, but not in any particularly tricky way.

You might also be interested in this math.se question : https://math.stackexchange.com/questions/1960778/approximating-the-log-of-the-modified-bessel-function-of-the-second-kind and some of the links therein.

I was not able to access this paper (paywall), but judging from the abstract it should provide what you are looking for: Computation of the Poisson-inverse Gaussian distribution

Recursion relations suitable for rapid computation are derived for the probabilities of the compound Poisson distribution when the compounder is the inverse-Gaussian distribution. Series representation of the probabilities are given. Asymptotic results as well as approximations for probabilities, compared with the exact values, are investigated.