Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: Using multiple equations with Optimize()


From   "Andrew Lynch" <aalhc5@mizzou.edu>
To   <statalist@hsphsun2.harvard.edu>
Subject   Re: st: Using multiple equations with Optimize()
Date   Mon, 12 Nov 2007 20:47:30 -0600

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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index