optimize issue with evaluator in a class

[Alfonso Sánchez-Peñalver]

Hi,

I have the following two functions defined as public in the class nelogit_fp.

// ML d2 Evaluator for a fixed parameters Logit
void nelogit_fp::MLEval(real scalar todo, real rowvector p, v, g, H)
{
    // y, X, and j must have been set. I'll write the check later
    
    // We need to get the parameters in the right place
    real colvector b, o, pij, pi
    real matrix XS
    
    p = p'
    if (rows(this.S) != 0)
    {
        // We have scale. Need to separate the coefficients, and scale the
        // independent variables.
        b = p[1..cols(this.X),.]
        o = p[cols(this.X)+1..rows(p),.]
        XS = this.X :/ exp(this.S*o)
    }
    else
    {
        // Don't have scale. The coefficients are all of the value function, and
        // we use the original X.
        b = p
        XS = this.X
    }
    
    // Calculations of probabilities
    pij = super.logitprob(XS, b, this.j)
    pi = rowsum(colshape(this.y :* pij, this.j))
    
    // Log-Likelihood
    this.ll = v = colsum(ln(pi))
    
    if(todo > 0)
    {
        real matrix px, dx
        real colvector wmk, dxk
        real scalar k
        // Calculation of the difference from the weighted mean of the explanatory
        // variables
        px = pij :* XS
        for(k=1; k<=cols(XS); k++)
        {
            wmk = rowsum(colshape(px[.,k], this.j))
            dxk = colshape(colshape(XS[.,k], this.j) :- wmk,1)
            if (k == 1) dx = dxk
            else dx = dx , dxk
        }
        
        g = cross(this.y, dx)
        
        if (rows(this.S) != 0)
        {
            // g = g, this.y'((-dx*b):*S)
            g = (g, cross(this.y,(-dx*b):*S))
        }
            
        
        if(todo > 1)
        {
            // Hessian
            H = cross(px,-dx)
            if (rows(this.S) != 0)
            {
                real matrix Hbo, Hoo
                // beta omega
                Hbo = cross((pij:* (1 :+ (dx*b)) :- this.y) :* XS,S)
                // omega omega
                Hoo = cross(((pij:* ((-dx*b) :- 1) :+ this.y) :* (XS*b)) :* S,S)
                // Form it
                H = ((H , Hbo) \ (Hbo' , Hoo))
            }
            _makesymmetric(H)
        }
    }
}
// Estimate Method to estimate the model
void nelogit_fp::Estimate()
{
    transmorphic SS
    SS = optimize_init()
    optimize_init_evaluator(SS, &(this.MLEval()))
    optimize_init_evaluatortype(SS, "d2")
    
    real rowvector b0
    b0 = (invsym(cross(this.X,this.X))*cross(this.X,this.y))'
    if(rows(this.S) != 0)
        b0 = b0, J(1, cols(this.S), 0)
    optimize_init_params(SS, b0)
    this.b = optimize(SS)
    this.b
}

I run the following

mata:
nl = nelogit_fp()
nl.setup(st_data(.,"choice"), st_data(.,"`xvars'"), 3, st_data(.,"`svars'"))
nl.Estimate()

end

and get the following error

. mata:
------------------------------------------------- mata (type end to exit) -------------------------------------------------------------
: nl = nelogit_fp()

: nl.setup(st_data(.,"choice"), st_data(.,"`xvars'"), 3, st_data(.,"`svars'"))

: nl.Estimate()
    nelogit_fp::MLEval():  3001  expected 6 arguments but received 1
  nelogit_fp::Estimate():     -  function returned error
                 <istmt>:     -  function returned error
(1 line skipped)
---------------------------------------------------------------------------------------------------------------------------------------
r(3001);

end of do-file

r(3001);

It is the command optimize_init_evaluator(SS, &(this.MLEval())) that throws the error. What is the problem? Is it that we cannot do this within a class?

[Hua Peng (StataCorp)]

The function is defined

 
 void nelogit_fp::MLEval(real scalar todo, real rowvector p, v, g, H)

and the function is called

 
 optimize_init_evaluator(SS, &(this.MLEval()))

and Mata gives an informative error message:

 
 nelogit_fp::MLEval():  3001  expected 6 arguments but received 1

[Alfonso Sánchez-Peñalver]

and the function is called

optimize_init_evaluator(SS, &(this.MLEval()))

Not quite. The function is pointed to in that command and called within. It's optimize_init_evaluator() who does the calling. Both the definition of MLEval() and the optimization setup follow the steps in the examples in the help files and in David M. Drukker (StataCorp)'s blog on how to use optimize() to estimate a Poisson model.

The informative message is not quite so informative. Since I'm not calling the function directly I don't know what gets passed or not to it. Now, reading Bill Gould (StataCorp)'s book (Gould (2018)), in section 12.7 he says that "pointers to class member functions are not allowed." I don't want to get into whether that should or should not be, but I suspect this is the issue, which is why I asked if we cannot do all this in a class in #1. Is this the problem? Is there a way around it? Because if there isn't and we cannot build a class that is an ML estimator and has a method that does the estimation using its members previously set... that is a very strong shortcoming of object programming within Mata!!!

Reference:

Gould, Willaim. 2018. The Mata Book: The Book for Serious Programmers and Those Who Want to Be. College Station, TX: Stata Press

[Alfonso Sánchez-Peñalver]

Sorry, I was reading #3 and what I wanted to say is that MLEval() is supposed to be pointed to when using optimize_init_evaluator(). The error that is thrown is the same as the error that is thrown when I call MLEval() without any arguments, which I guess is why Hua Peng (StataCorp) said that the function was being called there. However, it is not supposed to. It's supposed to be a pointer so that it is called from within.

In any case, this needs fixing. I believe the idea of having a class that is a ML estimator that does everything by allowing simple commands is a right one, and something that object-oriented programming should allow. We shouldn't have to go between Stata and Mata to do this, and a class provides a very simple interface for the user.

[Jeff Pitblado (StataCorp)]

optimize_init_evaluator() requires globally callable Mata functions. Class methods, even public ones, will not work.

Compose a call-thru Mata function that takes the evaluator arguments and an object from your class and calls the object's evaluator method. Register your object as a custom argument with optimize_init_argument() and this call-thru Mata function as the evaluator.

Based on Alfonso's original post, the Mata function would look something like

void nelogit_d2(real scalar todo, real rowvector p, class nelogit_fp scalar nl, v, g, H)
{
        nl.MLeval(todo, p, v, g, H)
}

The calls in nelogit_fp::Estimate() to register the custom evaluator argument and the above evaluator would look something like

        optimize_init_argument(SS, 1, nl)
        optimize_init_evaluator(SS, &nelogit_d2())

[Alfonso Sánchez-Peñalver]

Hi Jeff Pitblado (StataCorp) that is exactly the workaround that I have been working on. It is a very inelegant workaround, and the reason is that it exposes all the members that I didn't want exposing.

The variables are all private (protected really since they are declared in the uppermost class), so only class methods can manipulate them. If kept private, they can't be accessed and set from the outside evaluator.

I also have to expose an internal function that I use to calculate a probability.

Seriously, this beats the purpose of object-oriented programming, which is to provide a clean interface easy to use, do all the work within the object so that another programmer cannot mess up the process, and present with results to the user in an easy way to the user to handle.

[Alfonso Sánchez-Peñalver]

Actually, yours doesn't expose anything. Sorry about that. You gave me an easy way to go about it. Still a bit redundant, and it does expose the evaluator, which I wanted to keep private after testing. But it's a good workaround. Thanks.

Minor correction

optimize_init_argument(SS, 1, this)
optimize_init_evaluator(SS, &nelogit_d2())

You want to pass the actual instance of the object, not a general one, so that all the elements that have been set are there.

Thanks again.