st: too many variables from simple Mata iteration?

[László Sándor <sandorl@gmail.com>]

Hi all,

I have coded up an EM algorithm using -moptimize-, and it calls a
maximizing (-moptimize-) Mata function alternating with a Bayesian
updating of a variable (done in Mata) that captures the probability
that an observation is one of two classes. I myself surely don't
create new variables in this process. However, if I run this code for
a long-long time (convergence can be slow for EM), my code gets
aborted with a message that I reached the limit of the number of
variables.

Could you help me at what point are new variables created, and how
could I dropped them safely when they are never used again?

I paste my code below, though it is rather complicated to parse, so
don't feel obligated to do so. (Hopefully you can answer my question
even without understanding my code. And you're doing me a great favor
in any case.)

Thank you,

Laszlo

program emlatentclass, eclass
                version 11
                syntax varlist [if] [in], BY(varname) [CLASS(varname)
CLUSTER(varname)]
                marksample touse
                markout `touse' `varlist' `by' `class' `cluster'
                tokenize `varlist'
                local lhs `1'
                macro shift
                local rhs `*'
                local k: word count `rhs'
                qui sum `touse'
                local nobs = r(sum)
                cap drop prob_class1_`lhs'
                if "`class'" == "" {
                    gen prob_class1_`lhs' = 0.5
                }
                else {
                    gen prob_class1_`lhs' = `class'
                    replace prob_class1_`lhs' = 0.5 if prob_class1_`lhs' ==.
                }
                local mix "prob_class1_`lhs'"
                set matadebug on
                mata:
