2013-04-16 1 views
3

두 개의 RasterStack 객체로 작업하고 있습니다. 각 객체는 단일 시간 단계를 나타내는 10 개의 레이어로 구성되어 있습니다.RasterStack 객체의 단일 셀을 통한 속도 최적화 루프

# Mock data 
pred.rst.stck <- do.call("stack", lapply(seq(10), function(i) { 
    pred.rst <- raster(nrows = 15, ncols = 15, xmn= 0, xmx = 10, ymn = 0, ymx = 10) 
    pred.rst[] <- rnorm(225, 50, 10) 
    return(pred.rst) 
}) 
resp.rst.stck <- do.call("stack", lapply(seq(10), function(i) { 
    resp.rst <- raster(nrows = 10, ncols = 10, xmn = 0, xmx = 10, ymn = 0, ymx = 10) 
    resp.rst[] <- rnorm(100, 50, 10) 
    return(resp.rst) 
}) 

pred.rst.stck 응답 변수 세트와 예측 변수의 집합 resp.rst.stck로서 기능한다. 예측기 RasterStack의 각 단일 셀에 대해, 응답 RasterStack의 각 셀에 선형 모델을 맞추고, 맞는 모델 당 해당 R- 제곱을 추출하여이를 합산하고 싶습니다. 전체 과정을 가속화하는 우아한 방법이 있는지

# Parallelization 
library(parallel) 
n.cores <- detectCores() 
clstr <- makePSOCKcluster(n.cores) 

# Extract cell values from RasterStack objects 
pred.vals <- getValues(pred) 
resp.vals <- getValues(resp) 
clusterExport(clstr, c("pred.vals", "resp.vals")) 

# Loop through all predictor cells 
rsq.sums <- parLapply(clstr, seq(nrow(pred.vals)), function(i) { 
    # For each predictor cell, loop through all response cells, 
    # fit linear model, extract and sum up single R-squared 
    do.call("sum", lapply(seq(nrow(resp.vals)), function(j) { 
    summary(lm(resp.vals[j, ] ~ pred.vals[i, ]))$r.squared 
    })) 
}) 

parLapply 비록 수행 훨씬 더 일반 lapply에 비해, 내가 궁금 : 짧게 이야기하면, 여기에 지금까지 R의 parallel 패키지를 사용하여 내 가장 빠른 방법입니다. 어떤 제안?

건배,
플로리안

답변

5

당신이 시도 할 수있는 몇 가지 요령이 있습니다. 선형 모델을 만드는 방법을 이해할 수는 없지만 선형 모델에서 계산하는 r.squared은 Pearson의 상관 계수 (기본 인수가있는 cor)의 제곱과 같습니다. 선형 모델.

# Finding r.squared using a lm 
f1 <- function(n){summary(lm(resp.vals[n,] ~pred.vals[n,]))$r.squared} 
# Finding r.squared using Pearson's 
f2 <- function(n){ cor(resp.vals[n,],pred.vals[n,])^2} 

# Both give the same result 
f1(1) 
[1] 0.0009032986 
f2(1) 
[1] 0.0009032986 

그리고 이러한 기능의 단일 통화의 타이밍 측면에서

:

는 데이터를 사용하여 이러한 두 가지 기능을 비교

require(microbenchmark) 
microbenchmark(f1(1) , f2(1)) 
Unit: microseconds 
    expr  min  lq median  uq  max neval 
f1(1) 2009.328 2037.0680 2071.045 2124.9225 4102.079 100 
f2(2) 73.075 77.7365 84.690 89.7215 5485.368 100 

그래서 당신은 당신의 처리 시간을 절감 할 수 있어야 25- lm 대신 cor을 사용하여

cov()^2의 사용을 위해 원래의 기능을 교환하는 빠른 시스템 시간 비교는이 경우로 보여줍니다

#Using cov() 
    user system elapsed 
    0.013 0.002 1.085 

#Using lm() 
    user system elapsed 
    0.159 0.028 26.334 
+0

믿을 수없는, 대단히 감사합니다! – fdetsch

+0

@flowla 당신을 환영합니다! :-) –