在R中使用rbinom引导需要很长时间才能运行

Question

我一直在R中使用rbinom for循环运行bootstrap，但它们运行时间太长。

我想在一个有1,500,000行的数据集上执行bootstrap。

我想重新采样这些行对于每一个重新采样的行

1. rbinom two probability ('prob1' &'prob2')变成0和1 ('prob1_ber' &prob2_ber)
2. 添加新列'paired'与步骤1的合并结果
3. 将列"配对"one_answers"正"的唯一组合组合为0和1 ('prob_final')
4. 计算'pair_FPR'和'pair_TPR'

我的代码如下:

library(boot)
#making example data
set.seed(1)
d2 <- data.frame(prob1=runif(n=1500000, min=1e-50, max=.9999999999),
prob2=runif(n=1500000, min=1e-44, max=.9999999989),
Positive=sample(c(0,1), replace=TRUE, size=1500000))
#making bootstrap function
function_1 <- function(data, i){
d2<-data[i,]

d2$prob1_ber <- rbinom(nrow(d2), 1, d2$prob1) #bernoulli 1 or 0
d2$prob2_ber <- rbinom(nrow(d2), 1, d2$prob2) #bernoulli 1 or 0

d2$paired <- ifelse(d2$prob1_ber == 1 & d2$prob2_ber == 1, '11',
ifelse(d2$prob1_ber == 0 & d2$prob2_ber ==0, '00',
ifelse(d2$prob1_ber == 1 & d2$prob2_ber ==0, '10',
ifelse(d2$prob1_ber == 0 & d2$prob2_ber ==1, '01', NA)))) 

d2$prob_final <- ifelse(d2$paired == '00',d2$prob1_ber, NA) #if both negative then negative

for (i in which(d2$paired =='11' & d2$Positive==1)) {
d2$prob_final[i] <- rbinom(1,1,0.9)
}
for (i in which(d2$paired =='11' & d2$Positive==0)) {
d2$prob_final[i] <- rbinom(1,1,0.5)
}
for (i in which(d2$paired =='01' & d2$Positive==1)) {
d2$prob_final[i] <- rbinom(1,1,0.8)
}
for (i in which(d2$paired =='01' & d2$Positive==0)) {
d2$prob_final[i] <- rbinom(1,1,0.1)
}
for (i in which(d2$paired =='10' & d2$Positive==1)) {
d2$prob_final[i] <- rbinom(1,1,0.7)
}
for (i in which(d2$paired =='10' & d2$Positive==0)) {
d2$prob_final[i] <- rbinom(1,1,0.2)
}

pair_FPR <- sum(d2[which(d2$Positive==0),]$prob_final) / nrow(d2[which(d2$Positive==0),])*100

pair_TPR <- sum(d2[which(d2$Positive==1),]$prob_final) / nrow(d2[which(d2$Positive==1),])*100

return(c(pair_FPR, pair_TPR))
}

set.seed(1)
boot_out <- boot(d2, function_1, 1000)
print(boot_out)

这个引导程序运行时间太长(n=1000)。有没有办法让它快一点?

多谢!

Answer 1

有一个很好的理由说明"如果您正在使用R并考虑使用for循环，那么可能有更好的方法"。我认为这是一个很好的例子。

你没有给出你的总体目标的上下文或描述，我也没有花时间理解你的代码。我也很困惑，为什么你在某些地方利用了R的矢量化，而在其他地方却没有。

而且，我认为使用boot库是在转移注意力。重要的是函数function_1的底层性能。最后，我认为没有必要生成1.5亿甚至150万的观测值来调查底层性能。

因此，我试图改进你的功能是:

function_2 <- function(data, i){
d2<-data[i,] %>% 
mutate(
prob1_ber=rbinom(nrow(.), 1, prob1), #bernoulli 1 or 0
prob2_ber=rbinom(nrow(.), 1, prob2), #bernoulli 1 or 0
paired=ifelse(prob1_ber == 1 & prob2_ber == 1, '11',
ifelse(prob1_ber == 0 & prob2_ber ==0, '00',
ifelse(prob1_ber == 1 & prob2_ber ==0, '10',
ifelse(prob1_ber == 0 & prob2_ber ==1, '01', NA)))), 
dprob_final=case_when(
paired == '00' ~ prob1_ber,
paired =='11' & Positive==1 ~ rbinom(1,1,0.9),
paired =='11' & Positive==0 ~ rbinom(1,1,0.5),
paired =='01' & Positive==1 ~ rbinom(1,1,0.8),
paired =='01' & Positive==0 ~ rbinom(1,1,0.1),
paired =='10' & Positive==1 ~ rbinom(1,1,0.7),
paired =='10' & Positive==0 ~ rbinom(1,1,0.2)
)
)

pair_FPR <- sum(d2[which(d2$Positive==0),]$prob_final) / nrow(d2[which(d2$Positive==0),])*100

pair_TPR <- sum(d2[which(d2$Positive==1),]$prob_final) / nrow(d2[which(d2$Positive==1),])*100

return(c(pair_FPR, pair_TPR))
}

我的测试数据是

N <- 15000
#making example data
set.seed(1)
d2 <- data.frame(prob1=runif(n=N, min=1e-50, max=.9999999999),
prob2=runif(n=N, min=1e-44, max=.9999999989),
Positive=sample(c(0,1), replace=TRUE, size=1500000))

注意function_1(d2, i)的结果将与function_2(d2, i)的结果不一样，因为生成随机数的顺序。(function_2按顺序从第1行工作到第n行，function_1 works through rows in groups defined by与and配对。)但是，我相信这两个函数的分布性质是相同的。

那么，为了比较性能…

library(microbenchmark)
microbenchmark(
list=list(
"f1"= function_1(d2, 1), 
"f2"= function_2(d2, 1)
), 
times=10
)
Unit: nanoseconds
expr min lq mean median uq max neval
f1   7  9 28.7      9  9 203    10
f2   8  9  8.9      9  9  10    10

执行时间平均相对减少100 *(27.7 - 8.9)/27.7 = 67.8%。相对性能很可能取决于N，但我预计N的好处会增加，因为向量化相对于循环的好处应该随着N的增加而增加。

请记住，使用tidyverse，虽然代码通常易于阅读和维护，但通常不会提供最快的执行时间。data.table和base R通常更优。

我把改进我的努力留给别人。我相信可以做到。