如何使这个R代码(for循环)更有效

Question

我正在做一个模拟研究，我写了下面的R代码。有没有办法在不使用两个for循环的情况下编写这段代码，或者使其更高效(运行得更快)?

S = 10000
n = 100
v = c(5,10,50,100)
beta0.mle = matrix(NA,S,length(v))  #creating 4 S by n NA matrix 
beta1.mle = matrix(NA,S,length(v))
beta0.lse = matrix(NA,S,length(v))
beta1.lse = matrix(NA,S,length(v))
for (j in 1:length(v)){
  for (i in 1:S){
    set.seed(i)
    beta0 = 50
    beta1 = 10
    x = rnorm(n)
    e.t = rt(n,v[j])
    y.t = e.t + beta0 + beta1*x
    func1 = function(betas){
      beta0 = betas[1]
      beta1 = betas[2]
      sum = sum(log(1+1/v[j]*(y.t-beta0-beta1*x)^2))
      return((v[j]+1)/2*sum)
    }
    beta0.mle[i,j] = nlm(func1,c(1,1),iterlim = 1000)$estimate[1]
    beta1.mle[i,j] = nlm(func1,c(1,1),iterlim = 1000)$estimate[2]
    beta0.lse[i,j] = lm(y.t~x)$coef[1]
    beta1.lse[i,j] = lm(y.t~x)$coef[2]
  }
}

第二个for环内的func1函数用于nlm函数(在误差为t分布时查找mle)。我想在R中使用parallel包，但我没有找到任何有用的功能。

Answer 1

让任何东西在R中运行得更快的关键是用矢量化函数(例如apply家族)替换for循环。此外，对于任何编程语言，您应该查找使用相同参数多次调用昂贵函数(如nlm)的地方，并查看可以将结果存储在何处，而不是每次都重新计算。

在这里，我像您一样开始定义参数。此外，由于beta0和beta1总是50和10，我也将在这里定义它们。

S <- 10000
n <- 100
v <- c(5,10,50,100)
beta0 <- 50
beta1 <- 10

接下来，我们将在循环之外定义func1，以避免每次都重新定义它。func1现在有两个额外的参数，v和y.t，所以它可以用新的值调用。

func1 <- function(betas, v, y.t, x){
  beta0 <- betas[1]
  beta1 <- betas[2]
  sum <- sum(log(1+1/v*(y.t-beta0-beta1*x)^2))
  return((v+1)/2*sum)
}

现在我们做真正的工作了。我们没有使用嵌套循环，而是使用嵌套apply语句。外部lapply将为v的每个值制作一个列表，内部vapply将为S的每个值制作一个矩阵，用于获取四个值(beta0.mle, beta1.mle, beta0.sle, beta1.lse)。

values <- lapply(v, function(j) vapply(1:S, function(s) {
  # This should look familiar, it is taken from your code
  set.seed(s)
  x <- rnorm(n)
  e.t <- rt(n,j)
  y.t <- e.t + beta0 + beta1*x
  # Rather than running `nlm` and `lm` twice, we run it once and store the results
  nlmmod <- nlm(func1,c(1,1), j, y.t, x, iterlim = 1000)
  lmmod <- lm(y.t~x)
  # now we return the four values of interest
  c(beta0.mle = nlmmod$estimate[1],
    beta1.mle = nlmmod$estimate[2],
    beta0.lse = lmmod$coef[1],
    beta1.lse = lmmod$coef[2])
}, numeric(4)) # this tells `vapply` what to expect out of the function
)

最后我们可以把所有的东西重新组织成四个矩阵。

beta0.mle <- vapply(values, function(x) x["beta0.mle", ], numeric(S))
beta1.mle <- vapply(values, function(x) x["beta1.mle", ], numeric(S))
beta0.lse <- vapply(values, function(x) x["beta0.lse.(Intercept)", ], numeric(S))
beta1.lse <- vapply(values, function(x) x["beta1.lse.x", ], numeric(S))

最后需要注意的是，根据您使用S索引设置种子的原因，可能需要重新组织它以更快地运行。如果知道用什么种子生成x和rnorm是很重要的，那么这可能是我能做的最好的。但是，如果您这样做只是为了确保所有v的值都在x的相同值上进行测试，那么我们可以做更多的重组，这样可以使用replicate产生更多的速度。