[R] Bessel function with large index value

[David Scott]

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