# incompleteBesselK: The Incomplete Bessel K Function In DistributionUtils: Distribution Utilities

## Description

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

## Usage

 ```1 2 3 4 5 6``` ```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.

## 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)=integral_1^infinity 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.

## Value

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

## Note

The problem of calculation of the incomplete Bessel K function is equivalent to the problem of calculation of the cumulative distribution function of the generalized inverse Gaussian distribution. See ```Generalized Inverse Gaussian.```

## Author(s)

David Scott d.scott@auckland.ac.nz, Thomas Tran, Richard Slevinsky, Hassan Safouhi.

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

`besselK`

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32``` ```### 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) ```

### Example output

```Loading required package: RUnit
 2.225310761266483
 0.003246798003148282
 1.224998798576096e-05
 4.150045864731306e-07
 0.0005285043252121944
 0.003747507782386040 0.003170736451340495 0.002737234415162217
 0.002401939898349073 0.002136211905082957 0.001921197711255536
 0.001744093999367867 0.001595960459265232 0.001470397577616868
 0.001362722145744567
[,1]                  [,2]                  [,3]
value        0.003747507794965931  0.003170736456762924  0.002737234415977192
abs.error    2.632198252527153e-07 6.420325315237256e-07 1.648957840166878e-07
subdivisions 1                     1                     1
message      "OK"                  "OK"                  "OK"
call         Expression            Expression            Expression
[,4]                  [,5]                  [,6]
value        0.002401939898248203  0.002136211905014395  0.001921197711256391
abs.error    3.516601189747829e-10 1.919910717446044e-08 5.817213851785149e-09
subdivisions 1                     1                     1
message      "OK"                  "OK"                  "OK"
call         Expression            Expression            Expression
[,7]                  [,8]                  [,9]
value        0.001744093999384008  0.001595960459294433  0.001470397577636734
abs.error    5.074635570395979e-10 1.535913138303264e-09 1.253679558192617e-10
subdivisions 1                     1                     1
message      "OK"                  "OK"                  "OK"
call         Expression            Expression            Expression
[,10]
value        0.001362722145757839
abs.error    3.704750934196496e-10
subdivisions 1
message      "OK"
call         Expression
 -1.257989112185554e-11 -5.422428998869133e-12 -8.149756910991179e-13
  1.008698333193614e-13  6.856234330276934e-14 -8.554355140910630e-16
 -1.614051774179170e-14 -2.920081711155209e-14 -1.986627008732089e-14
 -1.327171879339506e-14
 1.257989112185554e-11
```

DistributionUtils documentation built on May 1, 2019, 9:12 p.m.