我正在使用两个 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 平方并将其相加。长话短说,这是迄今为止我使用 R 的 parallel
包最快的方法:
# 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.squared
相当于 PIL 逊相关系数的平方 (cor
在 R 中使用默认参数),其计算速度比线性模型快得多。
使用您的数据比较这两个函数:
# 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
因此,通过使用 cor
而不是 lm
,您应该能够将处理时间缩短 25 倍。
将原始函数替换为使用 cov()^2
的快速系统时间比较表明情况确实如此:
#Using cov()
user system elapsed
0.013 0.002 1.085
#Using lm()
user system elapsed
0.159 0.028 26.334
关于r - 通过 RasterStack 对象的单个单元格进行速度优化的循环,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/16042154/