Nothing
R/I-fn.R
In Bessel: Computations and Approximations for Bessel Functions
Defines functions
besselIs
bI
besselIasym
besselI.ftrms
besselI.nuAsym
besselK.nuAsym
Documented in
besselIasym
besselI.ftrms
besselI.nuAsym
besselIs
besselK.nuAsym
bI
###--- This used to be part /u/maechler/R/MM/NUMERICS/bessel-fn.R
#### "Modified" Bessel Function I_nu(x)
#### ----------------------------------------------
## besselI() - Definition as infinite sum -- working for "mpfr" numbers, too!
## ---------
besselIs <-
function(x, nu, nterm = 800, expon.scaled = FALSE, log = FALSE,
Ceps = if(isNum) 8e-16 else 2^(- x@.Data[[1]]@prec))
{
## Purpose: besselI() primitively
## ----------------------------------------------------------------------
## Arguments: (x,nu) as besselI; nterm: number of terms for "infinite" sum
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 21 Oct 2002, 21:59
if(length(nu) > 1)
stop(" 'nu' must be scalar (length 1)!")
j <- (nterm-1):0 # sum smallest first!
n <- length(x)
if(n == 0) return(x)
if(is(nu, "mpfr")) x <- mpfr(x, precBits = max(64, .getPrec(nu)))
l.s.j <- outer(j, (x/2), function(X,Y) X*2*log(Y))##-> {nterm x n} matrix
##
isNum <- is.numeric(x) || is.complex(x)
## improve accuracy for lgamma(j+1) for "mpfr" numbers
## -- this is very important [evidence: e.g. bI(10000, 1)]
if(is(l.s.j, "mpfr"))
j <- mpfr(j, precBits = max(sapply(l.s.j@.Data, slot, "prec")))
else if(!isNum) j <- as(j, class(x))
## underflow (64-bit AMD) for x > 745.1332
## for large x, this already overflows to Inf :
## s.j <- outer(j, (x^2/4), function(X,Y) Y^X) ##-> {nterm x n} matrix
## s.j <- s.j / (gamma(j+1) * gamma(nu+1 + j)) but without overflow :
log.s.j <-
if(expon.scaled)
l.s.j - rep(x, each=nterm) - lgamma(j+1) - lgamma(nu+1 + j)
else
l.s.j - lgamma(j+1) - lgamma(nu+1 + j)
s.j <-
if(log) # NB: lsum() works on whole matrix
lsum(log.s.j) # == log(sum_{j} exp(log.s.j) )
else ## log I_nu(x) -- trying to avoid overflow/underflow for large x OR large nu
## log(s.j) ; e..x <- exp(-x) # subnormal for x > 1024*log(2) ~ 710;
exp(log.s.j)
if(log) {
if(any(i0 <- x == 0)) s.j[i0] <- 0
if(any(lrgS <- log.s.j[1,] > log(Ceps) + s.j))
lapply(x[lrgS], function(x)
warning(sprintf("bI(x=%g): 'nterm' may be too small", x),
call.=FALSE))
nu*log(x/2) + s.j
} else { ## !log
s <- colSums(s.j)
if(any(i0 <- x == 0)) s[i0] <- 0
if(!all(iFin <- is.finite(s)))
stop(sprintf("infinite s for x=%g", x[!iFin][1]))
if(any(lrgS <- s.j[1,] > Ceps * s))
lapply(x[lrgS], function(x)
warning(sprintf("bI(x=%g): 'nterm' may be too small", x),
call.=FALSE))
(x/2)^nu * s
}
}
## old name [back compatibility]:
bI <- function(x, nu, nterm = 800, expon.scaled = FALSE, log = FALSE,
Ceps = if(isNum) 8e-16 else 2^(- x@.Data[[1]]@prec))
{
.Deprecated("besselIs")
isNum <- is.numeric(x) || is.complex(x)
besselIs(x, nu, nterm=nterm, expon.scaled=expon.scaled, log=log, Ceps=Ceps)
}
###--------------- besselI() for large x or also large nu ---
###
### see --> ~/R/MM/NUMERICS/bessel-large-x.R for code usage
### ================================
besselIasym <- function(x, nu, k.max = 10, expon.scaled=FALSE, log=FALSE)
{
## Purpose: Asymptotic expansion of Bessel I_nu(x) function x -> oo
## by Abramowitz & Stegun (9.7.1), p.377 :
##
## I_a(z) = exp(z) / sqrt(2*pi*z) * f(z,..) where
## f(z,..) = 1 - (mu-1)/ (8*z) + (mu-1)(mu-9)/(2! (8z)^2) - ...
## = 1- (mu-1)/(8z)*(1- (mu-9)/(2(8z))*(1-(mu-25)/(3(8z))*..))
## where mu = 4*a^2 *and* |arg(z)| < pi/2
## ----------------------------------------------------------------------
## Arguments: x, nu, expon.scaled: as besselI()
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 11 Apr, 21 Nov 2008
## Note that for each x the series eventually *DIVERGES*
## it should be "stopped" as soon as it *converged* (:-)
stopifnot(k.max == round(k.max))
## I_\nu(z) = exp(z) / sqrt(2*pi*z) * f(z, nu)
## First compute f(x, nu) to order k.max
## f <- 1
d <- 0 ## d = 1 - f <==> f = 1 - d
if(k.max >= 1) {
## m. <- 4*nu^2
x8 <- 8*x
for(k in k.max:1) {
## mu <- 4*nu^2; d <- (1 - d)*(mu - (2*k-1)^2)/(k*x8)
## mu - (2k-1)^2 = (2nu - (2k-1)) (2nu + (2k-1)) =
## = (2(nu-k)+1)(2(nu+k)-1)
d <- (1 - d)*((2*(nu-k)+1)*(2*(nu+k)-1))/(k*x8)
}
}
if(log) {
## f = 1 - d ==> log(f) = log1p(-d) :
(if(expon.scaled) log1p(-d) else x + log1p(-d)) - log(2*pi*x) / 2
} else {
(if(expon.scaled) (1-d) else exp(x)*(1-d)) / sqrt(2*pi*x)
}
}
besselI.ftrms <- function(x, nu, K = 20)
{
## Purpose: all *Terms* in besselIasym()
## by Abramowitz & Stegun (9.7.1), p.377 :
## f(x,..) = 1 - (mu-1)/ (8*x) + (mu-1)(mu-9)/(2! (8x)^2) - ...
## = 1- (mu-1)/(8x)*(1- (mu-9)/(2(8x))*(1-(mu-25)/(3(8x))*..))
## where mu = 4*nu^2
## ----------------------------------------------------------------------
## Arguments: x, nu, expon.scaled: as besselI()
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 11 Apr, 21 Nov 2008
stopifnot(length(x) == 1, length(nu) == 1, length(K) == 1,
x > 0, K == round(K), K >= 0)
## I_\nu(x) = exp(x) / sqrt(2*pi*x) * f(x, nu)
kk <- seq_len(K)
mu <- 4*nu^2
x8 <- 8*x
multf <- - (mu - (2*kk-1)^2)/(kk*x8)
ser <- cumprod(multf)
## now, for k-term approx. f_k of f(x,nu) we have
## f_k = 1 + sum(ser[1:k]), i.e. ## (f_k)_k = 1 + cumsum(c(0, ser))[k+1] for k= 0,1,...,K
ser
}
###----------------------------------------------------------------------------
### When BOTH nu and x are large :
besselI.nuAsym <- function(x, nu, k.max, expon.scaled=FALSE, log=FALSE)
{
## Purpose: Asymptotic expansion of Bessel I_nu(x) function
## when BOTH nu and x are large
## Abramowitz & Stegun , p.378, __ 9.7.7. __
##
## I_nu(nu * z) ~ 1/sqrt(2*pi*nu) * exp(nu*eta)/(1+z^2)^(1/4) *
## * {1 + u_1(t)/nu + u_2(t)/nu^2 + ... }
## where t = 1 / sqrt(1 + z^2),
## eta = sqrt(1 + z^2) + log(z / (1 + sqrt(1+z^2)))
##
##
## and u_k(t) from p.366 __ 9.3.9 __
## u0(t) = 1
## u1(t) = (3*t - 5*t^3)/24
## u2(t) = (81*t^2 - 462*t^4 + 385*t^6)/1152
## ...
## with recursion 9.3.10 for k = 0, 1, .... :
##
## u_{k+1}(t) = t^2/2 * (1 - t^2) * u'_k(t) +
## 1/8 \int_0^t (1 - 5*s^2)* u_k(s) ds
## ----------------------------------------------------------------------
## Arguments: x, nu, expon.scaled: as besselI()
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 22 Nov 2008, 15:22
stopifnot(k.max == round(k.max),
0 <= k.max, k.max <= 4)
z <- x/nu # -> computing I_nu(z * nu)
sz <- sqrt(1 + z^2) ## << "FIXME": use hypot(.)
t <- 1/sz
eta <- (if(expon.scaled) ## <==> * exp(-x) == exp(- z * nu) <==> "eta - z"
## sz - z = sqrt(1 + z^2) - z = 1/(sqrt(..) + z)
1/(z + sz) else sz) +
log(z / (1 + sz))
## I_nu(nu * z) ~ 1/sqrt(2*pi*nu) * exp(nu*eta)/(1+z^2)^(1/4) *
## * {1 + u_1(t)/nu + u_2(t)/nu^2 + ... }
if(k.max == 0)
d <- 0
else { ## k.max >= 1 --- Find the Debye polynomials u_j(t) in
## /usr/local/app/R/R_local/src/QRMlib/src/bessel.c
t2 <- t^2
u1.t <- (t*(3 - 5*t2))/24
d <-
if(k.max == 1) {
u1.t/nu
} else { ## k.max >= 2
u2.t <- ##(81*t^2 - 462*t^4 + 385*t^6)/1152
t2*(81 + t2*(-462 + t2*385)) / 1152
if(k.max == 2)
(u1.t + u2.t/nu)/nu
else { ## k.max >= 3
u3.t <- t*t2*(30375 +
t2*(-369603 +
t2*(765765 - t2*425425)))/ 414720
if(k.max == 3)
(u1.t + (u2.t + u3.t/nu)/nu)/nu
else { ## k.max >= 4
u4.t <- t2*t2*(4465125 +
t2*(-94121676 +
t2*(349922430 +
t2*(-446185740 + t2*185910725))))/39813120
if(k.max == 4)
(u1.t + (u2.t + (u3.t + u4.t/nu)/nu)/nu)/nu
}
}
}
}
if(log) {
log1p(d) + nu*eta - (log(sz) + log(2*pi*nu))/2
} else {
(1+d) * exp(nu*eta)/sqrt(2*pi*nu*sz)
}
} ## {besselI.nuAsym}
besselK.nuAsym <- function(x, nu, k.max, expon.scaled=FALSE, log=FALSE)
{
## Purpose: Asymptotic expansion of Bessel K_nu(x) function
## when BOTH nu and x are large
## Abramowitz & Stegun , p.378, __ 9.7.8. __
##
## K_nu(nu * z) ~ sqrt(pi/(2*nu)) * exp(-nu*eta)/(1+z^2)^(1/4) *
## * {1 - u_1(t)/nu + u_2(t)/nu^2 - ... }
## { see besselI.nuAsym() above, for t, eta, u_k ...}
## ----------------------------------------------------------------------
## Arguments: x, nu, expon.scaled: as besselK()
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 23 Nov 2009, 13:39
stopifnot(k.max == round(k.max),
0 <= k.max, k.max <= 4)
z <- x/nu # -> computing K_nu(z * nu)
sz <- sqrt(1 + z^2) ## << "FIXME": use hypot(.)
t <- 1/sz
eta <- (if(expon.scaled) ## <==> * exp(-x) == exp(- z * nu) <==> "eta - z"
## sz - z = sqrt(1 + z^2) - z = 1/(sqrt(..) + z)
1/(z + sz) else sz) +
log(z / (1 + sz))
## K_nu(nu * z) ~ [see above]
if(k.max == 0)
d <- 0
else { ## k.max >= 1 --- Find the Debye polynomials u_j(t) in
## /usr/local/app/R/R_local/src/QRMlib/src/bessel.c
t2 <- t^2
u1.t <- (t*(3 - 5*t2))/24
d <-
if(k.max == 1) {
- u1.t/nu
} else { ## k.max >= 2
u2.t <- ##(81*t^2 - 462*t^4 + 385*t^6)/1152
t2*(81 + t2*(-462 + t2*385)) / 1152
if(k.max == 2)
(- u1.t + u2.t/nu)/nu
else { ## k.max >= 3
u3.t <- t*t2*(30375 +
t2*(-369603 +
t2*(765765 - t2*425425)))/ 414720
if(k.max == 3)
(- u1.t + (u2.t - u3.t/nu)/nu)/nu
else { ## k.max >= 4
u4.t <- t2*t2*(4465125 +
t2*(-94121676 +
t2*(349922430 +
t2*(-446185740 + t2*185910725))))/39813120
if(k.max == 4)
(- u1.t + (u2.t + (-u3.t + u4.t/nu)/nu)/nu)/nu
else stop("k.max >= 5 not yet implemented")
}
}
}
}
if(log) {
log1p(d) - nu*eta - (log(sz) - log(pi/(2*nu)))/2
} else {
(1+d) * exp(-nu*eta)*sqrt(pi/(2*nu * sz))
}
} ## {besselK.nuAsym}
Try the Bessel package in your browser
Any scripts or data that you put into this service are public.
Bessel documentation built on May 2, 2019, 4:38 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions besselIs bI besselIasym besselI.ftrms besselI.nuAsym besselK.nuAsym
Documented in besselIasym besselI.ftrms besselI.nuAsym besselIs besselK.nuAsym bI
###--- This used to be part /u/maechler/R/MM/NUMERICS/bessel-fn.R
#### "Modified" Bessel Function I_nu(x)
#### ----------------------------------------------
## besselI() - Definition as infinite sum -- working for "mpfr" numbers, too!
## ---------
besselIs <-
function(x, nu, nterm = 800, expon.scaled = FALSE, log = FALSE,
Ceps = if(isNum) 8e-16 else 2^(- x@.Data[[1]]@prec))
{
## Purpose: besselI() primitively
## ----------------------------------------------------------------------
## Arguments: (x,nu) as besselI; nterm: number of terms for "infinite" sum
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 21 Oct 2002, 21:59
if(length(nu) > 1)
stop(" 'nu' must be scalar (length 1)!")
j <- (nterm-1):0 # sum smallest first!
n <- length(x)
if(n == 0) return(x)
if(is(nu, "mpfr")) x <- mpfr(x, precBits = max(64, .getPrec(nu)))
l.s.j <- outer(j, (x/2), function(X,Y) X*2*log(Y))##-> {nterm x n} matrix
##
isNum <- is.numeric(x) || is.complex(x)
## improve accuracy for lgamma(j+1) for "mpfr" numbers
## -- this is very important [evidence: e.g. bI(10000, 1)]
if(is(l.s.j, "mpfr"))
j <- mpfr(j, precBits = max(sapply(l.s.j@.Data, slot, "prec")))
else if(!isNum) j <- as(j, class(x))
## underflow (64-bit AMD) for x > 745.1332
## for large x, this already overflows to Inf :
## s.j <- outer(j, (x^2/4), function(X,Y) Y^X) ##-> {nterm x n} matrix
## s.j <- s.j / (gamma(j+1) * gamma(nu+1 + j)) but without overflow :
log.s.j <-
if(expon.scaled)
l.s.j - rep(x, each=nterm) - lgamma(j+1) - lgamma(nu+1 + j)
else
l.s.j - lgamma(j+1) - lgamma(nu+1 + j)
s.j <-
if(log) # NB: lsum() works on whole matrix
lsum(log.s.j) # == log(sum_{j} exp(log.s.j) )
else ## log I_nu(x) -- trying to avoid overflow/underflow for large x OR large nu
## log(s.j) ; e..x <- exp(-x) # subnormal for x > 1024*log(2) ~ 710;
exp(log.s.j)
if(log) {
if(any(i0 <- x == 0)) s.j[i0] <- 0
if(any(lrgS <- log.s.j[1,] > log(Ceps) + s.j))
lapply(x[lrgS], function(x)
warning(sprintf("bI(x=%g): 'nterm' may be too small", x),
call.=FALSE))
nu*log(x/2) + s.j
} else { ## !log
s <- colSums(s.j)
if(any(i0 <- x == 0)) s[i0] <- 0
if(!all(iFin <- is.finite(s)))
stop(sprintf("infinite s for x=%g", x[!iFin][1]))
if(any(lrgS <- s.j[1,] > Ceps * s))
lapply(x[lrgS], function(x)
warning(sprintf("bI(x=%g): 'nterm' may be too small", x),
call.=FALSE))
(x/2)^nu * s
}
}
## old name [back compatibility]:
bI <- function(x, nu, nterm = 800, expon.scaled = FALSE, log = FALSE,
Ceps = if(isNum) 8e-16 else 2^(- x@.Data[[1]]@prec))
{
.Deprecated("besselIs")
isNum <- is.numeric(x) || is.complex(x)
besselIs(x, nu, nterm=nterm, expon.scaled=expon.scaled, log=log, Ceps=Ceps)
}
###--------------- besselI() for large x or also large nu ---
###
### see --> ~/R/MM/NUMERICS/bessel-large-x.R for code usage
### ================================
besselIasym <- function(x, nu, k.max = 10, expon.scaled=FALSE, log=FALSE)
{
## Purpose: Asymptotic expansion of Bessel I_nu(x) function x -> oo
## by Abramowitz & Stegun (9.7.1), p.377 :
##
## I_a(z) = exp(z) / sqrt(2*pi*z) * f(z,..) where
## f(z,..) = 1 - (mu-1)/ (8*z) + (mu-1)(mu-9)/(2! (8z)^2) - ...
## = 1- (mu-1)/(8z)*(1- (mu-9)/(2(8z))*(1-(mu-25)/(3(8z))*..))
## where mu = 4*a^2 *and* |arg(z)| < pi/2
## ----------------------------------------------------------------------
## Arguments: x, nu, expon.scaled: as besselI()
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 11 Apr, 21 Nov 2008
## Note that for each x the series eventually *DIVERGES*
## it should be "stopped" as soon as it *converged* (:-)
stopifnot(k.max == round(k.max))
## I_\nu(z) = exp(z) / sqrt(2*pi*z) * f(z, nu)
## First compute f(x, nu) to order k.max
## f <- 1
d <- 0 ## d = 1 - f <==> f = 1 - d
if(k.max >= 1) {
## m. <- 4*nu^2
x8 <- 8*x
for(k in k.max:1) {
## mu <- 4*nu^2; d <- (1 - d)*(mu - (2*k-1)^2)/(k*x8)
## mu - (2k-1)^2 = (2nu - (2k-1)) (2nu + (2k-1)) =
## = (2(nu-k)+1)(2(nu+k)-1)
d <- (1 - d)*((2*(nu-k)+1)*(2*(nu+k)-1))/(k*x8)
}
}
if(log) {
## f = 1 - d ==> log(f) = log1p(-d) :
(if(expon.scaled) log1p(-d) else x + log1p(-d)) - log(2*pi*x) / 2
} else {
(if(expon.scaled) (1-d) else exp(x)*(1-d)) / sqrt(2*pi*x)
}
}
besselI.ftrms <- function(x, nu, K = 20)
{
## Purpose: all *Terms* in besselIasym()
## by Abramowitz & Stegun (9.7.1), p.377 :
## f(x,..) = 1 - (mu-1)/ (8*x) + (mu-1)(mu-9)/(2! (8x)^2) - ...
## = 1- (mu-1)/(8x)*(1- (mu-9)/(2(8x))*(1-(mu-25)/(3(8x))*..))
## where mu = 4*nu^2
## ----------------------------------------------------------------------
## Arguments: x, nu, expon.scaled: as besselI()
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 11 Apr, 21 Nov 2008
stopifnot(length(x) == 1, length(nu) == 1, length(K) == 1,
x > 0, K == round(K), K >= 0)
## I_\nu(x) = exp(x) / sqrt(2*pi*x) * f(x, nu)
kk <- seq_len(K)
mu <- 4*nu^2
x8 <- 8*x
multf <- - (mu - (2*kk-1)^2)/(kk*x8)
ser <- cumprod(multf)
## now, for k-term approx. f_k of f(x,nu) we have
## f_k = 1 + sum(ser[1:k]), i.e. ## (f_k)_k = 1 + cumsum(c(0, ser))[k+1] for k= 0,1,...,K
ser
}
###----------------------------------------------------------------------------
### When BOTH nu and x are large :
besselI.nuAsym <- function(x, nu, k.max, expon.scaled=FALSE, log=FALSE)
{
## Purpose: Asymptotic expansion of Bessel I_nu(x) function
## when BOTH nu and x are large
## Abramowitz & Stegun , p.378, __ 9.7.7. __
##
## I_nu(nu * z) ~ 1/sqrt(2*pi*nu) * exp(nu*eta)/(1+z^2)^(1/4) *
## * {1 + u_1(t)/nu + u_2(t)/nu^2 + ... }
## where t = 1 / sqrt(1 + z^2),
## eta = sqrt(1 + z^2) + log(z / (1 + sqrt(1+z^2)))
##
##
## and u_k(t) from p.366 __ 9.3.9 __
## u0(t) = 1
## u1(t) = (3*t - 5*t^3)/24
## u2(t) = (81*t^2 - 462*t^4 + 385*t^6)/1152
## ...
## with recursion 9.3.10 for k = 0, 1, .... :
##
## u_{k+1}(t) = t^2/2 * (1 - t^2) * u'_k(t) +
## 1/8 \int_0^t (1 - 5*s^2)* u_k(s) ds
## ----------------------------------------------------------------------
## Arguments: x, nu, expon.scaled: as besselI()
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 22 Nov 2008, 15:22
stopifnot(k.max == round(k.max),
0 <= k.max, k.max <= 4)
z <- x/nu # -> computing I_nu(z * nu)
sz <- sqrt(1 + z^2) ## << "FIXME": use hypot(.)
t <- 1/sz
eta <- (if(expon.scaled) ## <==> * exp(-x) == exp(- z * nu) <==> "eta - z"
## sz - z = sqrt(1 + z^2) - z = 1/(sqrt(..) + z)
1/(z + sz) else sz) +
log(z / (1 + sz))
## I_nu(nu * z) ~ 1/sqrt(2*pi*nu) * exp(nu*eta)/(1+z^2)^(1/4) *
## * {1 + u_1(t)/nu + u_2(t)/nu^2 + ... }
if(k.max == 0)
d <- 0
else { ## k.max >= 1 --- Find the Debye polynomials u_j(t) in
## /usr/local/app/R/R_local/src/QRMlib/src/bessel.c
t2 <- t^2
u1.t <- (t*(3 - 5*t2))/24
d <-
if(k.max == 1) {
u1.t/nu
} else { ## k.max >= 2
u2.t <- ##(81*t^2 - 462*t^4 + 385*t^6)/1152
t2*(81 + t2*(-462 + t2*385)) / 1152
if(k.max == 2)
(u1.t + u2.t/nu)/nu
else { ## k.max >= 3
u3.t <- t*t2*(30375 +
t2*(-369603 +
t2*(765765 - t2*425425)))/ 414720
if(k.max == 3)
(u1.t + (u2.t + u3.t/nu)/nu)/nu
else { ## k.max >= 4
u4.t <- t2*t2*(4465125 +
t2*(-94121676 +
t2*(349922430 +
t2*(-446185740 + t2*185910725))))/39813120
if(k.max == 4)
(u1.t + (u2.t + (u3.t + u4.t/nu)/nu)/nu)/nu
}
}
}
}
if(log) {
log1p(d) + nu*eta - (log(sz) + log(2*pi*nu))/2
} else {
(1+d) * exp(nu*eta)/sqrt(2*pi*nu*sz)
}
} ## {besselI.nuAsym}
besselK.nuAsym <- function(x, nu, k.max, expon.scaled=FALSE, log=FALSE)
{
## Purpose: Asymptotic expansion of Bessel K_nu(x) function
## when BOTH nu and x are large
## Abramowitz & Stegun , p.378, __ 9.7.8. __
##
## K_nu(nu * z) ~ sqrt(pi/(2*nu)) * exp(-nu*eta)/(1+z^2)^(1/4) *
## * {1 - u_1(t)/nu + u_2(t)/nu^2 - ... }
## { see besselI.nuAsym() above, for t, eta, u_k ...}
## ----------------------------------------------------------------------
## Arguments: x, nu, expon.scaled: as besselK()
## ----------------------------------------------------------------------
## Author: Martin Maechler, Date: 23 Nov 2009, 13:39
stopifnot(k.max == round(k.max),
0 <= k.max, k.max <= 4)
z <- x/nu # -> computing K_nu(z * nu)
sz <- sqrt(1 + z^2) ## << "FIXME": use hypot(.)
t <- 1/sz
eta <- (if(expon.scaled) ## <==> * exp(-x) == exp(- z * nu) <==> "eta - z"
## sz - z = sqrt(1 + z^2) - z = 1/(sqrt(..) + z)
1/(z + sz) else sz) +
log(z / (1 + sz))
## K_nu(nu * z) ~ [see above]
if(k.max == 0)
d <- 0
else { ## k.max >= 1 --- Find the Debye polynomials u_j(t) in
## /usr/local/app/R/R_local/src/QRMlib/src/bessel.c
t2 <- t^2
u1.t <- (t*(3 - 5*t2))/24
d <-
if(k.max == 1) {
- u1.t/nu
} else { ## k.max >= 2
u2.t <- ##(81*t^2 - 462*t^4 + 385*t^6)/1152
t2*(81 + t2*(-462 + t2*385)) / 1152
if(k.max == 2)
(- u1.t + u2.t/nu)/nu
else { ## k.max >= 3
u3.t <- t*t2*(30375 +
t2*(-369603 +
t2*(765765 - t2*425425)))/ 414720
if(k.max == 3)
(- u1.t + (u2.t - u3.t/nu)/nu)/nu
else { ## k.max >= 4
u4.t <- t2*t2*(4465125 +
t2*(-94121676 +
t2*(349922430 +
t2*(-446185740 + t2*185910725))))/39813120
if(k.max == 4)
(- u1.t + (u2.t + (-u3.t + u4.t/nu)/nu)/nu)/nu
else stop("k.max >= 5 not yet implemented")
}
}
}
}
if(log) {
log1p(d) - nu*eta - (log(sz) - log(pi/(2*nu)))/2
} else {
(1+d) * exp(-nu*eta)*sqrt(pi/(2*nu * sz))
}
} ## {besselK.nuAsym}
Try the Bessel package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.