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

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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.