Bessel functions appear in surprising places. Although originally studied in physics as the solution to differential equations that describe waves and heat diffusion, Bessel functions also have applications in probability, statistics, and signal processing. SAS software provides the JBESSEL function, which is a Bessel function of the first kind, and the IBESSEL function, which is a modified Bessel function of the first kind. I have previously written about these functions.
A SAS customer recently asked me about how to implement the modified Bessel function of the second kind in SAS. This Bessel function, which is often denoted by Kν(x), appears in the formulas for the generalized inverse Gaussian (GIG) distribution. The modified Bessel function of the second kind is also used in survival analysis, spatial statistics, financial modeling, and to model variance parameters in Bayesian models. Although SAS does not have a built-in "KBESSEL" function, you can define the Kν(x) function in SAS IML.
The definition of the modified Bessel function of the second kind
For the purpose of this article, the modified Bessel function of the second kind,
Kν(x), where the parameter ν is a real number and the domain of the function is x > 0.
A good way to compute the function is to use the integral definition from Equation 10.32.10 of the
NIST Digital Library of Mathematical Functions:
\(
K_\nu(x) = \frac{1}{2} \left( \frac{x}{2} \right)^\nu \int_0^\infty t^{-\nu-1} \exp\left( -t - \frac{x^2}{4t} \right) \, dt
\)
In this formula, x and ν are parameters to the integrand. The dummy variable, t, is integrated from 0 to infinity.
The integrand has a peak that can be computed analytically by using calculus. The peak occurs at
\(
t_{peak} = \frac{\sqrt{(\nu + 1)^2 + x^2} - (\nu + 1)}{2}
\)
In the SAS IML Language, the QUAD subroutine is used to compute integrals. Knowing the location of the peak is helpful because the QUAD subroutine in IML contains a PEAK= option that you can specify to assist the numerical integration routine. (Alternatively, you can perform a change of variables u = t / tpeak so that the peak is always at u=1. I leave that as an exercise.)
The optional exponential scaling factor
The documentation for the the IBESSEL function, which is the modified Bessel function of the first kind, contains an optional binary argument. By default, the call IBESSEL(nu, x) computes the function Iν(x). However, the function decays rapidly to zero as x increases, which makes the numerical computation difficult when x is moderately large (like x ≥ 50). Accordingly, a numerical trick is to compute the expression exp(x)*Iν(x), which you can do by using the syntax IBESSEL(nu, x, 1).
A similar situation occurs for the Bessel function of the second kind, Kν(x), which also decays quickly to zero for moderately large x. (Equation 10.40.2 in the NIST Digital Library shows that the function decays asymptotically like exp(-x)/sqrt(x).) When x is moderately large, it becomes numerically challenging to evaluate the function by using a standard numerical quadrature routine like CALL QUAD. The integration routine can prematurely terminate because the value of the integral is less than the algorithm's absolute error tolerance.
Accordingly, it is common in numerical software for the syntax to support an optional binary flag that causes the function to return a "scaled value," which is exp(x)*Kν(x). The exponential scaling factor is evaluated as part of the integrand. This keeps the peak of the integrand near 1 and permits the quadrature routine to evaluate the integral accurately.
Not only does this trick enable you to compute the scaled value exp(x)*Kν(x) to high precision, but it also enables you to find the value of log(Kν(x)) accurately. Taking the natural log of the scaled value gives log(exp(x)*Kν(x)) = x + log(Kν(x)), which you can solve for log(Kν(x)). This computation is shown in the next section. The log of the Bessel function is used in applications such as loglikelihood computations.
Implement modified Bessel function of the second kind in SAS
Because SAS already supports the JBESSEL and IBESSEL functions, I will reserve the name KBESSEL for a future built-in SAS implementation. Instead, I will define a SAS IML function named BesselK, which takes three parameters: the parameter ν, the variable x, and an optional binary flag. If the binary flag is 0 (the default), the function will return Kν(x), otherwise the function returns the scaled value exp(x)*Kν(x).
proc iml; /* Define the integrand function for the unscaled modified Bessel function of the second kind, K_{nu}(x) */ start BesselK_integ(t) global(g_nu, g_x); return( t**(-g_nu - 1) * exp( -t - (g_x**2)/(4*t)) ); finish; /* Define the integrand function that includes the EXP term: exp(x)*K_{nu}(x). Use this scaled function when x is moderately large or you want log(K_{nu}(x)). */ start BesselK_Scale_integ(t) global(g_nu, g_x); return( t**(-g_nu - 1) * exp(g_x -t - (g_x**2)/(4*t)) ); finish; /* Define the wrapper function */ start BesselK(nu, x, scale=0) global(g_nu, g_x); g_nu = nu; if scale = 0 then funcName = "BesselK_integ"; else funcName = "BesselK_Scale_integ"; n = nrow(x); K = j(n, 1, .); /* the result vector K_nu(x). Use missing for domain error such as x <= 0 */ tau = nu+1; do i = 1 to n; if x[i] > 0 then do; g_x = x[i]; /* Calculate the exact peak of the integrand */ pk = (sqrt(tau**2 + g_x**2) - tau) / 2; call QUAD(integral_at_x, funcName, {0 .I}) PEAK=pk; K[i] = integral_at_x; end; end; /* Multiply the integral by the leading constants */ K = 0.5 * (x / 2)##nu # K; return(K); finish; store module=(BesselK_integ BesselK_Scale_integ BesselK); QUIT; |
Let's test the function. First, let's visualize the function Kν(x) for ν=1 and x in (0, 3].
proc iml; load module=_all_; x = T(do(0,3,0.05)); nu = 1; title "Modified Bessel Function of the Second Kind (nu=1)"; K = BesselK(nu, x); call series(x, K) grid={x y}; |
The graph shows that the Kν(x) decays to 0 rapidly. Let's turn on the binary flag and visualize the function exp(x)*Kν(x) on the same domain:
title "Scaled Modified Bessel Function exp(x)*K(x) (nu=1)"; expx_K = BesselK(nu, x, 1); call series(x, expx_K) grid={x y}; |
The graph shows that the exponentially scaled function still decays, but not as rapidly. For example, whereas K1(30) ~ 2.5E-14, the scaled function at x=30 is much larger: exp(30)*K1(30) ~ 0.23.
By using the scaled version, you can compute log(Kν(x)), which is useful for various applications:
/* What good is the exp scaling? 1. Computations are more accurate for large values of x. 2. You can easily compute log(K_{nu}(x)) */ x = 50; expx_K = BesselK(nu, x, 1); logK = log(expx_K) - x; print x expx_K logK; |
If you plan to compute log(Kν(x)) regularly, you could encapsulate this logic into a self-contained module.
Summary
This article shows how to implement the modified Bessel function of the second kind, Kν(x), in SAS. The method involves integrating a complicated expression over an infinite domain. In addition, the article shows how you can compute the exponentially scaled expression exp(x)*Kν(x), which enables you to evaluate log(Kν(x)).
In a subsequent article, I will show how you can use the Kν Bessel function to evaluate the PDF and log-PDF of the generalized inverse Gaussian distribution in SAS.
