이 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] 
    } 
} 

nlm 기능을 위해 사용되는 제 2 내부 루프 for 함수 func1 (t 오류가 분산 될 때 찾을 MLE). R에서 parallel 패키지를 사용하고 싶었지만 유용한 기능을 찾지 못했습니다.

Answer 1

R에서 더 빠르게 실행하는 데있어 핵심은 for 루프를 벡터화 된 기능 (예 : apply 제품군)으로 대체하는 것입니다. 또한 모든 프로그래밍 언어에 대해 동일한 매개 변수를 사용하여 값 비싼 함수 (예 : 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) 
} 

이제 실제 작업을 수행합니다. 중첩 된 루프가 아니라 중첩 된 적용 문을 사용합니다. 외부 lapply는 S의 각 값에 대해 (beta0.mle, beta1.mle, beta0.sle, beta1.lse) 싶어 네 개의 값에 대한 매트릭스를 만들 것입니다 v과 내부 vapply의 각 값에 대한 목록을 만들 것입니다.

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 
) 

마지막으로 모든 것을 4 개의 행렬로 재구성 할 수 있습니다.

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을 사용하여 속도를 높일 수있는 더 많은 재구성이 필요할 수 있습니다.