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.pi...@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 <lig...@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 <lig...@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.pi...@precheza.cz > >>> > >>> ______________________________________________ > >>> R-help@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. > > ______________________________________________ R-help@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.