Re: st: Using multiple equations with Optimize()

["Andrew Lynch" <aalhc5@mizzou.edu>]

Thanks Jeff!

I should have been more thorough, as I had initialized the variables (I 
didn't include the full do file, as it is kind of long).  Though I hadn't 
tested the do loop to make sure it worked.  I tested just the loop (and it 
works) but still seem to have a problem.  Just to make sure, can I include a 
do loop in an optimize procedure evaluation function?   If I do will it 
re-evaluate the do loop at each update of the parameters?  I have tried 
several simplistic optimization simulations just to make sure it works, and 
even though the do loop I include works on its own I get the same error when 
I run the optimization (r3201 vector required).  The example below is 
trivial, but I am just trying to find an example that works with a do loop.


clear all
mata:
mata clear

/* Data */
    y=invnormal(uniform(252,1))
    xt=invnormal(uniform(252,1))

/* Setup */
   n=rows(y)
   c=J(n,1,.)
   x=(c,xt)

/* Test Optimization */
    void maxlik(todo, b, x, y, v, g, h)
    {
         b0=b[1]
         b1=b[2]
         i=1
       do {
          x[i,1]=i
          i=i+1
       } while (i<(n+1))

     d = y-x*(b0,b1)'
     v = d'*d
    }

S = optimize_init()
optimize_init_evaluator(S, &maxlik())
optimize_init_evaluatortype(S, "d0")
optimize_init_params(S, J(1,2,0))
optimize_init_which(S, "min")
optimize_init_argument(S, 1, x)
optimize_init_argument(S, 2, y)
b=optimize(S)

b

end



----- Original Message ----- 
From: "Jeff Pitblado, StataCorp LP" <jpitblado@stata.com>
To: <statalist@hsphsun2.harvard.edu>
Sent: Monday, November 12, 2007 4:40 PM
Subject: Re: st: Using multiple equations with Optimize()



> Andrew Lynch <aalhc5@mizzou.edu> is using -optimize()- fit a model:
>
> > I have run into a bit of a dilemma with programing a Kalman filter with 
> > the
> > optimize() procedure.  I have the full documentation and have not been 
> > able
> > to clarify where exactly the optimize() command looks when it optimizes. 
> > I
> > am minimizing the process below.  My goal is to minimize the Log 
> > Likelihood
> > function, but while I have observed data yt and xt, a and at are 
> > unobserved
> > and my goal is to generate them.  To do so, I need Stata to rerun the do
> > loop every time it updates one of the four parameters (b1, b2, q and h). 
> > I
> > say all that to ask the question, "Does optimize() only look to the v=
> > equation to optimize with the parameters or does it look to the entire
> > evaluation and move through it in some ordered fashion (thereby honoring 
> > the
> > do loop)?"  I have been unable to get this to run, and I am unsure 
> > whether
> > it is my dumb coding error at some point or just that optimize() won't do
> > what I want it to do.
> >
> > Any enlightenment would be greatly appreciated!
> -optimize()- will call the evaluator function each time it changes the
> parameter values.  The function value returned by the evaluator is used to
> compute numerical derivatives (in the case of type d0 and v0 evaluators) 
> in
> addition to the Newton-Raphson step and determining convergence.
>
> Unfortunately, Andrew's evaluator function is referring to undefined 
> objects.
>
> For example, the -n- object (which appears to be a -real scalar-) is not
> defined within -maxlik()-.  It appears that -n- is probably the number of
> elements in -xt- and -yt-; say
>
> n = rows(xt)
>
> The -do- loop in -maxlik()- is building intermediate calculations that
> ultimately go into -L-, which appears to be a -real vector-.  -L- must be
> initialized before it can be "filled-in":
>
> L = J(n,1,.)
>
> Finally, the following vectors need to be initialized before their 
> elements
> can be referenced:
>
> at
> ft
> pt
>
> a
> f
> p
> y
>
> Maybe some of them should be initialized to the column vector of zeros, 
> such
> as
>
> at = J(n+1,1,0)
>
> My best suggestion to Andrew is to work with -maxlik()- directly before 
> using
> it with -optimize()-.  Build a do-file that verifies that -maxlik()- 
> returns
> valid -v- values given predetermined parameters values with his -xt-/-yt-
> data.  Andrew will have to address the issues I mention above, so the
> following do-file could be used as a starting point; it runs without a
> run-time error, so all Andrew needs to do now is modify -maxlik()- (from
> within check.do) until it produces the correct values given -b-, -xt-, and
> -yt-.
>
> ***** BEGIN: check.do
> clear all
>
> // NOTE: I'm just generating some random values here, but you should use
> // inputs that will allow you to verify that -mylik()- is producing valid
> // results.
> set seed 1234
>
> mata:
>
> xt = invnormal(uniform(100,1))
> yt = invnormal(uniform(100,1))
>
> // current parameter values
> b = J(1,4,0)
>
> void maxlik(todo, b, xt, yt, v, g, h)
> {
>        b1      = b[1]
>        b2      = b[2]
>        q       = b[3]
>        // NOTE: this seems like a scale parameter, so you might consider
>        // treating it as if it were from a log space (as I've done here) 
> to
>        // ensure that -h- is positive
>        h       = exp(b[4])
>
>        n       = rows(xt)
>        L       = J(n,1,0)
>
>        // NOTE: figure out how the following are supposed to be 
> initialized
>        at      = J(n+1, 1, 0)
>        ft      = J(n+1, 1, 0)
>        pt      = J(n+1, 1, 0)
>        a       = J(n, 1, 0)
>        f       = J(n, 1, 0)
>        p       = J(n, 1, 0)
>        y       = J(n, 1, 0)
>
>        i       = 1
>        do {
>                a[i]=b1*at[i] + b2*xt[i]
>                p[i]=(b1)^2*pt[i] + q
>                p[i]=(b1)^2*pt[i] + q
>                y[i]=a[i]
>                f[i]=ft[i] + h
>                at[i+1]=a[i] + ft[i+1]*p[i]*(yt[i]-y[i])
>                pt[i+1]=p[i] - ft[i+1]*p[i]*p[i]
>                L[i]=(.5)/f[i] + (.5)*(yt[i]-y[i])*(yt[i]-y[i])/f[i]
>                i=i+1
>        } while(i<(n+1))
>
>        v=sum(L)
> }
>
> // compute the function value
> maxlik(., b, xt, yt, v=., ., .)
>
> // display the function value, if you know what the value should be, then
> // maybe show the -reldif()- between what you got and what you expect
> v
>
> end
> ***** END: check.do
>
> --Jeff
> jpitblado@stata.com
>
> For reference, here is the Mata code that Andrew submitted:
>
>
> >  /* Optimization Procedure */
> >
> >  void maxlik(todo, b, xt, yt, v, g, h)
> >  {
> >       b1=b[1]
> >       b2=b[2]
> >       q=b[3]
> >       h=b[4]
> >       i=1
> >
> >       do {
> >            a[i]=b1*at[i] + b2*xt[i]
> >            p[i]=(b1)^2*pt[i] + q
> >            y[i]=a[i]
> >            f[i]=ft[i] + h
> >            at[i+1]=a[i] + ft[i+1]*p[i]*(yt[i]-y[i])
> >            pt[i+1]=p[i] - ft[i+1]*p[i]*p[i]
> >            L[i]=(.5)/f[i] + (.5)*(yt[i]-y[i])*(yt[i]-y[i])/f[i]
> >            i=i+1
> >
> >       } while(i<(n+1))
> >
> >   v=sum(L)
> >  }
> >
> >  S = optimize_init()
> >  optimize_init_evaluator(S, &maxlik())
> >  optimize_init_evaluatortype(S, "d0")
> >  optimize_init_params(S, (0,0,0,0))
> >  optimize_init_which(S, "min")
> >  optimize_init_argument(S, 1, xt)
> >  optimize_init_argument(S, 2, yt)
> >  b=optimize(S)
> >
> > *
> > *   For searches and help try:
> > *   http://www.stata.com/support/faqs/res/findit.html
> > *   http://www.stata.com/support/statalist/faq
> > *   http://www.ats.ucla.edu/stat/stata/
> *
> *   For searches and help try:
> *   http://www.stata.com/support/faqs/res/findit.html
> *   http://www.stata.com/support/statalist/faq
> *   http://www.ats.ucla.edu/stat/stata/ 
*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/