This is a reply to my own question. I thought I had found an answer but
it seems not so (some analysis follows below). Maybe Martin Maechler or
Robin Hankin or Duncan Murdoch may have some ideas---I know the question
is a bit specialized.
David Scott wrote:> I am looking for a method of dealing with the modified Bessel function
> K_\nu(x) for large \nu.
>
> The besselK function implementation of this allows for dealing with
> large values of x by allowing for exponential scaling, but there is no
> facility for dealing with large \nu.
>
> What would work for me would be an lbesselK function in the manner of
> lgamma which returned the log of K_\nu(x) for large \nu. Does anybody
> have any leads on this?
>
> Note that I have trawled through Abramowitz and Stegun and found 9.7.8
> which doesn't work for me because of the complication in the definition
> of the x argument. I have also seen a result of Ismail (1977) reported
> by Barndorff-Nielsen and Blaesild which has the other problem, the
> treatment of the x argument is too simple.
>
> To do the calculation I am attempting, I need to have the Bessel
> function in a form that will allow a cancellation with a Gamma function
> of high order in the denominator.
>
> David Scott
>
>
After posting I checked the GNU Scientific Library
(http://www.gnu.org/software/gsl/) and found:
********************
? Function: double gsl_sf_bessel_lnKnu (double nu, double x)
? Function: int gsl_sf_bessel_lnKnu_e (double nu, double x,
gsl_sf_result * result)
These routines compute the logarithm of the irregular modified
Bessel function of fractional order \nu, \ln(K_\nu(x)) for x>0, \nu>0.
********************
I then recalled that Robin Hankin and Duncan Murdoch had made the GSL
available. I installed the package gsl and investigated the function
bessel_lnKnu.
Unfortunately, it appears that this function has *smaller* range than
besselK when it comes to the index. The following shows it:
library(plyr)
library(gsl)
### Check calculations using both methods
lnKnu <- maply(expand.grid(x = 100*(1:7), nu = 10*(1:100)), bessel_lnKnu)
lnKnu
Knu <- maply(expand.grid(x = 100*(1:7), nu = 10*(1:100)), besselK)
Knu
lnKnu/log(Knu)
I was expecting what happens with gamma and lgamma
### Compare gamma function
lgam <- lgamma(100*(1:7))
lgam
gam <- gamma(100*(1:7))
gam
lgam/log(gam)
It seems that bessel_lnKnu is set up to protect against infinity when x
becomes small:
### Does lnKnu protect against Inf when x goes to zero?
lnnear0 <- maply(expand.grid(x = 0.00000001*(1:10), nu = 10*(0:5)),
bessel_lnKnu)
lnnear0
near0 <- maply(expand.grid(x = 0.00000001*(1:10), nu = 10*(0:5)), besselK)
near0
lnnear0/log(near0)
So, I am still in need of a solution: an implementation of log of Bessel
K which protects against the index nu becoming large.
David Scott
--
_________________________________________________________________
David Scott Department of Statistics
The University of Auckland, PB 92019
Auckland 1142, NEW ZEALAND
Phone: +64 9 923 5055, or +64 9 373 7599 ext 85055
Email: d.scott at auckland.ac.nz, Fax: +64 9 373 7018
Director of Consulting, Department of Statistics