em("`lhs'","`rhs'","`mix'","`by'","`touse'",`k',"`cluster'")
                ereturn post
                ereturn display
end

version 11
mata:
mata set matalnum on
mata set matastrict on


real scalar function estep(real rowvector b1, real rowvector b2, real
scalar s1, real scalar s2, string scalar lhs, string scalar rhs,
string scalar mix, string scalar by, string scalar touse) {
            real scalar logli_l, logli_h, sup
            real vector y, id, pi, old_pi, pi_i
            real matrix X
            st_view(id,.,by,touse)
            info = panelsetup(id, 1)
            st_view(X,.,rhs,touse)
            X = (X,J(rows(X),1,1))
            st_view(y,.,lhs,touse)
            st_view(pi,.,mix,touse)
            old_pi = pi
            index_hi = normalden(y:-X*b1', 0, s1)
            index_li = normalden(y:-X*b2', 0, s2)
            for (i=1;i<=rows(info);i++) {
                index_l = panelsubmatrix(index_li,i,info)
                index_h = panelsubmatrix(index_hi,i,info)
                panelsubview(pi_i,pi,i,info)
                logli_l = quadcolsum(log(index_l))
                logli_h = quadcolsum(log(index_h))
                if (logli_l-logli_h > 45) pi_i[.] = J(rows(pi_i),1,0)
                else    if (logli_l-logli_h < -45) pi_i[.] = J(rows(pi_i),1,1)
                        else pi_i[.] =
J(rows(pi_i),1,1/(1+exp(logli_l-logli_h)))
            }
            if (hasmissing(pi)) printf("The minimum of the first index
is %9.6g, the maximum is %9.6g (s.e. %9.6g). The minimum of the second
index is %9.6g, the maximum is %9.6g (s.e.
%9.6g).",min(y:-X*b1'),max(y:-X*b1'),s1,min(y:-X*b2'),max(y:-X*b2'),s2)
            if (hasmissing(pi)) _error(3351, "E-step produced missing value.")
            sup = colmaxabs(pi-old_pi)
            return(sup)
}
mata mosave estep(), replace // dir(PERSONAL)

void function mstep(transmorphic M,
                    real scalar todo,
                    real rowvector b, fv, S, H)
{
                      real colvector  p1, p2, s1, s2, h11_1, h11_2,
h22_1, h22_2, h12_1, h12_2, h1121, h2212, h1112, h2122, h1122, h1221
                      real matrix H11_1, H11_2, H22_1, H22_2, H12_1,
H12_2, H1121, H2212, H1112, H2122, H1122, H1221
                      real colvector  y
                      real colvector mix
                      real scalar t1, t14, t15, t52
                      real colvector t2, t3, t5, t6, t8, t9, t13, t16,
t19, t20, t21, t22, t23, t24, t26, t29, t30, t33, t34, t35, t38, t39,
t42, t43, t44, t46, t51, t55, t60, t65, t66, t70, t80, t82, t83, t84,
t88, t89, t90, t92, t93, t94, t97, t101, t104, t105, t107, t112, t117,
t120, t132, t134, t145, t147
//                      real scalar ll

                      p1 = moptimize_util_xb(M, b, 1)
                      p2 = moptimize_util_xb(M, b, 2)
                      s1 = moptimize_util_xb(M, b, 3)
                      s2 = moptimize_util_xb(M, b, 4)
                      y = moptimize_util_depvar(M, 1)

st_view(mix,.,moptimize_util_userinfo(M,1),moptimize_util_userinfo(M,2))

                        if (hasmissing(mix)) _error(3351, "mix has
missing value.")
                        if (hasmissing(y)) _error(3351, "y has missing value.")
                        if (hasmissing(p1)) _error(3351, "p1 has
missing value.")
                        if (hasmissing(p2)) _error(3351, "p2 has
missing value.")
                        if (hasmissing(s1)) _error(3351, "s1 has
missing value.")
                        if (hasmissing(s2)) _error(3351, "s2 has
missing value.")
                                t1 = sqrt(2)
                                t2 = t1:*mix
                                t3 = exp(s1)
//                                t4 = t3:^2
                                t5 = t3:^(-2)
                                t6 = t5:/t3
                                t8 = y:-p1
                                t9 = t8:^2
                                t13 = exp(-1*t5:*t9:/2)
                                t14 = sqrt(pi())
                                t15 = 1/t14
                                t16 = t15:*t13
                                t19 = t16:*t2:/t3
                                t20 = -1*(mix:-1)
                                t21 = t1:*t20
                                t22 = y:-p2
                                t23 = t22:^2
                                t24 = exp(s2)
//                                t25 = t24:^2
                                t26 = t24^(-2)
                                t29 = exp(-1*t26:*t23:/2)
                                t30 = t15:*t29
                                t33 = 1:/t24:*t30:*t21
                                t34 = t19:+t33
                                t35 = 2*t34:^(-1)
                                t38 = t9:*t2
                                t39 = t3:^4
                                t42 = t13:/t39:/t3
                                t43 = t35:*t15
                                t44 = t43:*t42
                                t46 = mix:^2
                                t51 = t13:^2
                                t52 = 1/pi()
                                t55 = 4*t34:^(-2)
                                t60 = t13:*t6
                                t65 = t26:/t24
                                t66 = t65:*t22
                                t70 = t29:*t66:*t20:*t55:*t52:*t60:*t8:*mix:/2
                                t80 = t6:*t8:*t2
                                t82 = t15:*t60:*t38
                                t83 = t82:-t19
                                t84 = t83:*t55:/2
                                t88 =
-3/2*t43:*t60:*t8:*t2:+t44:*t9:*t8:*t2:/2:-t84:*t16:*t80:/2
                                t89 = t23:*t21
                                t90 = t29:*t65
                                t92 = t15:*t90:*t89
                                t93 = t92:-t33
                                t94 = t93:*t55:/2
                                t97 = t94:*t16:*t80:/2
                                t101 = t24:^4
                                t104 = t29:/t101:/t24
                                t105 = t43:*t104
                                t107 = t20:^2
                                t112 = t29:^2
                                t117 = t66:*t21
                                t120 = t84:*t30:*t117:/2
                                t132 =
-3/2*t43:*t90:*t22:*t21:+t105:*t23:*t22:*t21:/2:-t94:*t30:*t117:/2
                                t134 = t9:^2
                                t145 = t93:*t84:/2
                                t147 = t23:^2

                      fv = ln(mix:*(normalden(t8, 0,
t3)-normalden(t22, 0, t24))+normalden(t22, 0, t24))

//                      ll = moptimize_util_sum(M,fv)

                      if (todo>=1) {

                        S   = (t43:*t60:*t8:*t2:/2,
t43:*t90:*t22:*t21:/2, t35:*(t15:*t60:*t9:*t2:-t19):/2,
t35:*(t15:*t90:*t23:*t21:-t33):/2)

                        if (todo==2) {

                                h11_1 =
-1*t35:*t16:*t6:*t2:/2:+t44:*t38:/2:-t55:*t52:*t51:/t39:*t5:*t9:*t46:/2
                                h1121 = -1*t70
                                h12_1 = t88
                                h1122 = -1*t97
                                h11_2 =
-1*t35:*t30:*t65:*t21:/2:+t105:*t89:/2:-t55:*t52:*t112:/t101:*t26:*t23:*t107:/2
                                h1221 = -1*t120
                                h12_2 = t132
                                h22_1 =
t35:*(-2*t82:+t15:*t42:*t134:*t2:/2:+t19:/2):-t55:*t83:^2:/4
                                h2212 = -1*t145
                                h22_2 =
t35:*(-2*t92:+t15:*t104:*t147:*t21:/2:+t33:/2):-t55:*t93:^2:/4

                                H11_1 = moptimize_util_matsum(M, 1,1, h11_1, 1)
                                H22_1 = moptimize_util_matsum(M, 3,3, h22_1, 1)
                                H12_1 = moptimize_util_matsum(M, 1,3, h12_1, 1)
                                H11_2 = moptimize_util_matsum(M, 2,2, h11_2, 1)
                                H22_2 = moptimize_util_matsum(M, 4,4, h22_2, 1)
                                H12_2 = moptimize_util_matsum(M, 2,4, h12_2, 1)
                                H1121 = moptimize_util_matsum(M, 1,2, h1121, 1)
                                H2212 = moptimize_util_matsum(M, 3,4, h2212, 1)
                                H1122 = moptimize_util_matsum(M, 1,4, h1122, 1)
                                H1221 = moptimize_util_matsum(M, 2,3, h1221, 1)

                                H     = (   H11_1,   H1121,   H12_1,   H1122 \
                                            H1121',   H11_2,   H1221,   H12_2 \
                                            H12_1',   H1221',   H22_1,   H2212 \
                                            H1122',   H12_2',   H2212',   H22_2)
                        }
                }
        }
mata mosave mstep(), replace // dir(PERSONAL)

void function em(string scalar lhs, string scalar rhs, string scalar
mix, string scalar by, string scalar touse, real scalar k, string
scalar cluster) {
    transmorphic M
    real vector b1o, b1n, b2o, b2n, s1o, s1n, s2o, s2n
    real scalar d, sup
    b1o = J(1,k+1,1)
    b2o = J(1,k+1,1)
    s1o = 0.5*colmaxabs(st_data(.,lhs,touse))
    s2o = 0.5*colmaxabs(st_data(.,lhs,touse))
    d = 1
    sup = 1
while (d>epsilon(10^8)) {
// M-step
    M = moptimize_init()
    moptimize_init_evaluator(M, &mstep())
    moptimize_init_evaluatortype(M,"lf2")
    moptimize_init_touse(M, touse)
    moptimize_init_ndepvars(M, 1)
    moptimize_init_depvar(M, 1, lhs)
    moptimize_init_eq_n(M, 4)
    moptimize_init_eq_indepvars(M, 1, rhs)
    moptimize_init_eq_cons(M, 1, "on")
    moptimize_init_eq_indepvars(M, 2, rhs)
    moptimize_init_eq_cons(M, 2, "on")
    moptimize_init_eq_indepvars(M, 3, "")
    moptimize_init_eq_cons(M, 3, "on")
    moptimize_init_eq_indepvars(M, 4, "")
    moptimize_init_eq_cons(M, 4, "on")
    moptimize_init_eq_name(M, 1, "class1")
    moptimize_init_eq_name(M, 2, "class2")
    moptimize_init_eq_name(M, 3, "se1")
    moptimize_init_eq_name(M, 4, "se2")
    moptimize_init_cluster(M, cluster)
    moptimize_init_eq_coefs(M, 1, b1o)
    moptimize_init_eq_coefs(M, 2, b2o)
    moptimize_init_eq_coefs(M, 3, s1o)
    moptimize_init_eq_coefs(M, 4, s2o)
    moptimize_init_search(M, "on")
    moptimize_init_search_random(M, "on")
//    moptimize_init_search_bounds(M, 3, (0,.))
//    moptimize_init_search_bounds(M, 4, (0,.))
    moptimize_init_valueid(M, "LogL with given mixing")
    moptimize_init_nuserinfo(M, 2)
    moptimize_init_userinfo(M, 1, mix)
    moptimize_init_userinfo(M, 2, touse)
    moptimize_init_singularHmethod(M, "hybrid")
//    moptimize_init_tracelevel(M,"hessian")
//    moptimize_init_conv_maxiter(M, 2)
//    moptimize_query(M)
    _moptimize(M)
//    printf("The moptimize return code is
%5.4g",moptimize_result_returncode(M))
    moptimize_result_post(M)
//    moptimize_result_display()
    b1n = moptimize_result_eq_coefs(M,1)
    b2n = moptimize_result_eq_coefs(M,2)
    s1n = moptimize_result_eq_coefs(M,3)
    s2n = moptimize_result_eq_coefs(M,4)
    d = rowmaxabs((b1n-b1o,b2n-b2o,s1n-s1o,s2n-s2o, sup))
    b1o = b1n
    b2o = b2n
    s1o = s1n
    s2o = s2n
// E-step
    sup = estep(b1o,b2o,s1o,s2o,lhs,rhs,mix,by,touse)
}
}
mata mosave em(), replace // dir(PERSONAL)

end
exit
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/