```
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

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

Author

Adrian Mander, MRC Biostatistics Unit, Cambridge, UK.