[R] Bessel function with large index value

[Martin Maechler]

>>>>> "DS" == David Scott <d.scott at auckland.ac.nz>
>>>>>     on Sat, 21 Nov 2009 02:29:38 +1300 writes:

    DS> This is a reply to my own question. I thought I had found an answer but 
    DS> it seems not so (some analysis follows below). Maybe Martin Maechler or 
    DS> Robin Hankin or Duncan Murdoch may have some ideas---I know the question 
    DS> is a bit specialized.

Indeed.

The good news is that last winter (when in the Swiss Alps, typically
after enjoying the sun and snow) I've written an experimental and
unpublished R package called  'Bessel'
where I've started to collect bessel implementations /
approximations that I had collected and partly tested in the
past,  also looking at what other packages on CRAN had.
My 'Bessel' package depends on my new package  'Rmpfr'
(for arbitrary precision arithmetic), as I really want to
explore the quality of the different Bessel computations, and
for that  have wanted access to "infinite" precision arithmetic.

Apropos 'Rmpfr': This has been on R-forge for about a year, and
	now become a CRAN package a couple of weeks ago.
 Thanks to Brian Ripley and Stefan Theussl, the package is
 available for Windows from R-forge, but not yet from CRAN.
 For MacOSX, it has not yet been made available in precompiled
 form, but there you should be able to build it from the
 sources, after installing the (GNU) GMP library, and the MPFR
 library which builds on GMP.

Back to Bessel and my experimental code:
In there, I have a function, currently only for the Bessel
I_[\nu] that uses the Debye polynomial series needed for these
cases,
and indeed my function has a  'log' argument.

    DS> 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.

As my implementation uses A_&_S 9.3.7 for I_nu,
I wonder why you say that that  the completely analogous 9.3.8
for  K_nu is not appropriate.

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

so, working on the log scale saves "all problems", right?

I think we should further correspond offline on the topic;
given your interest, I guess I should probably make my 'Bessel' available
from R-forge ....

Martin Maechler, ETH Zurich

    >> David Scott
    >> 
    >> 

    DS> After posting I checked the GNU Scientific Library 
    DS> (http://www.gnu.org/software/gsl/) and found:

    DS> ********************
    DS> — Function: double gsl_sf_bessel_lnKnu (double nu, double x)
    DS> — Function: int gsl_sf_bessel_lnKnu_e (double nu, double x, 
    DS> gsl_sf_result * result)

    DS> These routines compute the logarithm of the irregular modified 
    DS> Bessel function of fractional order \nu, \ln(K_\nu(x)) for x>0, \nu>0.
    DS> ********************
    DS> I then recalled that Robin Hankin and Duncan Murdoch had made the GSL 
    DS> available. I installed the package gsl and investigated the function
    DS> bessel_lnKnu.

    DS> Unfortunately, it appears that this function has *smaller* range than 
    DS> besselK when it comes to the index. The following shows it:

    DS> library(plyr)
    DS> library(gsl)
    DS> ### Check calculations using both methods
    DS> lnKnu <- maply(expand.grid(x = 100*(1:7), nu = 10*(1:100)), bessel_lnKnu)
    DS> lnKnu
    DS> Knu <- maply(expand.grid(x = 100*(1:7), nu = 10*(1:100)), besselK)
    DS> Knu
    DS> lnKnu/log(Knu)

    DS> I was expecting what happens with gamma and lgamma

    DS> ### Compare gamma function
    DS> lgam <- lgamma(100*(1:7))
    DS> lgam
    DS> gam <- gamma(100*(1:7))
    DS> gam
    DS> lgam/log(gam)

    DS> It seems that bessel_lnKnu is set up to protect against infinity when x 
    DS> becomes small:

    DS> ### Does lnKnu protect against Inf when x goes to zero?
    DS> lnnear0 <- maply(expand.grid(x = 0.00000001*(1:10), nu = 10*(0:5)), 
    DS> bessel_lnKnu)
    DS> lnnear0
    DS> near0 <- maply(expand.grid(x = 0.00000001*(1:10), nu = 10*(0:5)), besselK)
    DS> near0
    DS> lnnear0/log(near0)

    DS> So, I am still in need of a solution: an implementation of log of Bessel 
    DS> K which protects against the index nu becoming large.

    DS> David Scott

    DS> -- 
    DS> _________________________________________________________________
    DS> David Scott	Department of Statistics
    DS> The University of Auckland, PB 92019
    DS> Auckland 1142,    NEW ZEALAND
    DS> Phone: +64 9 923 5055, or +64 9 373 7599 ext 85055
    DS> Email:	d.scott at auckland.ac.nz,  Fax: +64 9 373 7018

    DS> Director of Consulting, Department of Statistics

    DS> ______________________________________________
    DS> R-help at r-project.org mailing list
    DS> https://stat.ethz.ch/mailman/listinfo/r-help
    DS> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
    DS> and provide commented, minimal, self-contained, reproducible code.