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$ – Noam D. Elkies Nov 13 '18 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$ – AccidentalFourierTransform Nov 13 '18 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$ – Stefano Vespucci Nov 13 '18 at 3:33
  • $\begingroup$ What do you mean by "high input positive integer values"? $\endgroup$ – Amir Sagiv Nov 13 '18 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$ – 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.

  • $\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$ – Stefano Vespucci Nov 15 '18 at 2:10
  • $\begingroup$ I would up-vote your answer! But I am not allowed :( $\endgroup$ – 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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.