Numerical Integration Mata

[Han Ng]

Hi all,

I have a question regarding numerical integration in mata. My problem is something of this sort. I am trying to integrate f(x,y,z) = x+2y+2z over the range of x from 0 to 1, fixing the value of y and z. However I am struggling to implement this in mata. Here are my codes

real rowvector fxy2(real rowvector x, real rowvector y, real rowvector z)
{
return(x:+(2:*y)+(2:*z))
}

real rowvector f_inner2(real rowvector y,real rowvector z)
{
f = integrate(&fxy2(),0,1,40,(y,z))
return (f)
}

: f_inner2(3,3)
fxy2(): 3001 expected 3 arguments but received 2
integrate_main(): - function returned error
integrate(): - function returned error
f_inner2(): - function returned error
<istmt>: - function returned error

The correct solution should be 12.5, but I am not getting that. Can someone please point out where is the problem in my code and how can I fix it? Thanks.

[William Lisowski]

For the benefit of others, this relies on on having installed the community contributed integrate package from SSC. Although it presents as a Stata command, it also includes an underlying Mata function

ssc describe integrate
ssc install integrate
ssc help mf_integrate

This seems to do what I understand you to want, but I suspect you might have presented an overly simplified version of the problem.

. mata:
------------------------------------------------- mata (type end to exit) ----------------------
: real rowvector fxy2(real rowvector x, real rowvector yz)
> {
>         return(x:+(2*yz[1])+(2*yz[2]))
> }

: real rowvector f_inner2(real scalar y,real scalar z)
> {
>     f = integrate(&fxy2(),0,1,40,(y,z))
>     return (f)
> }

: end
------------------------------------------------------------------------------------------------

. 
. mata: f_inner2(3,3)
  12.5

[Han Ng]

Hi William,

Thanks for the help. You are correct that the stuff I presented here is an overly simplified version of the problem I have. This is my actual problem, I tried the stuff you mentioned, but it doesn't seem to be working in my context.

h = 0.5
S = 2
beta = ( -.7862975501 \  .5862119147 )

