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

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)
``` |

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

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.

`incompleteBesselK`

and `incompleteBesselKR`

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

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

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

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`

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)
``` |

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.