Computes modified and exponentially scaled modified Bessels of the first kind of the nth order.


y = mbesseli(x, n, alpha)
y = mbesseli0(x)
y = mbesseli1(x)
y = mbesselei(x, n, alpha)
y = mbesselei0(x)
y = mbesselei1(x)
  • x (Kx1 vector) – abscissae.

  • n (scalar) – highest order.

  • alpha (scalar) – \(0 <= alpha < 1\).


y (KxN matrix) – evaluations of the modified Bessel or the exponentially scaled modified Bessel of the first kind of the nth order.


This example produces estimates for the “circular” response regression model (Fisher, N.I. Statistical Analysis of Circular Data. NY: Cambridge University Press, 1993.), where the dependent variable varies between \(-π\) and \(π\) in a circular manner. The model is

\[y = \mu + G(X\beta)\]

where \(\beta\) is a vector of regression coefficients, X a matrix of independent variables with a column of 1’s included for a constant, and y a vector of “circular” dependent variables, and where \(G()\) is a function mapping \(X\beta\) onto the \([ -\pi, \pi ]\) interval.

The log-likelihood for this model is from Fisher, N.I. … 1993, 159:

\[log⁡L=−N \times ln⁡(I_0(\kappa)) + \kappa ⁢\sum_i^N ⁢cos⁡(y_i−\mu−G(X_i\beta))\]

To generate estimates it is necessary to maximize this function using an iterative method. QNewton() is used here.

\(κ\) is required to be nonnegative and therefore in the example below, the exponential of this parameter is estimated instead. Also, the exponentially scaled modified Bessel is used to improve numerical properties of the calculations.

The atan() function is used in \(G()\) to map \(X\beta\) to the \([ -\pi, \pi ]\) interval as suggested by Fisher, N.I. … 1993, 158.

// Define function
proc G(u);

proc lpr(b);
   local dev;
   ** b[1] - kappa
   ** b[2] - mu
   ** b[3] - constant
   ** b[4:rows(b)] - coefficients
   dev = y - b[2]- G(b[3] + x * b[4:rows(b)]);
   retp(rows(dev)*ln(mbesselei0(exp(b[1])) -

// Load data matrix
loadm data;

// Separate x0 and y0
y0 = data[., 1];
x0 = data[., 2:cols(data)];

// Create constant
b0 = 2*ones(cols(x0), 1);

// Call optimization procedure
{ b, fct, grd, ret } = QNewton(&lpr, b0);

cov = invpd(hessp(&lpr, b));

print "estimates standard errors";
print b~sqrt(diag(cov));


For the functions that permit you to specify the order, the returned matrix contains a sequence of modified or exponentially scaled modified Bessel values of different orders. For the ith row of y:

\[y[i,.] = I_{\alpha}(x[i])\; I_{\alpha+1}(x[i])\;...\;I_{\alpha+n-1}(x[i])\]

The remaining functions generate modified Bessels of only the specified order.

The exponentially scaled modified Bessels are related to the unscaled modifed Bessels in the following way:

mbesselei0(x) = exp(-x) * mbesseli0(x)

The use of the scaled versions of the modified Bessel can improve the numerical properties of some calculations by keeping the intermediate numbers small in size.



See also

Functions besselj(), besselk(), bessely()