[R] Parameter estimation in nls
spencer.graves at pdf.com
Tue Nov 25 16:43:32 CET 2003
Since x <- 1:26 and your y's are all positive, your model,
ignoring the error term, is mathematically isomorphic to the following:
> x <- 1:26
> (fit <- lm(y~x+log(x)))
lm(formula = y ~ x + log(x))
(Intercept) x log(x)
35802074 -371008 -8222922
With reasonable starting values, I would expect "a" to converge to
roughly exp(35802074), "k" to (-8222922), and "b" to exp(-371008). With
values of these magnitudes for "a" and "b", the "nls" optimization is
highly ill conditioned.
What do these numbers represent? By using "nls" you are assuming
implicitly the following:
y = a*x^k*b^x + e, where the e's are independent normal errors
with mean 0 and constant variance.
Meanwhile, the linear model I fit above assumes a different noise
log(y) = log(a) + k*log(x) + x*log(b) + e, where these e's are
also independent normal, mean 0, constant variance.
If you have no preference for one noise model over the other, I
suggest you use the linear model I estimated, declare victory and worry
about something else. If you insist on estimating the multiplicative
model, you should start by dividing y by some number like 1e6 or 1e7 and
consider reparameterizing the problem if that is not adequate. Have you
consulted a good book on nonlinear regression? The two references cited
in "?nls" are both excellent:
Bates, D.M. and Watts, D.G. (1988) _Nonlinear Regression Analysis
and Its Applications_, Wiley
Bates, D. M. and Chambers, J. M. (1992) _Nonlinear models._
Chapter 10 of _Statistical Models in S_ eds J. M. Chambers and T.
J. Hastie, Wadsworth & Brooks/Cole.
hope this helps. spencer graves
Dr Andrew Wilson wrote:
>I am trying to fit a rank-frequency distribution with 3 unknowns (a, b
>and k) to a set of data.
>This is my data set:
>y <- c(37047647,27083970,23944887,22536157,20133224,
>and this is the fit I'm trying to do:
>nlsfit <- nls(y ~ a * x^k * b^x, start=list(a=5,k=1,b=3))
>(It's a Yule distribution.)
>However, I keep getting:
>"Error in nls(y ~ a * x^k * b^x, start = list(a = 5, k = 1, b = 3)) :
>I guess this has something to do with the parameter start values.
>I was wondering, is there a fully automated way of estimating parameters
>which doesn't need start values close to the final estimates? I know
>other programs do it, so is it possible in R?
>R-help at stat.math.ethz.ch mailing list
More information about the R-help