[Rd] Using the nls package
Douglas Bates
bates@stat.wisc.edu
29 Jul 2000 18:59:53 -0500
Telford Tendys <telford@progsoc.uts.edu.au> writes:
> I'm a bit confused about the nls package, I'm trying to use it for curve fitting.
>
> First off, the documentation for nls says ``see `nlsControl' for the
> names of the settable control values'' -- this is wrong, it should be
> nls.control (minor point but had me confused for a moment).
Thanks for pointing that out. We will correct that for the next
release.
> Now I'll try something very simple (maybe too simple):
> ----------------------------------------------------------------------
> > q <- data.frame( x=1:10 )
> > q$y <- eval( expression( 4 + 5 * x ), q )
> > library( nls )
> > w <- nls( y ~ P1 + P2 * x, q, start=list( P1=1, P2=1 ))
> Error in nls(y ~ P1 + P2 * x, q, start = list(P1 = 1, P2 = 1)) :
> maximum number of iterations exceeded
> > w <- nls( y ~ P1 + P2 * x, q, start=list( P1=1, P2=1 ), trace=T, control=nls.control( maxiter=10 ))
> 7570 : 1 1
> 7.888609e-29 : 4 5
> 0 : 4 5
> 0 : 4 5
> 0 : 4 5
> 0 : 4 5
> 0 : 4 5
> 0 : 4 5
> 0 : 4 5
> 0 : 4 5
> 0 : 4 5
> ----------------------------------------------------------------------
>
> There is something I really don't understand here -- it seems to
> go right to the answer but then cannot detect that it should
> terminate the search. The default tolerance of 1e-5 seems reasonable.
> Same happens with a slightly harder problem:
The convergence criterion for nls is the relative offset orthogonality
convergence criterion as described in
@Article{bate:watt:1981:tech,
journal = "Technometrics",
volume = 23,
pages = "179--183",
keywords = "Regression",
author = "Douglas M. Bates and Donald G. Watts",
title = "A Relative Offset Orthogonality Convergence Criterion
for Nonlinear Least Squares",
year = 1981
}
This criterion fails for artificial examples where the residual
sum-of-squares can be made exactly zero.
> ----------------------------------------------------------------------
> > q$y <- eval( expression( 4 + 5 / x + 7 * x ), q )> w <- nls( y ~ P1 + P2 / x + P3 * x, q, start=list( P1=1, P2=1, P3=1 ), trace=T, control=nls.control( maxiter=10 ))
> 16505.09 : 1 1 1
> 3.578718e-14 : 4 5 7
> 2.650573e-28 : 4 5 7
> 1.262177e-29 : 4 5 7
> 0 : 4 5 7
> 0 : 4 5 7
> 0 : 4 5 7
> 0 : 4 5 7
> 0 : 4 5 7
> 0 : 4 5 7
> 0 : 4 5 7
> ----------------------------------------------------------------------
>
> It knows what the answer should be but does not know that it has the
> answer, try injecting some noise in the sample points:
What is unusual here is that you get an exact zero for the residual
sum-of-squares. Often there would be a small, non-zero RSS even at
the parameter values used to simulate the response.
> ----------------------------------------------------------------------
> > q$y <- eval( expression( 4 + 5 / x + 7 * x + runif( 10 )), q )
> > w <- nls( y ~ P1 + P2 / x + P3 * x, q, start=list( P1=1, P2=1, P3=1 ), trace=T, control=nls.control( maxiter=10 ))
> 16882.09 : 1 1 1
> 0.4923209 : 4.796474 5.290280 6.947685
> ----------------------------------------------------------------------
>
> OK, now it stops right away with an approximate solution...
> must be alergic to getting the right answer!
> I guess with all that noise, this answer is probably the best possible.
>
> Well, I'm sure that isn't such a big bug to fix, probably
> a test for >0 when it should be >=0 or something. Maybe just
> detect exceptionally small sum of squares (i.e. zero) and bail
> out right there. Do statisticians have an instinct to disbelieve that
> their sample data could ever fit the model?
The problem with putting in a check for an exceptionally small sum of
squares is deciding when it is exceptionally small. There is no
absolute scale on the residual sum of squares.
We could perhaps detect an exact zero but even that would not catch
all the cases where "exact" values were simulated because generally
they will give a small but non-zero sum of squares.
Since such problems will not occur in practice but only in artificial
examples, we have not bothered detecting them.
> Now a bigger problem (for me), I want to be able to use complex
> numbers in the calculation:
>
> ----------------------------------------------------------------------
> > q$y <- eval( expression( 4 + x*(0+1i)), q )
> > w <- nls( y ~ P1 + P2 * x, q, start=list( P1=1, P2=1 ), trace=T, control=nls.control( maxiter=10 ))
> Error in model.frame(formula, rownames, variables, varnames, extras, extranames, :
> invalid variable type
> > traceback()
> [1] "model.frame.default(formula = ~y + x, data = q)"
> [2] "model.frame(formula = ~y + x, data = q)"
> [3] "eval(expr, envir, enclos)"
> [4] "eval(mf, sys.frame(sys.parent()))"
> [5] "as.list(eval(mf, sys.frame(sys.parent())))"
> [6] "nls(y ~ P1 + P2 * x, q, start = list(P1 = 1, P2 = 1), trace = T, "
> [7] " control = nls.control(maxiter = 10))"
> > model.frame.default( y ~ x, q )
> Error in model.frame(formula, rownames, variables, varnames, extras, extranames, :
> invalid variable type
> ----------------------------------------------------------------------
>
> Is there any chance of getting this to work? Should I split up the
> complex numbers into components and write my equations are equivalent
> expressions in reals (ugly but I could do it if I have to)?
That's the best short-term solution. In general, how do you want to
define the residual sum of squares for a complex-valued response?
Should it be the sum of the modulus of the residuals?
We could put complex-valued nonlinear least squares on the wishlist
but I would have to warn you that there are many other items in front
of it.
> Seems like something in the .Internal(model.frame(...)) is barfing
> at a complex number. Is that necessary or just an accident? Surely
> R has enough complex number capabilities to handle them here too...
Easiest thing is to test it.
> data(cars)
> model.frame(dist ~ speed, data = cars)
dist speed
1 2 4
...
49 120 24
50 85 25
> cars$speed <- cars$speed + rnorm(50)*(0+1i)
> is.complex(cars$speed)
[1] TRUE
> model.frame(dist ~ speed, data = cars)
Error in model.frame(formula, rownames, variables, varnames, extras, extranames, :
invalid variable type
So, yes, the model.frame function does not allow complex data.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._