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.
