Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
st: too many variables from simple Mata iteration?
From
László Sándor <[email protected]>
To
[email protected]
Subject
st: too many variables from simple Mata iteration?
Date
Sun, 7 Mar 2010 14:00:05 -0500
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/