László Sándor <[email protected]>:
I haven't looked through your program, or any prior posts this week,
but the Mata stuff needs to go outside the program stuff. I.e.
program x
stuff
end
mata:
stuff
end
2009/11/6 László Sándor <[email protected]>:
> 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/
