R中的嵌套for循环：查找拟合参数答案

【问题标题】：Nested for loops in R: find fit parametersR中的嵌套for循环：查找拟合参数
【发布时间】：2017-04-25 08:39:48
【问题描述】：

我有实验数据点，我想与我的模型相匹配，该模型有 3 个参数（p、q 和 R）。因此，我有三个嵌套的 for 循环来遍历三个参数的所有可能组合，然后确定最佳拟合（残差的最小二乘）。我就是这样做的（权重 = 1/exp(y)）：

tgreparpq <- function(x, p, q, R){exp(-(p + q)*x)*(1 + x*(q + R*p) + x^2*((R*q^2)/2 + R*q*p) + x^3*R^2*q^2*p/2)}

pvalues <- seq(1e-10,1.01,0.01)
qvalues <- seq(1e-10,1.01,0.01)
repair <- seq(1e-10,1.01,0.01)

bestl <- array(dim = c(length(pvalues), length(qvalues)))
sigmapq <- array(dim = c(length(pvalues), length(qvalues)))
bestp <- vector(length = klength)
bestq <- vector(length = klength)
bestr <- array(dim = c(length(pvalues), length(qvalues)))
bestrk <- vector(length = klength)

#number of dose rates = klength
klength = 2

for (k in 1:klength){
  x <- c(0, dat$Dose[dat$Rate%in%rle(dat$Rate)$values[k]])
  SR <- vector(length = length(repair))
  for (i in 1:length(pvalues)){
    for (j in 1:length(qvalues)){
      for (l in 1:length(repair)){

    y <- tgreparpq(x,pvalues[i], qvalues[j], repair[l])
    SR[l] <- sum(((y-c(1,dat$Survival[dat$Rate%in%rle(dat$Rate)$values[k]]))*1/exp(y))^2)

      }
     bestl[i,j] <- which(SR %in% min(SR))
     sigmapq[i,j] <- SR[bestl[i,j]]
     bestr[i,j] <- repair[bestl[i,j]]
    }
  }
  besti <- which(sigmapq == min(sigmapq), arr.ind = TRUE)[1]
  bestp[k] <- pvalues[besti]
  bestj <- which(sigmapq == min(sigmapq), arr.ind = TRUE)[2]
  bestq[k] <- qvalues[bestj]
  bestrk[k] <- bestr[besti,bestj]
}

这是一个相当缓慢的过程，我知道您不应该在 r 中使用那么多 for 循环。因此我的问题是：有没有更好的方法来确定合适的参数（即有没有办法替换 for 循环）？

编辑： 以下是一些示例数据：

dat =

Dose        Survival     Rate

1.9163E+00  6.42870E-01  3.0000E+01    
3.9713E+00  3.68150E-01  3.0000E+01    
5.9857E+00  1.76050E-01  3.0000E+01    
7.9572E+00  8.27670E-02  3.0000E+01    
1.0013E+01  2.01370E-02  3.0000E+01    
1.2015E+01  1.09200E-02  3.0000E+01
1.9683E+00  6.42530E-01  7.6800E+01    
2.9740E+00  4.86220E-01  7.6800E+01    
4.0354E+00  3.09730E-01  7.6800E+01    
5.0276E+00  2.13930E-01  7.6800E+01    
6.0851E+00  1.67200E-01  7.6800E+01    
7.0223E+00  1.04640E-01  7.6800E+01    
8.0531E+00  6.79020E-02  7.6800E+01    
9.0841E+00  3.14080E-02  7.6800E+01    
1.0135E+01  2.12510E-02  7.6800E+01    
1.1176E+01  8.39810E-03  7.6800E+01    
1.2168E+01  4.92070E-03  7.6800E+01    
1.4169E+01  1.53690E-03  7.6800E+01

【问题讨论】：

请提供reproducible example。当没有 dat 提供给测试时，显示你的整个代码是没有意义的。
tgreparpq的功能是什么？
我添加了它。这是一个将生存描述为剂量函数的函数（使用三个参数）。
所以您实际上是在对代码进行非线性拟合并扫描所有可能的参数值以最小化 MSE？我得到正确了吗？你可以检查这个stat.ethz.ch/R-manual/R-devel/library/stats/html/nls.html。也许它可以帮助你......
什么是klength？

标签： r performance for-loop

【解决方案1】：

