DistributionUtils (version 0.6-0)

incompleteBesselK: The Incomplete Bessel K Function

Description

Calculates the incomplete Bessel K function using the algorithm and code provided by Slavinsky and Safouhi (2009).

Usage

incompleteBesselK(x, y, nu, tol = (.Machine$double.eps)^(0.85), nmax = 120)
incompleteBesselKR(x, y, nu, tol = (.Machine$double.eps)^(0.85), nmax = 120)
SSFcoef(nmax, nu)
combinatorial(nu)
GDENOM(n, x, y, nu, An, nmax, Cnp)
GNUM(n, x, y, nu, Am, An, nmax, Cnp, GM, GN)

Arguments

x

Numeric. Value of the first argument of the incomplete Bessel K function.

y

Numeric. Value of the second argument of the incomplete Bessel K function.

nu

Numeric. The order of the incomplete Bessel K function.

tol

Numeric. The tolerance for the difference between successive approximations of the incomplete Bessel K function.

nmax

Integer. The maximum order allowed for the approximation of the incomplete Bessel K function.

n

Integer. Current order of the approximation. Not required to be specified by users.

An

Matrix of coefficients. Not required to be specified by users.

Am

Matrix of coefficients. Not required to be specified by users.

Cnp

Vector of elements of Pascal's triangle. Not required to be specified by users.

GN

Vector of denominators used for approximation. Not required to be specified by users.

GM

Vector of numerators used for approximation. Not required to be specified by users.

Value

incompleteBesselK and incompleteBesselKR both return an approximation to the incomplete Bessel K function as defined above.

Details

The function incompleteBesselK implements the algorithm proposed by Slavinsky and Safouhi (2010) and uses code provided by them.

The incomplete Bessel K function is defined by

$$K_\nu(x,y)=\int_1^\infty t^{-nu-1}\exp(-xt-y/t)\,dt$$

see Slavinsky and Safouhi (2010), or Harris (2008).

incompleteBesselK calls a Fortran routine to carry out the calculations. incompleteBesselKR is a pure R version of the routine for computing the incomplete Bessel K function.

The functions SSFcoef, combinatorial, GDENOM, and GNUM are subroutines used in the function incompleteBesselKR. They are not expected to be called by the user and the user is not required to specify input values for these functions.

The approximation to the incomplete Bessel K function returned by incompleteBesselK is highly accurate. The default value of tol is about 10^(-14) on a 32-bit computer. It appears that even higher accuracy is possible when x > y. Then the tolerance can be taken as .Machine$double.eps and the number of correct figures essentially coincides with the number of figures representable in the machine being used.

incompleteBesselKR is very slow compared to the Fortran version and is only included for those who wish to see the algorithm in R rather than Fortran.

References

Harris, Frank E. (2008) Incomplete Bessel, generalized incomplete gamma, or leaky aquifer functions. J. Comp. Appl. Math., 215, 260--269.

Slevinsky, Richard M., and Safouhi, Hassan (2009) New formulae for higher order derivatives and applications. J. Comp. Appl. Math. 233, 405--419.

Slevinsky, Richard M., and Safouhi, Hassan (2010) A recursive algorithm for the G transformation and accurate computation of incomplete Bessel functions. Appl. Numer. Math., In press.

See Also

besselK

Examples

# NOT RUN {
### Harris (2008) gives accurate values (16 figures) for
### x = 0.01, y = 4, and nu = 0:9
### nu = 0, Harris value is 2.22531 07612 66469
options(digits = 16)
incompleteBesselK(0.01, 4, 0)
### nu = 9, Harris value is 0.00324 67980 03149
incompleteBesselK(0.01, 4, 9)

### Other values given in Harris (2008)
### x = 4.95, y = 5.00, nu = 2
incompleteBesselK(4.95, 5, 2) ## 0.00001 22499 87981
### x = 10, y = 2, nu = 6
### Slevinsky and Safouhi (2010) suggest Harris (2008) value
### is incorrect, give value 0.00000 04150 01064 21228
incompleteBesselK(10, 2, 6)
### x = 3.1, y = 2.6, nu = 5
incompleteBesselK(3.1, 2.6, 5) ## 0.00052 85043 25244

### Check values when x > y using numeric integration
(numIBF <- sapply(0:9, incompleteBesselK, x = 4, y = 0.01))

besselFn <- function(t, x, y, nu) {
  (t^(-nu - 1))*exp(-x*t - y/t)
}

(intIBF <- sapply(0:9, integrate, f = besselFn, lower = 1, upper = Inf,
                 x = 4, y = 0.01))
intIBF <- as.numeric(intIBF[1, ])
numIBF - intIBF
max(abs(numIBF - intIBF)) ## 1.256649992398273e-11

options(digits = 7)
# }