如何使此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來提高速度。