mbesseli
==============================================

Purpose
----------------

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

Format
----------------
.. function:: y = mbesseli(x, n, alpha)
              y = mbesseli0(x)
              y = mbesseli1(x)
              y = mbesselei(x, n, alpha)
              y = mbesselei0(x)
              y = mbesselei1(x)

    :param x: abscissae.
    :type x: Kx1 vector

    :param n: highest order.
    :type n: scalar

    :param alpha: :math:`0 <= alpha < 1`.
    :type alpha: scalar

    :return y: evaluations of the modified Bessel or the exponentially scaled modified Bessel
        of the first kind of the nth order.

    :rtype y: KxN matrix

Examples
----------------
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
:math:`-π` and :math:`π` in a circular manner. The model is

.. math::

    y = \mu + G(X\beta)

where :math:`\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 :math:`G()` is a
function mapping :math:`X\beta` onto the :math:`[ -\pi, \pi ]` interval.

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

.. math::

    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. :func:`QNewton` is used here.

:math:`κ` 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 :func:`atan` function is used in :math:`G()` to map :math:`X\beta` to the :math:`[ -\pi, \pi ]` interval
as suggested by Fisher, N.I. ... 1993, 158.

::

    // Define function
    proc G(u);
       retp(2*atan(u));
    endp;

    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])) -
          sumc(exp(b[1])*(cos(dev)-1))));
    endp;

    // 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;
    print b~sqrt(diag(cov));

Remarks
-------

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*:

.. math::

   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.

Source
------

ribesl.src

.. seealso:: Functions :func:`besselj`, :func:`besselk`, :func:`bessely`