如何使此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回路用于nlm功能内部的功能func1（查找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来提高速度。