如何使此 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函数（在错误分布时查找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 将为您想要获得的四个值创建一个矩阵 (beta0.mle, beta1.mle, beta0.sle, beta1.lse) 对于 S.

的每个值

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.

如何使此 R 代码（for 循环）更高效？

How to make this R code (for loop) more efficient?

r

statistics

simulation

mle