Announcement

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

  • Numerical derivatives: how to pass a Stata local macro to the evaluator?

    I need to calculate numerical derivatives in mata. The function of the parameters is however stored in a local macro in Stata. I don't know how to pass this local to the evaluator.

    The following code outlines what I would like to do. Not surprisingly, this code won't work because the line code v = a is illegal. However, I hope it's clear what I'm trying to do. I also tried to use user-defined arguments for the evaluator with no luck.

    Code:
    local fun "b[1]^2/b[2]"
    
    clear mata
    mata:
    a = st_local("fun")
    a
    
    void myeval(b, v)
    {
        v = a // This should evaluate to v = b[1]^2/b[2] somehow
    }
    
    D = deriv_init()
    deriv_init_evaluator(D, &myeval())
    deriv_init_params(D, (0.5, 0.6))
     
    G = deriv(D, 1)
    G
    end
    Thank you in advance.
    Last edited by Andrea Discacciati; 20 Mar 2017, 07:29.

  • #2
    Since you're running Mata in interactive mode the easiest way should be the following

    Code:
     
        local fun "b[1]^2/b[2]"
    
    clear mata
    mata:
    
    
    void myeval(b, v)
    {
        v = `fun' // This should evaluate to v = b[1]^2/b[2] somehow
    }
    
    D = deriv_init()
    deriv_init_evaluator(D, &myeval())
    deriv_init_params(D, (0.5, 0.6))
     
    G = deriv(D, 1)
    G
    end

    Now if you want to be a bit more sophisticated and work with pointers to functions

    Code:
    local fun foo()
    
    clear mata
    
    mata:
    
    // create a function you pass to myeval which takes b as input and returns b[1]^2 /b[2]
    real scalar foo(real vector b)
    {
            return(b[1]^2/b[2])
    }
    
    
    a = st_local("fun")
    // find the adress of the function foo() which name is referred to in the local fun
    v = findexternal(a) // alternatively you could write v = findexternal("foo()")
    
    // (*v)((1,2))
    
    // Now you pass a pointer to a function to myeval, this pointer will be v and will point to foo()
    void myeval(b, pointer(real scalar function) scalar f, v) // See help [m-2] ftof for details
    {
            v = (*f)(b)
    }
    
    D = deriv_init()
    deriv_init_evaluator(D, &myeval())
    deriv_init_argument(D, 1, v)  // -> tell deriv_init that pointer v is the first additional argument of myeval()
    deriv_init_params(D, (0.5, 0.6))
     
    G = deriv(D, 1)
    G
    end

    Comment


    • #3
      Christophe: Thank you so much for your help!

      Unfortunately, I am not running Mata in interactive mode. I am sorry if my example misled you.
      I am actually writing a Stata command and inside this program I want to calculate a number of numerical derivatives.

      Something along these lines:

      Code:
      webuse auto, clear
      
      capture program drop silly
      program silly
      
      regress weight length width
      
      //funcitons I want to calculate the derivatives of
      local fun1 "b[1]^2/b[2]" // where b[1] refers to _b[length] and b[2] to _b[width]
      local fun2 "b[1]*b[2]"
      
      mata: myDeriv("`fun1'")
      mata: myDeriv("`fun2'")
      end
      
      clear mata
      mata
      void function myDeriv(string scalar fun)
      {
          b0 = st_matrix("e(b)")
          
          D = deriv_init()
          deriv_init_evaluator(D, &eval())
          deriv_init_params(D, b0)
          G = deriv(D, 1)
          G
      }
      
      void eval(b, v)
      {
          v = ???
      }
      end
      I guess my question remains: how do I pass the function I want to take the derivatives of to the evaluator now.

      Thank you once again.
      Last edited by Andrea Discacciati; 20 Mar 2017, 09:02.

      Comment


      • #4
        The answer is that you can't solve your problem the way you want, i.e. by passing the function as a Stata macro, in the non-interactive environnement. Mata does not understand macros an understands only variables and functions. The second solution I provided should solve this problem though. You have to wrap your different expressions as mata functions which are then passed to the derivatives evaluator in the same way as shown in my previous post.

        Comment


        • #5
          Thank you once again, Christophe. That's unfortunate, because it's difficult to wrap the expressions in pre-defined mata functions, as these expressions depend on the number of covariates included in the linear regression model (and this is decided by the user in the Stata command).

          Thanks!

          Comment


          • #6
            Originally posted by Andrea Discacciati View Post
            . . . these expressions depend on the number of covariates included in the linear regression model (and this is decided by the user in the Stata command).
            In what way? In other words: in the same manner that you would (presumably programmatically) create the expressions in Stata depending upon the number of covariates, you would programmatically build the same expressions in the Mata function.

            Comment


            • #7
              Thank you for your reply, Joseph.

              Yes, you are absolutely right: I should build the functions directly inside mata. The thing is, I am adapting previous code where the functions are built in Stata as, for example ([eq1]_b[x]^2/[eq2]_b[z]^3) (the actual equations are actually *much* longer, see below). So, I was thinking of substituting the different pieces with a series of subinstr to get something like (b[1]^2/b[5]^3) and pass everything to mata. But then again: I agree with you, the "right" thing to do would be to build these expressions directly inside mata.

              Side question: I was really surprised when mata complained when I ran this:

              Code:
              . mata
              ------------------------------------------------- mata (type end to exit) ----------------------------------------------------------------
              : real scalar foo(real vector b) {
              > return((exp(b[1]*(1))*(1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*1))*(1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4
              > 979))/((1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))*(1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979)))-(1+exp(b
              > [12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))*(1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979))/((1+exp(b[12]+b[8]*0+b[9
              > ]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))*(1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979)))-exp(b[1]*(1))*(1+exp(b[12]+b[8]*0+b[9]
              > *.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*1))/(1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))+1)/((exp(b[1]*1)*(1+exp(b[12
              > ]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979))*(1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*1))/(exp(b[1]*0)*(1+exp(b[12]+b
              > [8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979))*(1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))))-1))
              too many tokens
              r(3000);
              
              : }
              expression invalid
              r(3000);
              
              : 
              : end
              ------------------------------------------------------------------------------------------------------------------------------------------
              This expression (one of the many I'm dealing with) does not seem particularly long to me. -testnl-'s subroutine Deriv handles it without complaining (but it's quite slow). What's the maximum number of tokens that mata can handle? If I need to split these expressions up, it would be nice to know what the limit is.

              Thank you once again!

              Comment


              • #8
                Originally posted by Andrea Discacciati View Post
                What's the maximum number of tokens that mata can handle? If I need to split these expressions up, it would be nice to know what the limit is.
                This has come up on the list a few times in the past, and to my knowledge has never been answered. My hunch is that there is no one value, but rather the limit depends upon the complexity (i.e., nestedness) of the expression.

                For example, in the code below, I count that your expression has around 174 tokens, but running the code shows that Mata will handle a straightforward serial addition of 174 tokens without any trouble.

                So, my suggestion is to unravel the nesting rather than to try to limit the total number of tokens. You've got a lot of subexpressions of the nature of 1 + exp(b[12] + b[8] * 1 + b[9] * .5084 + b[10] * .5 + b[11] * .4979), and you might want to start with a companion function that flexibly evaluates those types of subexpression and then substitute the returned scalars into the main expression.
                Code:
                version 14.2
                
                clear *
                set more off
                
                #delimit ;
                local expression (exp(b[1]*(1))*(1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*
                .4979+b[2]+b[3]*1))*(1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979))/
                ((1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))*(1+exp(b[12]+
                b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979)))-(1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*
                .5+b[11]*.4979+b[2]+b[3]*0))*(1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*
                .4979))/((1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))*
                (1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979)))-exp(b[1]*(1))*
                (1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*1))/
                (1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))+1)/
                ((exp(b[1]*1)*(1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979))*
                (1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*1))/
                (exp(b[1]*0)*(1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979))*
                (1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))))-1));
                #delimit cr
                
                local tokens : subinstr local expression " " "", all
                foreach operator in "*" "/" "+" "-" {
                    local tokens : subinstr local tokens "`operator'" " ", all
                }
                local token_tally : word count `tokens'
                display in smcl as text `token_tally'
                
                local expression 0
                forvalues i = 1/`token_tally' {
                    local expression `expression' + 1
                }
                
                mata:
                `expression'
                end
                
                exit

                Comment


                • #9
                  And it might depend upon indexing vectors, too. For example, Mata will run this
                  Code:
                  version 14.2
                  
                  clear *
                  set more off
                  
                  local token_tally 174
                  
                  local expression 0
                  forvalues i = 1/`token_tally' {
                      local expression `expression' + b
                  }
                  
                  mata:
                  b = 1
                  `expression'
                  end
                  
                  exit
                  without a hitch, but not this
                  Code:
                  version 14.2
                  
                  clear *
                  set more off
                  
                  local token_tally 174
                  
                  local expression B[1]
                  forvalues i = 2/`token_tally' {
                      local expression `expression' + B[`i']
                  }
                  
                  mata:
                  B = J(1, `token_tally', 1)
                  `expression'
                  end
                  
                  exit
                  The numerous vector elements in your expression might be a contributing factor to the problem. You could try as a first step to assign the vector elements to simple scalars.

                  Toward that end, this runs
                  Code:
                  version 14.2
                  
                  clear *
                  set more off
                  
                  #delimit ;
                  local expression (exp(b[1]*(1))*(1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*
                  .4979+b[2]+b[3]*1))*(1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979))/
                  ((1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))*(1+exp(b[12]+
                  b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979)))-(1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*
                  .5+b[11]*.4979+b[2]+b[3]*0))*(1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*
                  .4979))/((1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))*
                  (1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979)))-exp(b[1]*(1))*
                  (1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*1))/
                  (1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))+1)/
                  ((exp(b[1]*1)*(1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979))*
                  (1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*1))/
                  (exp(b[1]*0)*(1+exp(b[12]+b[8]*1+b[9]*.5084+b[10]*.5+b[11]*.4979))*
                  (1+exp(b[12]+b[8]*0+b[9]*.5084+b[10]*.5+b[11]*.4979+b[2]+b[3]*0))))-1);
                  #delimit cr
                  
                  local expression : subinstr local expression "[" "", all
                  local expression : subinstr local expression "]" "", all
                  
                  foreach operator in "*" "/" "+" "-" {
                      local expression : subinstr local expression "`operator'" " `operator' ", all
                  }
                  
                  forvalues i = 1/12 {
                      local setup `setup' b`i' = 1;
                  }
                  
                  mata:
                  `setup'
                  `expression'
                  end
                  
                  exit
                  Last edited by Joseph Coveney; 21 Mar 2017, 21:14.

                  Comment


                  • #10
                    Thank you once again, Joseph. Very helpful!

                    Comment

                    Working...
                    X