>>>>> "DS" == David Scott <d.sc...@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.sc...@auckland.ac.nz, Fax: +64 9 373 7018 DS> Director of Consulting, Department of Statistics DS> ______________________________________________ DS> R-help@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. ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.