Title

    integrate() -- Numerical integration of Mata functions

Syntax

    real scalar integrate(&function(), real scalar lower, real scalar upper
        |, real scalar quadpts, real rowvector xarg)


    where

       &function() is the function that is integrated, the first argument of
           this function will be the unit of integraion.

       xarg is an optional argument that is passed to the integrand.

       lower is the lower limit of the integral, a "." value represents
           -infinity

       upper is the upper limit of the integral, a "." value represents
           +infinity

       quadpts specifies the number of quadrature points used in the
           numerical integration; the default is 60.

Description

    integrate is an implementation of three numerical integration algorithms:
    Gauss-Legendre quadrature; Gauss-Hermite quadrature ; and Gauss-Laguerre.
    Gauss-Legendre quadrature is used for the definite integrals,
    Gauss-Hermite quadrature is used for the indefinite integral between
    -infinity and +infinity; and Gauss-Laguerre quadrature is used for the
    indefinite integral is between 0 and +infinity. Any limits can be chosen
    and the command will select a combination of quadrature techniques to
    calculate the result. The current command does not attempt to look at
    numerical errors but the user can alter the number of quadrature points
    to inspect any numerical instabilities.

    The number of quadrature points can be chosen to be any number above 1
    but the larger this number the slower the algorithm.  There is no upper
    limit because the quadrature points are chosen by calculating the
    eigenvalues and eigenvectors of a companion matrix.

Examples

    Below we use the integrate() function to evaluate the integration of x
    between -1 and 1:


real rowvector f1(real rowvector x)
{
  return(x)
}


    Then to integrate this function type


: integrate(&f1(), -1, 1)
2.27336e-15


    The analytical solution is 0 and given this is a numerical evaluation the
    answer is a very small number. In fact this accuracy cannot be improved
    by increasing the quadrature points as this is the accuracy of adding
    numbers.

    The integrate function can handle double integrals but will become
    increasingly more computer intensive with each additional integral.
    Double integration will square the number of operations. As an example
    integrate the function x+y with respect to x and then y over the unit
    square.


real rowvector fin(real rowvector x, real rowvector y)
{
  return(x:+y)
}
real rowvector fout(real rowvector y)
{
  for(i=1; i<=cols(y);i++) {
    if (i==1) f=integrate(&fin(), -1, 1, 40, y[i])
    else f = f, integrate(&fin(), -1, 1, 40, y[i])
  }
  return(f)
}


    Note that the integrate function requires the integrand to return a
    rowvector, however, integrate() only returns a real scalar. The function
    fout() above loops over the elements of y and produces the integration
    with respect to x for each value of y and then returns these multiple
    elements as a rowvector. Note that because fin() requires two arguments
    the integrate() function requires five arguments including the number of
    quadrature points specified as 40.  To evaluate the integral type the
    following


: integrate(&fout(), -1, 1)
  -1.14925e-16


    Again the answer is very close to 0.

    The last example is one where the upper limit of the inner integral is a
    function of y. The integral to be evaluates is the integration from 0 to
    2 with respect to y of the integral from 0 to y^2 of 6xy with respect to
    x. The example follows the same as the first double integral with the
    exception in the fout2() function that the upper limit is y[i]^2.


real rowvector fin2(real rowvector x, real rowvector y)
{
  return(6:*x:*y)
}
real rowvector fout2(real rowvector y)
{
  for(i=1; i<=cols(y);i++) {
    if (i==1) f=integrate(&fin2(), 0, y[i]^2, 40, y[i])
    else f = f, integrate(&fin2(), 0, y[i]^2, 40, y[i])
  }
  return(f)
}


    This example was chosen because it is easy to confirm analytically that
    the answer is 32. The following command shows that the integrate gets the
    correct answer


: integrate(&fout2(), 0, 2)
  32



Author

Adrian Mander, MRC Biostatistics Unit, Cambridge, UK.

Email adrian.mander@mrc-bsu.cam.ac.uk