real rowvector expsum(real scalar u, real rowvector store)
{
    h = store[1]
    S = store[2]
    beta = store[3..length(store)]'
    
    j = (1..S)
    huj = ((h:^j):*(u:^j)):/factorial(j)
    output = exp(huj*beta)
    return(output)
}
expsum(0.5,(h,S,beta'))


real rowvector f_inner2(real rowvector store)
{
    f = integrate(&expsum(),0,1,40,store)
    return(f)
}

f_inner2((h,S,beta'))

f_inner2((h,S,beta'))
                expsum():  3204  <tmp>[1,40] found where scalar required
        integrate_main():     -  function returned error
             integrate():     -  function returned error
              f_inner2():     -  function returned error
                 <istmt>:     -  function returned error
r(3204);

[William Lisowski]

The function being integrated is only allowed to have a single argument, which is passed in as the fifth argument to the integrate() function. But that argument can be a vector or matrix containing multiple values that the function being integrated has been written to treat as if they were distinct arguments.

Your expsum() wants two arguments. You need to incorporate the u argument into the single rowvector argument store that currently includes h, S, and beta, so that your test command shown just after the definition of expsum() would look something like this

expsum((0.5,h,S,beta'))

[Han Ng]

Hi William,

Thanks for the reply. I am confused by your statement. I tried your suggestion, but to no avail.

real rowvector expsum(real rowvector store)
{
    u = store[1]
    h = store[2]
    S = store[3]
    beta = store[4..length(store)]'
    
    j = (1..S)
    huj = ((h:^j):*(u:^j)):/factorial(j)
    output = exp(huj*beta)
    return(output)
}
 expsum((0.5,h,S,beta'))

real rowvector f_inner2(real rowvector store)
{
    f = integrate(&expsum(),0,1,40,store)
    return(f)
}

f_inner2((h,S,beta'))

Is there a syntax mistake in my code?

[William Lisowski]

In post #4 I miswrote. Please ignore it, and lets return to the code of post #3.

As a reading of help integrate reminds us, the function to be integrated takes two arguments (as it did in post #2).

The first argument will receive a rowvector containing the values of the variable you are integrating across at which the function will be evaluated to compute the integral.

The second argument is as I described it - a single argument, perhaps a rowvector containing elements that are the actual arguments to the function.

So, starting again with your code from post #3,

real rowvector expsum(real scalar u, real rowvector store)

needs to become

real rowvector expsum(real rowvector u, real rowvector store)

But then the code defining expsum(u,store) needs to be made to function when u is a rowvector of elements between 0 and 1. It currently does not when u has more than 2 elements.

[Han Ng]

Hi William,

Thanks for the help. My codes work right now!

Best,
Han

[Han Ng]

Hi William

I have another issue regarding numerical integration. I have a string argument in my function, how can I incorporate this into my integration routine? This is my code:

function m(u,str){

if (str == "square"){
return(u:^2)
}

else if (str == "exponential"){
return(exp(u))
}

else{
return(2:*u :+ 3)
}
}

real rowvector expsum(real rowvector u, real rowvector store, string scalar str)
{
    h = store[1]
    S = store[2]
    beta = store[3..length(store)]'
    
    j = (1..S)
    u2 = (u'*colnonmissing(j)):^j
    h2 = (h:^j)
    huj = u2:*h2:/factorial(j)
    output = exp(huj*beta):*m(u,str)'
    return(output')
}

real rowvector f_inner2(real rowvector store, string scalar str, real scalar zl, real scalar zr)
{
    f = integrate(&expsum(),-1*zl,zr,40,(store, str))
    return(f)
}

I know the syntax for section highlighted in green is wrong. What is the correct way to incorporate the string argument or is it even possible? Thanks!

[William Lisowski]

I'm sure I can figure out something.

Could you post a complete version of your code so I don't have to piece together the missing parts, like the mata command and corresponding end and the commands you expect to use to call f_inner2()? And if you actually ran the code (which I guess is why you know the syntax you highlighted in green is wrong), the error message you got might be informative to me.

As the Statalist FAQ suggests in item 12 (underlining added by me)

Help us to help you by producing self-contained questions with reproducible examples that explain your data, your code, and your problem. This helps yet others too, as they will find it easier to learn from your questions and the answers to them.

Perhaps you imagine I know the definitive answer offhand, but I don't. I know enough of Mata to be able to read the help mata documentation with some understanding of the context, and then I experiment, adding commands to display intermediate results as needed to figure out what Stata is seeing at the point of failure.

The answer to your question is going to be to change store from a real rowvector to a pointer rowvector, but this will involve changes to the way expsum() is defined and integrate() is called.

[Han Ng]

Hi William,

Here is the code:

clear all 
sysuse auto 


mata
mata clear
*Parameters!

data = st_data(.,"mpg")
n = length(data)
x= 0.5 /*evaluation point*/ 
h = 0.25 /*bandwidth*/
S = 4 /*deg of polynomial approximation*/
power = (1::S)'*J(n,S,1)
minx = 0
maxx = 1
str = "epanechikov"

function m2(u,str){

    if (str=="epanechikov"){
        return(0.75:*(1:-u:^2):*(abs(u):<=1))
    }
    
    else if (str=="gaussian"){
        return(1:/sqrt(2*pi()):*exp(-(u:^2):/2))    
    }
    
    else{
        return(u:^2)
    }
}

data = data :- minx
x = x :- minx
maxx = maxx :- minx
zl = min((x/h,1))
zr = min(((maxx-x)/h,1))
u = (data:-x):/h
u = select(u, u:<=zr :& u:>=-zl)

function g(u,zl,zr,S)
{
    power = (1::S)'
    output = (((colnonmissing(power))'*(u:+zl)'):^power')':*((colnonmissing(power))'*(zr:-u)')'
    /*output = power:*((u:+zl):^(power:)):*(zr:-u)*/
    return(output)
}

function dg(u,zl,zr,S)
{
    power = (1::S)'
    /*output = power:*((u:+zl):^(power:-1)):*(zr:-u) - (u:+zl):^power*/
    output = (((colnonmissing(power:-1))'*(u:+zl)'):^(power:-1)')':*((colnonmissing(power))'*(zr:-u)')':*power - (((colnonmissing(power))'*(u:+zl)'):^power')'
    return(output)
}

gu = g(u,zl,zr,S)
dgu = dg(u,zl,zr,S)

store = J(S,S,0)

for(s=1;s<=S;s++){
    R = ((u:*h):^(s:-1)):/factorial(s-1)
    R = mean(gu:*R):/h
    store[s,.] = -R
}

/*Backout beta*/
dgu = mean(dgu):/(h:^2)
beta = qrinv(store')*dgu'
beta = lusolve(store',dgu')

/*Calculate fhat*/
fhat_m = colsum(m2(u)):/(n:*h)

real rowvector expsum(real rowvector u, real rowvector store, string scalar str)
{
    h = store[1]
    S = store[2]
    beta = store[3..length(store)]'
    
    j = (1..S)
    u2 = (u'*colnonmissing(j)):^j
    h2 = (h:^j)
    huj = u2:*h2:/factorial(j)
    output = exp(huj*beta):*m2(u,str)'
    return(output')
}
str
expsum((0.5,0.6),(h,S,beta'),str)

real rowvector f_inner2(real rowvector store,scalar string str,real scalar zl, real scalar zr)
{
    f = integrate(&expsum(),-1*zl,zr,40,(store,str))
    return(f)
}

f_inner2((h,S,beta'),zl,zr)
logfhat = log(fhat_m/f_inner2((h,S,beta'),zl,zr))
logfhat
beta

[William Lisowski]

When I run your code, several problems arise. I've stopped after the critical one.

. clear all

. sysuse auto
(1978 Automobile Data)

.
.
. mata
------------------------------------------------- mata (type end to exit) ----------------------
: mata clear

: *Parameters!
invalid expression
r(3000);

:
: data = st_data(.,"mpg")

: n = length(data)

: x= 0.5 /*evaluation point*/

: h = 0.25 /*bandwidth*/

: S = 4 /*deg of polynomial approximation*/

: power = (1::S)'*J(n,S,1)
                     *:  3200  conformability error
                 <istmt>:     -  function returned error
r(3200);

: minx = 0

: maxx = 1

: str = "epanechikov"

:
: function m2(u,str){
>
>     if (str=="epanechikov"){
>         return(0.75:*(1:-u:^2):*(abs(u):<=1))
>     }
>    
>     else if (str=="gaussian"){
>         return(1:/sqrt(2*pi()):*exp(-(u:^2):/2))    
>     }
>    
>     else{
>         return(u:^2)
>     }
> }

:
: data = data :- minx

: x = x :- minx

: maxx = maxx :- minx

: zl = min((x/h,1))

: zr = min(((maxx-x)/h,1))

: u = (data:-x):/h

: u = select(u, u:<=zr :& u:>=-zl)

:
: function g(u,zl,zr,S)
> {
>     power = (1::S)'
>     output = (((colnonmissing(power))'*(u:+zl)'):^power')':*((colnonmissing(power))'*(zr:-u)')
> '
>     /*output = power:*((u:+zl):^(power:)):*(zr:-u)*/
>     return(output)
> }

:
: function dg(u,zl,zr,S)
> {
>     power = (1::S)'
>     /*output = power:*((u:+zl):^(power:-1)):*(zr:-u) - (u:+zl):^power*/
>     output = (((colnonmissing(power:-1))'*(u:+zl)'):^(power:-1)')':*((colnonmissing(power))'*(
> zr:-u)')':*power - (((colnonmissing(power))'*(u:+zl)'):^power')'
>     return(output)
> }

:
: gu = g(u,zl,zr,S)

: dgu = dg(u,zl,zr,S)

:
: store = J(S,S,0)

:
: for(s=1;s<=S;s++){
>     R = ((u:*h):^(s:-1)):/factorial(s-1)
>     R = mean(gu:*R):/h
>     store[s,.] = -R
> }

:
: /*Backout beta*/
: dgu = mean(dgu):/(h:^2)

: beta = qrinv(store')*dgu'

: beta = lusolve(store',dgu')

:
: beta
       1
    +-----+
  1 |  .  |
  2 |  .  |
  3 |  .  |
  4 |  .  |
    +-----+

First, the comments section of the Mata help files tells us that coments in Mata are limited to

        /* enclosed comment */

        // rest-of-line comment

Second, the calculation of power fails, but that has no effect on any subsequent calculations, which use local definitions of power within your functions g() and dg().

Third, beta is apparently not calculated as you expect it. I can't fix or ignore this.

Please fix this all three problems and post complete code with the corrections, including an end command following the mata code.

[Han Ng]

Hi William

Sorry for the mistakes. I am uploading my do files and data file for your reference. It should make things easier for you. Thanks and so sorry for the trouble. I have indicated in the do file the area where I need help. It should be in line 96 onwards.

[Sven-Kristjan Bormann]

Do you know which results you expect to see? I managed to get your code run with only two minor changes, but I am not sure if the results are correct.
My attempt at correcting the code is below:

real rowvector f_inner2(real rowvector store,string scalar str,real scalar zl, real scalar zr)
{
f = integrate(&expsum(),-1*zl,zr,60,store,str) //Removed the parenthesis to make it two arguments
return(f)
}

f_inner2((h,S,beta'),str,zl,zr) // Added the missing str-variable
   1.004796048 // Result of the integration
logfhat = log(fhat_m/f_inner2((h,S,beta'),str,zl,zr))
logfhat
  -.9589213956

[William Lisowski]

f = integrate(&expsum(),-1*zl,zr,60,store,str) //Removed the parenthesis to make it two arguments

The integrate() function in use here is part of the user-contributed integrate package at SSC, as mentioned in post #2, and help mf_integrate tells us it takes only 5 arguments, the fifth of which is passed along to the function to be integrated as its second argument. Because the integrate() function is provided in compiled form, I can't view the source to see what is becoming of the sixth argument passed to it in the code above.

Let me note, as I should for anyone who comes across this at a later date, that in Stata 16, Mata supports integration through the Quadrature() class. If I were to write a program that required integration, I'd spend the time to learn it and use it.

Added in edit: when I make the changes noted in post #13 to the code provided in post #12, the code does not run successfully, the integrate() function throws errors for receiving 6 arguments when just 5 were expected. Is there perhaps some other source for the integrate() function?

[Han Ng]

Hi Sven and William,

Thanks for the help. The code works right now.