The besselK() function knows about overflows/underflows internally;
there is a constant xmax_BESS_K in src/nmath/bessel.h (and referred to
only in bessel_k.c), equal to 705.342, which is checked if expon.scaled is
FALSE. (The equivalent number for bessel_i.c is 709, defined as
exparg_BESS in bessel.h.) However, besselK(x) silently returns +Inf if
x>705.342. This behavior is reasonable for besselI (where the problem is
an overflow), but seems wrong for besselK (since the problem here is an
*underflow*). The answer may be as simple as returning 0 rather than +Inf
in this case.
If this is appropriate, the following patch should do the job:
*** bessel_k.c Thu Oct 17 13:55:55 2002
--- bessel_k.c.dist Thu Oct 17 13:48:53 2002
***************
*** 216,224 ****
ML_ERROR(ME_RANGE);
*ncalc = *nb;
for(i=0; i < *nb; i++)
! if (ex<=0)
! bk[i] = ML_POSINF;
! else bk[i]=0.0;
return;
}
k = 0;
--- 216,222 ----
ML_ERROR(ME_RANGE);
*ncalc = *nb;
for(i=0; i < *nb; i++)
! bk[i] = ML_POSINF;
return;
}
k = 0;
(By the way, the test for x<=0 seems inappropriate; it looks like x<0 is
caught elsewhere and that this is really a test for x==0.)
[Apologies if this is not considered a bug, but my judgment is that it
is -- couldn't be bothered with a round of sending it to r-devel first.]
Ben Bolker
--please do not edit the information below--
Version:
platform = i686-pc-linux-gnu
arch = i686
os = linux-gnu
system = i686, linux-gnu
status =
major = 1
minor = 6.0
year = 2002
month = 10
day = 01
language = R
Search Path:
.GlobalEnv, package:ctest, Autoloads, package:base
--
318 Carr Hall bolker@zoo.ufl.edu
Zoology Department, University of Florida http://www.zoo.ufl.edu/bolker
Box 118525 (ph) 352-392-5697
Gainesville, FL 32611-8525 (fax) 352-392-3704
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To:
r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._