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.
 A: 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. 
A: 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.
A: 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.

