[R] Solving A System of Equations

[Spencer Graves]

The best method I know is as follows:

1.  Write a function to be minimized, e.g., sum of squares of errors in 
the equations.  If your functions are continuously differentiable, then 
sum of squares is better than sum of absolute values, because the sum of 
squares it is continuous at the minimum while the sum of absolute values 
is not.

2.  If there is any question that there might be multiple local minima, 
then I would implement one of Prof. Blackwell's suggestions, namely 
testing it over an appropriate grid of points.  With only one unknown, I 
make a plot.  With two, I make a contour plot.  With three or more, I 
might try some kind of grid or Monte Carlo.

3.  With some appropriate starting value(s), I then pass the function to 
"optim".  If I want confidence intervals with nonlinear least squares, I 
may pass the output of optim to "nls".  If my objective function is a 
log(likelihood), "optim" will output the Hessian, which is the negative 
of the observed information, whose invers in the approximate covariance 
matrix.

Comments?
Spencer Graves

Thomas W Blackwell wrote:
> Maybe there are other R tools, but what I do, repeatedly, is
> to work out first and second derivatives w.r.t. x on paper.
> Then I execute one step of Newton's method in R, line by line
> at the command line, however many lines it takes.  Then use
> up-arrow to repeat those command lines for maybe five or six
> iterations of Newton's method.  Quick, clear and cheap.
> 
> The point is:  start with a four dimensional grid of starting
> values for x,y,z,zz.  This can be a four-dimensional array in
> which each variable is initially constant wrt the other three.
> In subsequent iterations, particular x,y,z,zz quadruples will
> either converge toward a solution or diverge, and you can use
> standard tests to distinguish between the two.
> 
> After you've done this once or twice by hand, at the command
> line, then you'll be ready to code up a single Newton step
> as a function, then call it an appropriate (fixed) number of
> times inside an outer function.  I wouldn't worry about trying
> to be adaptive and economize on computer cycles.  All of that
> stuff from the older numerical analysis literature is misplaced
> effort, in my view, in this day and age.
> 
> Other users may know of much fancier methods.
> 
> -  tom blackwell  -  u michigan medical school  -  ann arbor  -
> 
> On Tue, 8 Apr 2003, David Tyler wrote:
> 
> 
>>I'm trying to solve a system of 3 equations
>>as part of a sub-routine in R, ie
>>first eqn a/x-b*sqrtx+c=log(1/dx+1/e(sqrtx);
>>snd eqn (f*y)/z-g/y-h=-log(2/x+(z/y)/(i*x)
>>and third eqn is of the form zz=x/(j-k(z/y)
>>where a..k inclusive are constants, x,y,z
>>and zz are inputs.
>>
>>How can this be done in R?
> 
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://www.stat.math.ethz.ch/mailman/listinfo/r-help