以下是具有 158K 个名为“sh_data”的观测值的大型数据框的子集。
Patient_ID Age_in_days DEMAdmNo
396076 28542 0
396076 28570 0
396076 28598 0
396076 28626 0
396076 28654 0
396076 28682 0
396076 28710 0
396076 28738 0
396076 28766 0
396076 28794 0
396076 28822 0
396076 28850 0
396076 28878 0
396076 28906 0
396076 28934 0
396076 28962 0
396076 28990 0
396076 29018 0
396076 29046 0
396076 29074 0
396076 29102 1
396076 29165 0
396076 29200 0
396076 29228 0
396076 29263 0
396076 29200 0
396076 29228 0
396076 29263 0
我正在尝试计算过去六个月中第 3 列为 1(表示为 LACE_E)的记录的实例数。因此，对于第一个记录，年龄最小的地方将为零。对于第二条记录，如果天数差异 <= 183 天并且第一条记录的第 3 列为零，则它将为 1，依此类推。
我在 R 中编写了以下查询:
LACE_E <- numeric(0)
for(i in 1:length(sh_data[,1]))
{
LACE_E[i] = 0
for(j in 1:length(sh_data[,1]))
{
if(sh_data$Patient_ID[i] == sh_data$Patient_ID[j] & sh_data$Age_in_days[i] > sh_data$Age_in_days[j] & (sh_data$Age_in_days[i]- sh_data$Age_in_days[j])<= 183 & sh_data$DEMAdmNo[j] == 1)
{LACE_E[i] = LACE_E[i] + 1}
}
}
此查询需要很长时间才能处理。在我的系统中处理 100 行需要 1 小时。请帮忙!!
最佳答案
您无需转到 Rcpp 或 data.table 即可获得非常显着的改进。
获取原始数据并复制它以获得更多可用的时间:
d <- read.table(head = TRUE, text =
"Patient_ID Age_in_days DEMAdmNo
396076 28542 0
396076 28570 0
396076 28598 0
396076 28626 0
396076 28654 0
396076 28682 0
396076 28710 0
396076 28738 0
396076 28766 0
396076 28794 0
396076 28822 0
396076 28850 0
396076 28878 0
396076 28906 0
396076 28934 0
396076 28962 0
396076 28990 0
396076 29018 0
396076 29046 0
396076 29074 0
396076 29102 1
396076 29165 0
396076 29200 0
396076 29228 0
396076 29263 0
396076 29200 0
396076 29228 0
396076 29263 0 ")
d <- rbind(d, d, d, d, d, d, d, d, d, d)
作为函数和计时运行的原始代码:
f0 <- function(sh_data) {
LACE_E <- numeric(0)
for(i in 1:length(sh_data[,1])) {
LACE_E[i] = 0
for(j in 1:length(sh_data[,1])) {
if(sh_data$Patient_ID[i] == sh_data$Patient_ID[j] &
sh_data$Age_in_days[i] > sh_data$Age_in_days[j] &
(sh_data$Age_in_days[i]- sh_data$Age_in_days[j])<= 183 &
sh_data$DEMAdmNo[j] == 1) {
LACE_E[i] = LACE_E[i] + 1
}
}
}
}
system.time(v0 <- f0(d))
## user system elapsed
## 4.803 0.007 4.812
分析显示大约 90% 的时间花在了在内循环中使用 $ 提取列上:
Rprof()
v0 <- f0(d)
Rprof(NULL)
head(summaryRprof()$by.total)
## "f0" 4.94 100.00 0.60 12.15
## "$" 4.24 85.83 0.72 14.57
## "$.data.frame" 3.52 71.26 0.36 7.29
## "[[" 3.16 63.97 0.46 9.31
## "[[.data.frame" 2.70 54.66 0.96 19.43
## "%in%" 0.92 18.62 0.22 4.45
将列提取移出循环可显着提高性能:
f1 <- function(sh_data) {
LACE_E <- numeric(0)
Patient_ID <- sh_data$Patient_ID
Age_in_days <- sh_data$Age_in_days
DEMAdmNo <- sh_data$DEMAdmNo
for(i in 1:length(sh_data[,1])) {
LACE_E[i] = 0
for(j in 1:length(sh_data[,1])) {
if(Patient_ID[i] == Patient_ID[j] &
Age_in_days[i] > Age_in_days[j] &
(Age_in_days[i]- Age_in_days[j])<= 183 &
DEMAdmNo[j] == 1) {
LACE_E[i] = LACE_E[i] + 1
}
}
}
}
system.time(v1 <- f1(d))
## user system elapsed
## 0.163 0.000 0.164
从空结果开始并增长它几乎总是一个坏主意；预先分配结果是更好的做法。在这种情况下，算法已经是 O(n^2)，因此您不会注意到太多，但它确实有所不同，尤其是在添加其他改进之后。 f2 预分配结果:
f2 <- function(sh_data) {
n <- nrow(sh_data)
LACE_E <- numeric(n)
Patient_ID <- sh_data$Patient_ID
Age_in_days <- sh_data$Age_in_days
DEMAdmNo <- sh_data$DEMAdmNo
for(i in 1:n) {
LACE_E[i] = 0
for(j in 1:n) {
if(Patient_ID[i] == Patient_ID[j] &
Age_in_days[i] > Age_in_days[j] &
(Age_in_days[i]- Age_in_days[j])<= 183 &
DEMAdmNo[j] == 1) {
LACE_E[i] = LACE_E[i] + 1
}
}
}
}
system.time(v2 <- f2(d))
## user system elapsed
## 0.147 0.000 0.148
使用正确的逻辑运算符 && 而不是 & 可以进一步改进:
f3 <- function(sh_data) {
n <- nrow(sh_data)
LACE_E <- numeric(n)
Patient_ID <- sh_data$Patient_ID
Age_in_days <- sh_data$Age_in_days
DEMAdmNo <- sh_data$DEMAdmNo
for(i in 1:n) {
LACE_E[i] = 0
for(j in 1:n) {
if(Patient_ID[i] == Patient_ID[j] &&
Age_in_days[i] > Age_in_days[j] &&
(Age_in_days[i] - Age_in_days[j]) <= 183 &&
DEMAdmNo[j] == 1) {
LACE_E[i] = LACE_E[i] + 1
}
}
}
}
system.time(v3 <- f3(d))
## user system elapsed
## 0.108 0.002 0.111
这些都是您需要转到Rcpp 的步骤，但您不需要转到Rcpp 来执行它们。
为了获得更快的速度，您可以进行字节编译:
f3c <- compiler::cmpfun(f3)
system.time(v3 <- f3c(d))
## user system elapsed
## 0.036 0.000 0.036
这些计算是在 R 3.1.3 中完成的。 微基准测试总结:
microbenchmark(f0(d), f1(d), f2(d), f3(d), f3c(d), times = 10)
## Unit: milliseconds
## expr min lq mean median uq max neval cld
## f0(d) 5909.39756 5924.8493 5963.63608 5947.23469 6011.94567 6048.03571 10 d
## f1(d) 196.16466 197.3252 200.22471 197.93345 202.49236 210.22011 10 c
## f2(d) 187.68169 190.5644 194.02454 192.47596 195.63821 204.27415 10 c
## f3(d) 109.17816 110.6695 112.55218 111.93915 114.43341 116.92342 10 b
## f3c(d) 37.37348 38.8757 39.34564 39.58563 40.50597 40.58568 10 a
R.version$version.string
## [1] "R version 3.1.3 Patched (2015-03-16 r68072)"
4月发布的R 3.2.0对解释器和字节码引擎进行了多项改进，进一步提升了性能:
## Unit: milliseconds
## expr min lq mean median uq max neval cld
## f0(d) 4351.33908 4429.71559 4471.32960 4479.13901 4499.39769 4601.05390 10 d
## f1(d) 183.57765 184.68961 190.10391 187.30951 199.56235 200.57238 10 c
## f2(d) 177.47063 181.09790 189.78291 185.58951 190.34782 233.90264 10 c
## f3(d) 105.79767 108.02553 114.48950 110.17056 112.85710 149.42474 10 b
## f3c(d) 14.41182 14.43227 14.70098 14.49289 14.84504 15.67709 10 a
R.version$version.string
## [1] "R Under development (unstable) (2015-03-24 r68072)"
良好的 R 编程实践和性能分析工具的使用可以让您走得更远。如果您想进一步改进，可以转到 Rcpp 但这可能足以满足您的目的。
关于r - 提高 R 中嵌套 For 循环的性能，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/29295504/