I'm attempting to find a precise as possible approximation to the logarithm of the Modified Bessel Function of the Second Kind:
$$\log K_{\alpha}(x) = \log\Big[\frac{1}{2}\int_0^{∞} t^{\alpha-1} \exp\{-\frac{1}{2}x(t+t^{-1})\}dt\Big].$$
The problem is that $K_{\alpha}(x)$ diverges both as $\alpha\rightarrow\infty$ and as $x\rightarrow0$, which is why I'm taking the logarithm in the first place. But of course the log trick only works if I can use it to break up the expression in some way, which I can't do in this case because I can't pass it through the integral.
There is an asymptotic approximation for large $\alpha$, given by:
$$K_{\alpha}(x)\sim\sqrt{\frac{\pi}{2\alpha}}\Big(\frac{ex}{2\alpha}\Big)^{-\alpha}.$$
Clearly this approximation can be broken up by the log, unfortunately it is not sufficiently precise (for my purposes) for values of $\alpha$ and $x$ which cause overflow.
Does anyone have another way to get a more precise estimate for $\log K_{\alpha}(x)$?
If you want you can assume that $\alpha=k-\frac{1}{2}$ for $k\in{\bf Z}_+$.
Best Answer
I am currently struggling with the same problem, getting a bit further but not all the way. If you solved it in the meantime, I'd be glad to hear about it.
EDIT 2: Improved answer some more as my technique improved some more.
EDIT 3: The pull request (linked below) now has somewhat better implementation then discussed here that reduces the error a bit further, but still breaks for some values.
C++ code (written for Stan math library) can be found at https://github.com/stan-dev/math/pull/1121
Case 0: $x$ small relative to $\alpha$
I use this when $\alpha > 15$ AND $10x < \sqrt{\alpha + 1}$. The formula can be found on Wikipedia on Bessel K and also as formula 1.10 of Temme, Journal of Computational Physics, vol 19, 324 (1975). It has:
$$ K_\alpha(x) \simeq \frac{\Gamma(\alpha)}{2} \left( \frac{2}{x} \right)^\alpha $$
Case 1: Small to medium $\alpha$ and $x$
My main workhorse is using the logarithm of Equation 26 of Rothwell: Computation of the logarithm of Bessel functions of complex argument and fractional order.
which has (maintaining their notation)
$$ K_\nu(z) = \frac{\sqrt{\pi}}{\Gamma(\nu + \frac{1}{2})}(2z)^{-\nu}e^{-z} \int_0^1 [\beta e^{-u^\beta}(2z + u^\beta)^{\nu-1/2} u^{n-1} + e^{-1/u} u^{-2\nu-1}(2zu + 1)^{\nu-1/2}] \mathrm{d}u $$
Where $\beta = \frac{2n}{2\nu+1}$ and the authors suggest using $n = 8$.
You still can't push the log through the integral, but for large areas of the parameter space, it is OK to just compute the integral and take the log afterwards.
The integral ceases to be numerically tractable around roughly $\nu > 50$ or when $\nu > 0.5$ AND $log(z) > \frac{300}{\nu - 0.5} - log(2)$ (one of the summands turn to infinity).
Fruther, the formula ceases to be very accurate for small $z$ and small $\nu$, where below roughly $z < 10^-4$ the relative error starts climbing from around $10^{-30}$ to around $10{^-2}$ for $z < 10^{-7}$
There is a C++ implementation of the Rothwell formula at https://github.com/stan-dev/stan/wiki/Stan-Development-Meeting-Agenda/0ca4e1be9f7fc800658bfbd97331e800a4f50011
Case 2: $\alpha \gg x$
Here the log of the approximate formula you've written works OK, but I got a bit better results with the same formula as in Case 0.
Case 3: $x \gg \alpha$
Asymptotic formula such as 10.4.2 here: https://dlmf.nist.gov/10.40 can be used. Once again, you need to compute the sum on non-log scale (it has negative elements), and take the log afterward, but this has worked great for me. Summing up to 10 first elements worked OK in my case.
Remaining cases For large $x$ with comparably large $\alpha$, I still haven't found a reliable solution.
Here is a plot of error with respect to recursive formula for consecutive neighbours with the approach I've described (
ratiois the relative error of the actual value, I also test the gradient of this function wrt. both parameters as computed by Stan's autodiff - the gradient is a tougher problem, but probably not relevant here):The biggest relative error shown is 3.8e-02 for $\alpha = 148$ and $x = 105$.