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

[Andrea Discacciati]

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.

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.

[Christophe Kolodziejczyk]

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

 
    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

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

[Andrea Discacciati]

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:

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.

[Christophe Kolodziejczyk]

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.

[Andrea Discacciati]

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!

[Joseph Coveney]

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

[Andrea Discacciati]

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:

. 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!

[Joseph Coveney]

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.

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

[Joseph Coveney]

And it might depend upon indexing vectors, too. For example, Mata will run this

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

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

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

[Andrea Discacciati]

Thank you once again, Joseph. Very helpful!