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