rfunctionderivativebessel-functions

Computing logarithmic derivatives of modified Bessel function of the second kind


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.


Solution

  • 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