I want to know how I can compute the first and second derivatives of logarithm of modified Bessel function of the second kind.
For instance, I'm interested in finding the following derivatives with respect to x:
The first derivative is equivalent to (ā/āx(K_š(x))) / K_š(x). I could use both R functions 'BesselK' (can be found in Bessel and fOptions packages) and 'BesselDK' (fOptions package) to calculate this. Is there a better alternative?
Particularly, how to calculate the second derivative above in R? I couldn't find anything specific anywhere.
Breaking out Python's sympy
library (because R doesn't have native symbolic computation facilities):
from sympy import *
nu, x = symbols('nu, x')
b = log(besselk(nu, x))
b1 = b.diff(x)
print(b1)
b2 = b1.diff(x)
print(b2)
Translating the results into R code (for convenience I'm defining a besselk
that swaps the order of the arguments to match Python)
besselk <- function(nu, x) besselK(x, nu)
blogd <- function(x, nu) {
(-besselk(nu - 1, x)/2 - besselk(nu + 1, x)/2)/besselk(nu, x)
}
blogd2 <- function(x, nu) {
(-besselk(nu - 1, x)/2 - besselk(nu + 1, x)/2)*(besselk(nu - 1, x)/2 +
besselk(nu + 1, x)/2)/besselk(nu, x)**2 + (besselk(nu, x)/2 +
besselk(nu - 2, x)/4 + besselk(nu + 2, x)/4)/besselk(nu, x)
}
Check results:
blogd(1.5, 2.5)
## [1] -2.051282
blogd2(1.5, 2.5)
## [1] 0.9375411
library(numDeriv)
lb <- function(x, nu) log(besselK(x, nu))
grad(lb, x = 1.5, nu = 2.5)
## [1] -2.051282
drop(hessian(lb, x = 1.5, nu = 2.5))
## 0.9375411
I was lazy and cut-and-pasted the output from sympy
. These functions would be made more efficient by vectorizing the calls to besselK
(e.g. bvec <- besselK(x, nu = nu + (-2:2))
) and plugging the values in to the formula (as it is, besselK
is called more times than necessary, especially in the second-derivative calculation; vectorization might not matter much, but not calling besselK
more times than necessary would make a big difference).
PS BesselK
, BesselDK
are in the now-removed-from-CRAN fAsianOptions package