不是答案，但我冒昧地准备了一个数据框供其他人使用您的数据。它可能会加快回答速度。

df = data.frame(Dose = c(   
  1.9163E+00, 
  3.9713E+00, 
  5.9857E+00, 
  7.9572E+00, 
  1.0013E+01, 
  1.2015E+01, 
  1.9683E+00, 
  2.9740E+00, 
  4.0354E+00, 
  5.0276E+00, 
  6.0851E+00, 
  7.0223E+00, 
  8.0531E+00, 
  9.0841E+00, 
  1.0135E+01, 
  1.1176E+01, 
  1.2168E+01, 
  1.4169E+01),
  Survival = c(
  6.42870E-01,
  3.68150E-01,
  1.76050E-01,
  8.27670E-02,
  2.01370E-02,
  1.09200E-02,
  6.42530E-01,
  4.86220E-01,
  3.09730E-01,
  2.13930E-01,
  1.67200E-01,
  1.04640E-01,
  6.79020E-02,
  3.14080E-02,
  2.12510E-02,
  8.39810E-03,
  4.92070E-03,
  1.53690E-03),
  Rate = c( 
  3.0000E+01,
  3.0000E+01,
  3.0000E+01,
  3.0000E+01,
  3.0000E+01,
  3.0000E+01,
  7.6800E+01,
  7.6800E+01,
  7.6800E+01,
  7.6800E+01,
  7.6800E+01,
  7.6800E+01,
  7.6800E+01,
  7.6800E+01,
  7.6800E+01,
  7.6800E+01,
  7.6800E+01,
  7.6800E+01))

希望它可以帮助某人。我正在尝试

nls(Survival ~ tgreparpq(Dose, p, q, R), data = df, start = list(p = 0.1, q = 0.1, R = 0.1))

但我得到了一个“奇异渐变”。我正在检查，但可能您的模型没有唯一的解决方案？

【讨论】：

谢谢。 nls 的问题在于，根据我给出的起始值，我得到的结果非常不同。这就是我想“手工”做的原因。
即拟合函数复杂或参数过多时的非线性拟合问题。相信我，我正在做一个项目，我们试图用超过 10 个参数拟合一个函数，这很痛苦。如果函数有很多最小值，那么你得到的结果将很大程度上取决于初始参数。您会发现最小值更接近您的初始参数。所以知道什么是好的开始是个好主意。

【解决方案2】：

您可以使用expand.grid 为要循环的变量的每个可能组合获取一个数据框，然后使用apply，这通常要快得多。

此代码不是您上述代码的完整解决方案，但它是该想法的一个最小示例...

dat <- read.table(header = T, text = "
Dose        Survival     Rate
1.9163E+00  6.42870E-01  3.0000E+01    
3.9713E+00  3.68150E-01  3.0000E+01    
5.9857E+00  1.76050E-01  3.0000E+01    
7.9572E+00  8.27670E-02  3.0000E+01    
1.0013E+01  2.01370E-02  3.0000E+01    
1.2015E+01  1.09200E-02  3.0000E+01
1.9683E+00  6.42530E-01  7.6800E+01    
2.9740E+00  4.86220E-01  7.6800E+01    
4.0354E+00  3.09730E-01  7.6800E+01    
5.0276E+00  2.13930E-01  7.6800E+01    
6.0851E+00  1.67200E-01  7.6800E+01    
7.0223E+00  1.04640E-01  7.6800E+01    
8.0531E+00  6.79020E-02  7.6800E+01    
9.0841E+00  3.14080E-02  7.6800E+01    
1.0135E+01  2.12510E-02  7.6800E+01    
1.1176E+01  8.39810E-03  7.6800E+01    
1.2168E+01  4.92070E-03  7.6800E+01    
1.4169E+01  1.53690E-03  7.6800E+01
")

tgreparpq <- function(x, p, q, R){exp(-(p + q)*x)*(1 + x*(q + R*p) + x^2*((R*q^2)/2 + R*q*p) + x^3*R^2*q^2*p/2)}

pvalues <- seq(1e-10,1.01,0.01)
qvalues <- seq(1e-10,1.01,0.01)
repair <- seq(1e-10,1.01,0.01)

x <- 1

results <- expand.grid(pvalues, qvalues, repair)
results$y <- apply(results, 1, function(result_n) {
  tgreparpq(x, result_n[1], result_n[2], result_n[3])
})

【讨论】：