Hi
I adopted an idea from Chemical Engineering where one guy used set of
functions to check which can fit data well and used Excel for it. So I
tried if I can do it in R.
It is very rough and needs data to be a two column data frame with x in
the first column together with list of formulas.
You can use it like
nls.find(modely1, data)
and you can try to give it some starting values but it can arrive to
results without them. You can also set start=T, which will use default
method for nls. The result is list of model parameters ordered by residual
sum of residual sqares. It does not help you to think about your data, it
does not use any type of diagnostics.
But you can use it with some care.
regards
Petr
Below is three set of models and the function itself.
"modely1" <- list(
y~1/(a-x),
y~x^2-a*x,
y~a^x,
y~1-a^x,
y~a^(1/x),
y~x^a,
y~1-1/(x^a),
y~a^(1-log(x)),
y~1-exp(-a*x),
y~a*x^2-log(x),
y~a/cosh(x),
y~a*cosh(x/a))
"modely2" <- list(
y~a*x+b,
y~1/(a+b*x),
y~a/(b+x),
y~a+b*log(x),
y~a*exp(b*x),
y~a*exp(b/x),
y~a*(1-exp(b*x)),
y~(a/(a-b))*(exp(-b*x)-exp(-a*x)),
y~1/(1+exp(-a*(x-b))),
y~a*x^b,
y~a*(1+x)^b,
y~(a+x)^b,
y~(a-x)^b,
y~1-1/(1+a*x)^b,
y~a*x^(b*x),
y~a*x^(b/x),
y~x/(a+b*x-(a+b)*x^2),
y~a*cos(x-b),
y~a/(sin(x-b))
)
"modely3" <- list(
y~a+b*x+c*x^2,
y~1/(a+b*x+c*x^2),
y~a+b/x+c/(x^2),
y~a+b*log(x-c),
y~a+b*log(1+c^x),
y~a+b*x^c,
y~a*(x-b)^c,
y~a*(x/b)^c,
y~a*b^x+x^c,
y~a*x^(b*x-c),
y~(a*x^c)/(1+b*x^c),
y~(a*x^c)/(b^c+x^c),
y~a+b*exp(c*x),
y~a-b*exp(-c*x),
y~a*(1-b*exp(c*x)),
y~a*(1-exp(-b*x))^c,
y~a+(b*(exp(c*x)-1))/c,
y~exp(a+b*x+c*x^2),
y~a/(1+exp(-b*(x-c))),
y~a*sin((pi*(x-b))/c)
)
nls.find <- function(formula, data, start=NULL, ...){
ll <- length(formula)
result2 <- vector("list", ll)
result1 <- rep(NA, ll)
names(data) <- c("x","y")
for(i in 1:ll) {
if (is.null(start)) {
n <- length(all.vars(formula[[i]]))-2
start <- as.list(rep(mean(data[,2], na.rm=T),n))
names(start) <- letters[1:n]
fit<-try(nls(formula[[i]], data, start=start, ...))
} else {
if (is.logical(start)) {
fit<-try(nls(formula[[i]], data, ...))
} else {
fit<-try(nls(formula[[i]], data, start, ...)) }
}
if( class(fit)=="try-error") result1[i] <- NA else result1[i]
<-
sum(resid(fit)^2, na.rm=T)
if( class(fit)=="try-error") result2[[i]] <- NA else result2[[i]]
<-
coef(fit)
}
ooo<-order(result1)
result <- mapply(c, "sq.resid" = result1, result2, SIMPLIFY=FALSE)
names(result) <- paste(1:ll, as.character(formula), sep="/")
result[ooo]
}
Petr Pikal
petr.pikal at precheza.cz
724008364, 581252140, 581252257
r-help-bounces at r-project.org napsal dne 02.04.2009 18:40:43:
> Hello,
> Is there a function in r to find the best fitting model for a set of
data?> I would like to know if my data are related exponentially,linearly or if
there> is a logarithmic correlation between my x and y values.
> To get a better imagination I've added the graphics at the end of this
mail> as an attachment.
> Thanks, Ben
> --
>
> ______________________________________________
> 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.