[R] Solving A System of Equations
spencer.graves at pdf.com
Wed Apr 9 17:57:17 CEST 2003
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
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
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
More information about the R-help