Thank you for all your help today.
I have another quick question. I saved my first ever user-written ado
file in my personal folder. The proof is the output:
. personal dir
your personal ado-directory is ~/ado/personal/
emlatentclass.ado
But when I would run it (under Stata 11), I get:
. emlatentclass d_fdeic treat, by(taxproid07)
last estimates not found
(error occurred while loading emlatentclass.ado)
r(301);
I tried to catch basic syntax errors, to no avail so far. I'd be
grateful for any quick guesses what are the usual suspects, or you
might even skim my draft of code below. Thank you!
Laszlo
program emlatentclass, eclass
version 11
syntax varlist [if] [in], BY(varname)
marksample touse
markout touse `varlist' `by'
tokenize `varlist'
local lhs `1'
macro shift
local rhs `*'
local k: word count `rhs'
qui sum touse
local nobs = r(sum)
gen prob_class1 = 0.5
local mix "prob_class1"
set matadebug on
mata:
mata set matalnum on
void 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) {
real scalar logli_l, logli_h
real vector y
real matrix X
st_view(id,.,by)
info = panelsetup(id, 1)
st_view(X,.,rhs)
st_view(y,.,lhs)
st_view(pi,.,mix)
index_hi = normalden(y:-X*b1', 0, sqrt(s1))
index_li = normalden(y:-X*b2', 0, sqrt(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))
pi_i[.] =
J(rows(nu_i),1,exp(logli_h)/(exp(logli_h)+exp(logli_l)))
}
}
mata mosave estep, replace
function mstep(transmorphic M,
real rowvector b,
real colvector lnf)
{
real colvector p1, p2
real scalar s1, s2
real colvector y1
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))
lnf = ln(mix:*normalden(y:-p1, 0,
sqrt(s1))+(1:-mix):*normalden(y:-p2, 0, sqrt(s2)))
}
mata mosave mstep, replace
b1 = J(1,`k',1)
b2 = J(1,`k',1)
s1 = 1
s2 = 1
d = 1
while d>10^(-3)) {
// E-step
estep(b1,b2,s1,s2,"`lhs'","`rhs'","`mix'","`by'")
// M-step
M = moptimize_init()
moptimize_init_evaluator(M, &mstep())
moptimize_init_evaluatortype(M,"lf")
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_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_by(M, "`by'")
moptimize_init_eq_coefs(M, 1, b1)
moptimize_init_eq_coefs(M, 2, b2)
moptimize_init_eq_coefs(M, 3, s1)
moptimize_init_eq_coefs(M, 4, s2)
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, 1)
moptimize_init_userinfo(M, 1, "`mix'")
_moptimize(M)
moptimize_result_post(M)
b1n = moptimize_result_coefs(M,1)
b2n = moptimize_result_coefs(M,2)
s1n = moptimize_result_coefs(M,3)
s2n = moptimize_result_coefs(M,4)
d = sum(abs((b1n-b1,b2n-b2,s1n-s1,s2n-s2)))
b1 = b1n
b2 = b2n
s1 = s1n
s2 = s2n
}
end
// ereturn post b V, esample(touse) o(`nobs')
ereturn display
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/
