Nonlinear MLE with Many Variables

[Samuel Van Gorden]

I'm trying to estimate the coefficients and variance of the following model:

y = X*W + ε(0,σ²)

Where W is a vector of transformed coefficients given by:

w_i = exp(b_i) / sum(exp(b_i)).

This constrains the coefficients to be between 0 and 1 and to sum to 1 (i.e. they represent proportions of each variable).

I have tried doing this in MATA using moptimize() with the d0 model and transforming the coefficients manually, as well as using the constraint function for the equality constraint and a penalty function for the inequality constraints. Since this didn't seem to work I have moved on to STATA and using the various ml*/nl* methods. nlsur seemed promising except that it didn't seem possible to specify a vector of parameters. In other words, I would have to create a new b_i for every variable in my dataset which is unwieldy since we're usually using about 50 variables.

Can anyone provide some help as to what I'm looking for? I've been hammering away at this for weeks now with no results and it just seems like it should be simpler than this, especially considering it can be very easily achieved in MATLAB just using fminunc.

[Joro Kolev]

General advice is that if you know how to do this in Matlab, just do it in Matlab.

If you really want to do this in Mata, show your attempts, and explain how your attempts fail.

Finally if you spell out your model in scalar notation, I might be able to show you how to do this in -nlsur- , -nl-.

[Samuel Van Gorden]

Here's my MATA solution:

function linregeval_d0(transmorphic M, real scalar todo, real rowvector b, fv, g, H)
{
    x = moptimize_util_userinfo(M, 1);
    y = moptimize_util_depvar(M, 1);
    num_p = cols(b);
    w = exp(b[1..(num_p-1)]) / sum(exp(b[1..(num_p-1)]));
    sig = sqrt(b[num_p]);
    ll = lnnormalden(y, x*w', sig);
    fv = moptimize_util_sum(M, ll);
}

real colvector synthOptimize3(real matrix x1, real matrix y1)
{
    M = moptimize_init();
    moptimize_init_evaluator(M, &linregeval_d0());
    moptimize_init_evaluatortype(M, "d0");
    moptimize_init_userinfo(M, 1, x1);
    moptimize_init_depvar(M, 1, y1);
    moptimize_init_eq_indepvars(M, 1, x1);
    moptimize_init_eq_indepvars(M, 2, "");
    moptimize_init_eq_cons(M, 1, "off");
    moptimize_init_eq_coefs(M, 1, J(1,cols(x1),10));
    moptimize_init_eq_coefs(M, 2, 0.7061);
    moptimize_init_search(M, "off");
    
    moptimize(M);
    return((moptimize_result_eq_coefs(M))');
}

This seems to work some of the time and even gives results that are fairly close to what I get in Matlab. However, I always seem to get "convergence not achieved" on the passing runs after the max # of iterations - based on the likelihood values it seems like the likelihood function is asymptotic rather than actually achieving a local max. I understand this could be due to the data I'm using which is fine as I get a similar message in Matlab.

The failing runs fail for "could not calculate numerical derivatives -- discontinuous region with missing values encountered" before hitting the max # of iterations.


The scalar model would be the following:



where

and

for all observations.

However, I want to know if it's possible to specify w in vector form instead of having to write out an equation with 50 terms in the nlsum command line. It seems like this is how the normal ml function and moptimize() works (i.e. the parameters are passed to and from the evaluating function in the vector b).

[Samuel Van Gorden]

Joro Kolev does the above give you a better idea of the solution I may need? I'm trying to create a STATA package so using Matlab is not an option.

[FernandoRios]

Hi Samuel
If you are using "lf" estimators, you will need add each transformation individually. So each estimated coefficient will be transformed in your w_i's
I think, however, that "d0" family indicators may allow you to do this kind of transformation all at once, Because one needs to evaluate the combination by hand.

For example if you do:
ml model d1 (a=x1 x2 x3)

Internally, there will be a vector a with coefficients for x1 x2 x3 and the constant. You can transform them before evaluating the combination of interest (w0+w1x1+w2x2+w3x3)

You may want to check the manual to see if that fits your needs.
F

[Samuel Van Gorden]

FernandoRios thanks for your reply. Would you be able to explain a bit more how I would actually do the transformation? Is it different from what I have in line 4 of linregeval_d0? It seems like there's something weird going on because for certain datasets I obtain coefficient estimates fairly close to (but not exactly the same as) those I get in Matlab. But then for others the algorithm just errors out for "discontinuous region with missing values", whereas it completes successfully in Matlab. Is there any chance this is a bug in Mata's moptimize()?

By the way, I finally came up with a (somewhat hacky) solution in STATA:

program myeval_d0
    args todo b lnf
    tempvar xw sig
    mata: b = st_matrix("`b'")
    mata: w = exp(b[1..50]) / sum(exp(b[1..50]))
    mata: st_matrix("w", w')
    mkmat x*, matrix(X)
    mat xw = X*w
    svmat double xw, name(xw)
    mleval `sig' = `b', eq(2) scalar
    mlsum `lnf' = lnnormalden($ML_y1, xw, `sig')
    drop xw
end

mat a = J(1,50,10)
ml model d0 myeval_d0 (y=x*, noconst) /sigma, init(a 0.7061, copy) collinear
ml maximize

The results are farther off from Matlab than those obtained from Mata's moptimize(), but they still seem to give more weight to the highest weighted regressors found in Matlab than to the other regressors. However, the STATA method seems non-feasible due to the immense amount of time each iteration takes compared to moptimize().

[Samuel Van Gorden]

Hi Joro Kolev FernandoRios, any advice you can provide based on the information I gave? Sorry to bother you and thanks in advance!

[FernandoRios]

Hi Samuel
What you are doing is what I have in mind. Making the transformation in mata, and passing that to Stata.
Now i think these lines are not needed

mkmat x*, matrix(X)
    mat xw = X*w
    svmat double xw, name(xw)

You just need to add this one
mleval `xb' = `w', eq(1) scalar although i think you need to check if the variable names/colnames need to coincide or not. The rest looks fine to me. Perhaps, as extra precaution, you need to do mleval double, but I'm not sure if that is necessary.

HTH

[Samuel Van Gorden]

Hi FernandoRios. Unfortunately, I have tried doing what you mentioned and keep getting errors. It seems that mleval only understands `b', not `w'. When I try what you suggest I get the error "= invalid name" (when using `xb') and ", invalid name" (when using xb). Thank you for your responses.