Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Numerical Integration Mata

    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.

  • #2
    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
    Code:
    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.
    Code:
    . 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

    Comment


    • #3
      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.

      Code:
      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);

      Comment


      • #4
        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
        Code:
        expsum((0.5,h,S,beta'))

        Comment


        • #5
          Hi William,

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

          Code:
          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?

          Comment


          • #6
            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,
            Code:
            real rowvector expsum(real scalar u, real rowvector store)
            needs to become
            Code:
            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.

            Comment


            • #7
              Hi William,

              Thanks for the help. My codes work right now!

              Best,
              Han

              Comment


              • #8
                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:

                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!

                Comment


                • #9
                  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.

                  Comment


                  • #10
                    Hi William,

                    Here is the code:

                    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

                    Comment


                    • #11
                      When I run your code, several problems arise. I've stopped after the critical one.
                      Code:
                      . 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
                      Code:
                              /* 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.
                      Last edited by William Lisowski; 03 Jun 2020, 12:45.

                      Comment


                      • #12
                        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.
                        Attached Files

                        Comment


                        • #13
                          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:
                          Code:
                          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

                          Comment


                          • #14
                            Code:
                            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?
                            Last edited by William Lisowski; 03 Jun 2020, 15:42.

                            Comment


                            • #15
                              Hi Sven and William,

                              Thanks for the help. The code works right now.
                              Last edited by Han Ng; 03 Jun 2020, 16:17.

                              Comment

                              Working...
                              X