用Loop在R中求解多方程组

vc9ivgsu  于 2023-03-15  发布在  其他
关注(0)|答案(1)|浏览(198)

我有一个非线性方程组的列表,并希望解决他们同时使用循环。我的目标是找到每个方程的参数和r平方。

x<-c(1,2,3,4,5,6,7)
y<-c(0.12,0.15,0.19,0.29,0.52,0.60,0.66)

下面是方程的列表。

y~a*x+b
y~a*(x-b)
y~a+b*x+c*x^0.5
y~a*x/(b+x)

我可以单独求解每个方程并找到它的参数,但我有20多个方程,不知道是否有一种方法可以通过循环或其他方式求解所有方程。

fit1<-nls(y~a*x+b,start=list(a=1,b=1))
fit2<-nls(y~a*(x-b),start=list(a=1,b=1))
fit3<-nls(y~a+b*x+c*x^0.5,start=list(a=1,b=1,c=1)) 
fit4<-nls(y~a*x/(b+x),start=list(a=1,b=1))

r方形

Rsq_fit1<-cor(y,predict(fit1))^2
Rsq_fit2<-cor(y,predict(fit2))^2
Rsq_fit3<-cor(y,predict(fit3))^2
Rsq_fit4<-cor(y,predict(fit4))^2

系数和R平方

list(coef(fit1),coef(fit2),coef(fit3),coef(fit4))
list(Rsq_fit1,Rsq_fit2,Rsq_fit3,Rsq_fit4)
nwlls2ji

nwlls2ji1#

如果将xy放入 Dataframe 中,并将方程放入列表中:

df <- data.frame(x, y)

eq <- list(y~a*x+b,
           y~a*(x-b),
           y~a+b*x+c*x^0.5,
           y~a*x/(1+x))

然后你可以这样做:

do.call(rbind, lapply(eq, function(x) {
  v <- setdiff(all.vars(x), c("x", "y"))
  mod <- nls(x, df, start = setNames(as.list(rep(1, length(v))), v))
  rsq <- cor(df$y, predict(mod))^2
  EQ <- do.call("substitute", list(x, env = as.list(round(coef(mod), 3))))
  data.frame(equation = deparse(EQ), r_squared = rsq)
}))
#>                                 equation r_squared
#> 1                 y ~ 0.102 * x + -0.046 0.9397561
#> 2                y ~ 0.102 * (x - 0.449) 0.9397561
#> 3 y ~ 0.337 + 0.222 * x + -0.449 * x^0.5 0.9628876
#> 4                  y ~ 0.503 * x/(1 + x) 0.6932539

请注意,您的第四个方程不收敛,因此为了便于演示,我对它做了一些修改。
如果您希望以数值格式输出结果,并且即使nls失败也能输出结果,您可以执行以下操作:

p <- unique(unlist(sapply(eq, function(x) setdiff(all.vars(x), c("x", "y")))))

do.call(rbind, lapply(eq, function(x) {
  v <- setdiff(all.vars(x), c("x", "y"))
  coefdf <- do.call(data.frame, setNames(as.list(rep(NA, length(p))), p))
  mod <- tryCatch({
    nls(x, df, start = setNames(as.list(rep(1, length(v))), v))
  }, error = function(e) "error")
  if(class(mod) == "character") {
    return(cbind(
       data.frame(equation = deparse(x)), 
       converged = FALSE,
       coefdf,
       r_squared = NA))
  }
  coefdf[1, v] <- round(coef(mod), 3)
  cbind(equation = deparse(x), 
        converged = TRUE,
        coefdf,
        r_squared = cor(df$y, predict(mod))^2)
}))
#>                    equation converged     a      b      c r_squared
#> 1             y ~ a * x + b      TRUE 0.102 -0.046     NA 0.9397561
#> 2           y ~ a * (x - b)      TRUE 0.102  0.449     NA 0.9397561
#> 3 y ~ a + b * x + c * x^0.5      TRUE 0.337  0.222 -0.449 0.9628876
#> 4         y ~ a * x/(b + x)     FALSE    NA     NA     NA        NA

相关问题