Bessel (version 0.6-1)

bI: Bessel I() function Simple Series Representation

Description

Computes the modified Bessel \(I\) function, using one of its basic definitions as an infinite series. The implementation is pure R, working for numeric, complex, but also e.g., for objects of class "mpfr" from package Rmpfr.

Usage

besselIs(x, nu, nterm = 800, expon.scaled = FALSE, log = FALSE,
         Ceps = if (isNum) 8e-16 else 2^(-x@.Data[[1]]@prec))

Value

a “numeric” (or complex or "mpfr") vector of the same class and length as x.

Arguments

x

numeric or complex vector, or of another class for which arithmetic methods are defined, notably objects of class mpfr (package Rmpfr).

nu

non-negative numeric (scalar).

nterm

integer indicating the number of terms to be used. Should be in the order of abs(x), but can be smaller for large x. A warning is given, when nterm was chosen too small.

expon.scaled

logical indicating if the result should be scaled by \(exp(-abs(x))\).

log

logical indicating if the logarithm \(log I.()\) is required. This allows even more precision than expon.scaled=TRUE in some cases.

Ceps

a relative error tolerance for checking if nterm has been sufficient. The default is “correct” for double precision and also for multiprecision objects.

Author

Martin Maechler

References

Abramowitz, M., and Stegun, I. A. (1964,.., 1972). Handbook of mathematical functions (NBS AMS series 55, U.S. Dept. of Commerce).

See Also

This package BesselI, base besselI, etc

Examples

(nus <- c(outer((0:3)/4, 1:5, `+`)))
stopifnot(
  all.equal(besselIs(1:10, 1), # our R code
            besselI (1:10, 1)) # internal C code w/ different algorithm
  ,
  sapply(nus, function(nu)
   all.equal(besselIs(1:10, nu, expon.scale=TRUE), # our R code
             BesselI (1:10, nu, expon.scale=TRUE)) # TOMS644 code
   )
  ,
  ## complex argument [gives warnings  'nterm=800' may be too small]
  sapply(nus, function(nu)
   all.equal(besselIs((1:10)*(1+1i), nu, expon.scale=TRUE), # our R code
             BesselI ((1:10)*(1+1i), nu, expon.scale=TRUE)) # TOMS644 code
   )
)

## Large 'nu' ...
x <- (0:20)/4
(bx <- besselI(x, nu=200))# base R's -- gives (mostly wrong) warnings
if(require("Rmpfr")) { ## Use high precision (notably large exponent range) numbers:
  Bx <- besselIs(mpfr(x, 64), nu=200)
  all.equal(Bx, bx, tol = 1e-15)# TRUE -- warning were mostly wrong; specifically:
  cbind(bx, Bx)
  signif(asNumeric(1 - (bx/Bx)[19:21]), 4) # only [19] had lost accuracy

  ## With*out* mpfr numbers -- using log -- is accurate (here)
  (lbx <- besselIs(     x,      nu=200, log=TRUE))
  lBx <-  besselIs(mpfr(x, 64), nu=200, log=TRUE)
  stopifnot(all.equal(asNumeric(log(Bx)), lbx, tol=1e-15),
	    all.equal(lBx, lbx, tol=4e-16))
} # Rmpfr