R help - Nov 2009 - Bessel function with large index value

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


-- 
_________________________________________________________________
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

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