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.

enter image description here


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.

  • $\begingroup$ Maybe this will help? carma.newcastle.edu.au/jon/jmm.pdf $\endgroup$ Nov 13, 2018 at 3:15
  • $\begingroup$ Do you want a large-$x$ expansion, or a large-$\alpha$ expansion? In any case, see Spherical Bessel functions and their asymptotics. $\endgroup$ Nov 13, 2018 at 3:17
  • $\begingroup$ 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 $\endgroup$ Nov 13, 2018 at 3:33
  • $\begingroup$ What do you mean by "high input positive integer values"? $\endgroup$
    – Amir Sagiv
    Nov 13, 2018 at 3:47
  • $\begingroup$ @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 $\endgroup$ Nov 13, 2018 at 3:58

3 Answers 3


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.

  • $\begingroup$ 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). $\endgroup$ Nov 15, 2018 at 2:10
  • $\begingroup$ I would up-vote your answer! But I am not allowed :( $\endgroup$ Nov 21, 2018 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.


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.