[R] select nls starting values from list was: nls does not accept start values
Petr PIKAL
petr.pikal at precheza.cz
Tue Mar 3 11:55:59 CET 2009
Hallo
Here is my "final" solution, but i still wonder if it could not be
improved somehow.
fit <- modely(fol, data.frame(x=x, y=y1)
fit
or
fit <- modely(fol, data.frame(x=x, y=y2)
fit
shows that it arrives to some results. However selection starting values
as means of response variable seems to me suboptimal.
Is it possible to have some list of functions which will for each item in
formula list compute starting values for nls? Or is it this mean guess as
good as any other?
Petr
petr.pikal at precheza.cz
# formulas
fol <- list(y~1/(a-x), y~a*x^2+b*log(x))
set.seed(111)
x<-1:10
y1 <- 1/(0.5-x) +rnorm(length(x))/10
y2 <- 2*x^2 + 3*log(x) + rnorm(length(x))
# here is a function
modely <- function(formula, data, ...){
ll <- length(formula) # no of items in formula list
# initiation of vectors
result2 <- vector("list", ll)
result1 <- rep(NA, ll)
for(i in 1:ll) {
# some start values ***this shall be improved somehow***
n <- length(all.vars(fol[[i]])) - 2
start <- as.list(rep(mean(data[,2]),n))
names(start) <- letters[1:n]
# fit
fit<-try(nls(formula[[i]], data, start))
if( class(fit)=="try-error") result1[i] <- NA else result1[i] <-
sum(resid(fit)^2)
if( class(fit)=="try-error") result2[[i]] <- NA else result2[[i]] <-
coef(fit)
}
# ordering and combining results
ooo<-order(result1)
result <- mapply(c, "sq.resid" = result1, result2)
names(result) <- as.character(formula)
result[ooo]
}
Uwe Ligges <ligges at statistik.tu-dortmund.de> napsal dne 02.03.2009
10:54:16:
>
>
> Petr PIKAL wrote:
> > Thank you
> >
> > It was simplified version of my problem. I want to elaborate a
function
> > which can take predefined list of formulas, some data and evaluate
which
> > formulas can fit the data. I was inspired by some article in Chemical
> > engineering in which some guy used excel solver for such task. I was
> > curious if I can do it in R too. I am not sure if nls is appropriate
tool
> > for such task but I had to start somewhere.
> >
> > Here is a function which takes list of formulas and data and gives a
> > result for each formula.
> >
> > modely <- function(formula, data, ...){
> > ll <- length(formula) #no of items in formula list
> > result2 <- vector("list", ll) #prepare results
> > result1 <- rep(NA, ll)
> > for(i in 1:ll) {
> > fit<-try(nls(formula[[i]], data))
> > if( class(fit)=="try-error") result1[i] <- NA else result1[i] <-
> > sum(resid(fit)^2)
> > if( class(fit)=="try-error") result2[[i]] <- NA else result2[[i]] <-
> > coef(fit)
> > }
> >
> > ooo<-order(result1) #order results according to residual sum
> >
> > #combine results into one list together with functions used
> >
> > result <- mapply(c, "sq.resid" = result1, result2)
> > names(result) <- as.character(formula)
> > # output
> > result[ooo]
> > }
> >
> > # data
> > x <-1:10
> > y <-1/(.5-x)+rnorm(10)/100
> >
> > # list of formulas
> > fol <- structure(list(a = y ~ 1/(a - x), b = y ~ a * x^2 + b * log(x),
> > c = y ~ x^a), .Names = c("a", "b", "c"))
> >
> > modely(fol, data.frame(x=x, y=y)
> >
> > does not use "correct" model because when using default start values
it
> > results in
> >
> >> nls(fol[[1]], data.frame(x=x, y=y))
> > Error in numericDeriv(form[[3]], names(ind), env) :
> > Missing value or an infinity produced when evaluating the model
> >
> > however
> >
> > nls(fol[[1]], data.frame(x=x, y=y), start=list(a=mean(y)))
> >
> > gives correct result. Therefore I started think about how to add a
> > "better" starting value for some fits as a second part of my formula
list
> > to define structure like>
> >
> > list(a= formula1, start.formula1, b=formula2, start.formula2, ....)
> >
> > I wonder If you can push me to better direction.
>
>
> You can make up a list of lists (each containing one formula and its
> starting values) or specify formulas in one list and starting values in
> a corresponding second list.
> You need just the corresponding subsetting in your call to nls such as
> in the simple case I suggested already.
>
> Best,
> Uwe
>
>
>
> > Thanks again
> > Best regards
> > Petr
> >
> >
> >
> >
> > Uwe Ligges <ligges at statistik.tu-dortmund.de> napsal dne 02.03.2009
> > 09:41:45:
> >
> >> Petr PIKAL wrote:
> >>> Hi to all
> >>>
> >>> OK as I did not get any response and I really need some insight I
try
> >>> again with different subject line
> >>>
> >>> I have troubles with correct evaluating/structure of nls input
> >>>
> >>> Here is an example
> >>>
> >>> # data
> >>> x <-1:10
> >>> y <-1/(.5-x)+rnorm(10)/100
> >>>
> >>> # formula list
> >>> form <- structure(list(a = list(quote(y ~ 1/(a - x)),
> > "list(a=mean(y))")),
> >>> .Names = "a")
> >>>
> >>> # This gives me an error due to not suitable default starting value
> >>>
> >>> fit <- nls(form [[1]] [[1]], data.frame(x=x, y=y))
> >>>
> >>> # This works and gives me a result
> >>>
> >>> fit <- nls(form [[1]] [[1]], data.frame(x=x, y=y),
> > start=list(a=mean(y)))
> >>> *** How to organise list "form" and call to nls to enable to use
other
> >
> >>> then default starting values***.
> >>>
> >>> I thought about something like
> >>>
> >>> fit <- nls(form [[1]] [[1]], data.frame(x=x, y=y), start=get(form
> > [[1]]
> >>> [[2]]))
> >>> ^^^^^^^^^^^^^^^^^^^
> >>> but this gives me an error so it is not correct syntax. (BTW I tried
> > eval,
> >>> assign, sustitute, evalq and maybe some other options but did not
get
> > it
> >>> right.
> >>>
> >>> I know I can put starting values interactively but what if I want
them
> >
> >>> computed by some easy way which is specified by second part of a
list,
> >
> >>> like in above example.
> >> If you really want to orgnize it that way, why not simpler as in:
> >>
> >> form <- list(y ~ 1/(a - x), a = mean(y))
> >> fit <- nls(form[[1]], data.frame(x=x, y=y), start = form[2])
> >>
> >>
> >> Uwe Ligges
> >>
> >>
> >>> If it matters
> >>> WXP, R2.9.0 devel.
> >>>
> >>> Regards
> >>> Petr
> >>>
> >>> petr.pikal at precheza.cz
> >>>
> >>> ______________________________________________
> >>> R-help at r-project.org mailing list
> >>> https://stat.ethz.ch/mailman/listinfo/r-help
> >>> PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html
> >>> and provide commented, minimal, self-contained, reproducible code.
> >
More information about the R-help
mailing list