Solving A System of Equations

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  -

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?

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

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