{smcl} {* 29jan2014}{...} {cmd:help mata mm_integrate()} {hline} {title:Title} {p 4 4 2} {bf:mm_integrate() -- Numerical evaluation of a definite integral of a one dimensional function} {title:Syntax} {p 4 8 2} Integration using Simpson's rule (quadratic interpolation): {p 8 12 2} {it:real scalar} {cmd:mm_integrate_sr(}{it:f}{cmd:,} {it:a}{cmd:,} {it:b}{cmd:,} {it:n} [{cmd:,} {it:sc}, {it:...}]{cmd:)} {p 4 8 2} Integration using Simpson's 3/8 rule (cubic interpolation): {p 8 12 2} {it:real scalar} {cmd:mm_integrate_38(}{it:f}{cmd:,} {it:a}{cmd:,} {it:b}{cmd:,} {it:n} [{cmd:,} {it:sc}, {it:...}]{cmd:)} {p 4 8 2} where {p 12 16 2} {it:f}: {it:pointer scalar} containing address of function to be integrated; type {cmd:&}{it:funcname}{cmd:()} {p 12 16 2} {it:a}: {it:real scalar} containing lower limit of integral {p 12 16 2} {it:b}: {it:real scalar} containing upper limit of integral {p 12 16 2} {it:n}: {it:real scalar} containing number of intervals {p_end} {p 16 16 2} {it:n} must be a multiple of two for {cmd:mm_integrate_sr()} {p_end} {p 16 16 2} {it:n} must and a multiple of three for {cmd:mm_integrate_38()} {p 11 16 2} {it:sc}: {it:real scalar} indicating that function {it:f} is scalar {p 10 16 2} {it:...}: up to 10 additional arguments to pass on to function {it:f} {title:Description} {p 4 4 2} {cmd:mm_integrate_sr()} evaluates the integral of function {it:f} over a regular grid from {it:a} to {it:b} using Simpson's rule (quadratic approximation). The evaluation grid contains {it:n} equally-spaced intervals of width {it:d} = ({it:b}-{it:a})/{it:n}. Let {it:x}_{it:i} = {it:a} + {it:i}*{it:d}, then the integral is computed as {p 8 8 2} {it:d}/3 * [{it:f}({it:x}_0) + 4*{it:f}({it:x}_1) + 2*{it:f}({it:x}_2) + 4*{it:f}({it:x}_3) + 2*{it:f}({it:x}_4) + ... + 4*{it:f}({it:x}_({it:n}-1)) + {it:f}({it:x}_{it:n})] {p 4 4 2} Likewise, {cmd:mm_integrate_38()} evaluates the integral of function {it:f} using Simpson's 3/8 rule (cubic approximation), computed as {p 8 8 2} 3*{it:d}/8 * [{it:f}({it:x}_0) + 3*{it:f}({it:x}_1) + 3*{it:f}({it:x}_2) + 2*{it:f}({it:x}_3) + 3*{it:f}({it:x}_4) + 3*{it:f}({it:x}_5) + 2*{it:f}({it:x}_6) +... + 3*{it:f}({it:x}_({it:n}-2)) + 3*{it:f}({it:x}_({it:n}-1)) + {it:f}({it:x}_{it:n})] {p 4 4 2} {cmd:mm_integrate_sr()} and {cmd:mm_integrate_38()} assume {it:f} to accept a column vector of {it:x} values as first argument so that all {it:f}({it:x})'s can be computed in one call. Alternatively, if {it:f} is scalar, specify {it:sc}!=0 to call {it:f} for each {it:x} individually. Note that computing the {it:f}({it:x})'s in one call is usually much faster than using individual calls for each {it:f}({it:x}). {p 4 4 2} {it:f} cannot be a built-in function. To integrate a built-in function you need to create your own version of it. An example can be found below. Also see "Passing built-in functions" in {bf:{help m2_ftof:[M-2] ftof}}. {title:Remarks} {p 4 4 2} Example: Integrating the standard normal density from -1 to 1 {com}: function myf(x) return(normalden(x)) : I = mm_integrate_sr(&myf(), -1, 1, 1000) : I .6826894921 : reldif(I, normal(1)-normal(-1)) 5.11338e-14 : I = mm_integrate_38(&myf(), -1, 1, 999) : I .6826894921 : reldif(I, normal(1)-normal(-1)) 1.15463e-13{txt} {p 4 4 2}Of course there is no point in integrating {cmd:normalden()} as {cmd:normal()} already provides appropriate integrals. Nonetheless, let us further use {cmd:normalden()} to illustrate passing on arguments to {it:f}: {com}: function myf2(x, mean, sd) return(normalden(x, mean, sd)) : mean = 4; sd = 2.5 : mm_integrate_sr(&myf2(), 0, 3, 1000, 0, mean, sd) .2897789667 : mm_integrate_38(&myf2(), 0, 3, 999, 0, mean, sd) .2897789667{txt} {p 4 4 2}The fourth argument, {it:sc}=0, in the example above tells {cmd:mm_integrate_sr()} and {cmd:mm_integrate_38()} that it is ok to pass a vector of {it:x} values to {it:f} (the default). If the function you want to integrate does not allow element-by-element computations specify {it:sc}!=0 to call {it:f} individually for each x value. Example: {com}: real scalar myf3(real scalar x, real scalar k) > { > return(normalden(x) * ((1 - (x/k)^2) * (1 - 5*(x/k)^2))) > } : k = 1.5 : mm_integrate_sr(&myf3(), -k, k, 1000, 1, k) .1445150517 : mm_integrate_38(&myf3(), -k, k, 999, 1, k) .1445150517{txt} {p 4 4 2}Note that the above function could easily be revised to perform element-by-element computations, as is possible in many cases. Example: {com}: real matrix myf4(real matrix x, real scalar k) > { > return(normalden(x) :* ((1 :- (x/k):^2) :* (1 :- 5*(x/k):^2))) > } : k = 1.5 : mm_integrate_sr(&myf4(), -k, k, 1000, 0, k) .1445150517 : mm_integrate_38(&myf4(), -k, k, 999, 0, k) .1445150517{txt} {p 4 4 2}{cmd:mm_integrate_sr()} versus {cmd:mm_integrate_38()}: The literature claims that {cmd:mm_integrate_38()} has a smaller bias. However, as soon as {it:n} is moderately large both algorithms perform equally well. {p 4 4 2}Size of {it:n}: The larger {it:n}, the better the approximation but the slower the evaluation. In most situations, about 1000 should be sufficient. However, how accurate your results need to be depends on your application. {title:Conformability} {cmd:mm_integrate_sr(}{it:f}{cmd:,} {it:a}{cmd:,} {it:b}{cmd:,} {it:n}{cmd:,} {it:sc}{cmd:,} {it:...}{cmd:)} {it:f}: 1 {it:x} 1 {it:a}: 1 {it:x} 1 {it:b}: 1 {it:x} 1 {it:n}: 1 {it:x} 1 {it:sc}: 1 {it:x} 1 {it:...}: (depending on function {it:f}) {it:result}: 1 {it:x} 1. {cmd:mm_integrate_38(}{it:f}{cmd:,} {it:a}{cmd:,} {it:b}{cmd:,} {it:n}{cmd:,} {it:sc}{cmd:,} {it:...}{cmd:)} {it:f}: 1 {it:x} 1 {it:a}: 1 {it:x} 1 {it:b}: 1 {it:x} 1 {it:n}: 1 {it:x} 1 {it:sc}: 1 {it:x} 1 {it:...}: (depending on function {it:f}) {it:result}: 1 {it:x} 1. {title:Diagnostics} {p 4 4 2} None. {title:Source code} {p 4 4 2} {help moremata_source##mm_integrate:mm_integrate.mata} {title:Author} {p 4 4 2} Ben Jann, University of Bern, jann@soz.unibe.ch {title:Also see} {p 4 13 2} Online: {bf:{help m2_ftof:[M-2] ftof}}, {bf:{help moremata}}