Solving A System of Equations
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