Re: st: passing extra information to function evaluator program nl

[Stas Kolenikov <skolenik@gmail.com>]

You might find it more convenient to use -moremata-'s (kudos to Ben
Jann) -mm_root()- or -mm_nrroot()- functions. They don't need any
awkward data management; you still need to code your evaluator
functions, but you can pass anything to them in a very natural manner.

ssc install moremata, repalce

mata:

real scalar function mypolynom(real scalar x, real vector a) {
   K = cols(a)
   f = 0
   for(k=1;k<=K;k++) {
      f = f*x + a[k]
   }
  return( f )
}

A = (1, -2, -1 \ 2, 5, -4 \ 3, -8 ,-1 )

for(i=1;i<=3;i++) {
   rc = mm_root( x=., &mypolynom(), 0, 10, 0, 1000, A[i,.])
   if (!rc) x
}

end
// of mata

Note that -mm_root()- is pretty dumb, and does expect the user to
provide smart endpoints of an interval in which the root is known to
be located (so that the values of the function are of opposite sign).
-mm_nrroot()- is somewhat smarter, but I couldn't make it run, for
some reason (Mata version conflict of some kind?)

On Tue, Mar 23, 2010 at 11:14 AM, Ben Zipperer <benzipperer@gmail.com> wrote:
> I am interested in using Stata to find a real roots to a single
> polynomial several times, where the coefficients vary each time.  I
> basically want to apply Newton's method or something else to vectors
> of coefficients. I originally planned to restrict the solver at
> http://www.stata.com/support/faqs/lang/nl.html to a single polynomial
> and then iterate it, changing the coefficients at each iteration.
> However, I can't figure out how to pass the coefficients to the
> function evaluator program. How do I do this? Do I modify the syntax
> statement of the program? Or should I try a different strategy?
>
> For example, let's say the rows of the following matrix form the
> coefficients for a quadratic polynomial
>
> mat def M = (1,2,3\ 4,5,6\ 7,8,9)
>
> Then the nl program evaluator program is something like
>
> capture program drop nlpolysolve
> program nlpolysolve
>
>        syntax varlist(min=1 max=1) [if], at(name)
>
>        tempname x z
>        scalar `x' = `at'[1, 1]
>        scalar `z' = `at'[1, 2]
>
>        tempvar yh
>        // the real polynomial
>        gen double `yh' = a_`i' * x^2 + b_`i' * `x' + c_`i' + 1 in 1
>
>        replace `yh' = 1 - `z' in 2
>        replace `varlist' = `yh'
>
> end
>
> And my iterations are something like
>
> forvalues i=1/3 {
>  preserve
>  clear
>  set obs 2
>  gen y = 0
>  replace y = 1 in 1
>
>  local a_`i' = M[`i',1]
>  local b_`i' = M[`i',2]
>  local c_`i' = M[`i',3]
>
>  nl polysolve @ y, parameters(A B) initial(A 0.5 B 1)
> stuffiwanttopass(a_`i' b_`i' c_i')
>
>  restore
> }
>
> How do I pass the entries of the matrix to the nlpolysolve program above?
>
> thanks,
> Ben
> *
> *   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/
>



-- 
Stas Kolenikov, also found at http://stas.kolenikov.name
Small print: I use this email account for mailing lists only.

*
*   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/