问题描述
我有一个数据集，其中包含由多个模型预测的基于年份的数据，采用 data.table 格式。
library(data.table)
nYears = 20 # real data: 110
nMod = 3 # real data: ~ 100
nGrp = 45
dataset <- data.table(
group_code = rep(seq(1:nGrp ),times= 3*nYears ),Year = rep(seq(1:nYears ),each=nGrp ),value = rnorm(2700,mean = 10,sd = 2),var1 = rep (rnorm(nGrp,mean = nMod,sd = 1),times= nMod*nYears ),var2 = rep (rnorm(nGrp,mean = 1.5,sd = 0.5),model = as.character(rep(seq( from = 1,to = nMod ),each=nGrp *nYears ))
)
setkey(dataset,Year,model)
我需要根据这个 数据集 执行一组计算，这个向量名为 x，长度为 1001，包含在 seq(-2,8,by=0.01) 上。
为此，我创建了一个新的 data.table (dt)，其中包含 dataset 的重复版本以合并向量 x，相应地：
dt <- dataset[,lapply(.SD,function(x) rep(x,1001))]
dt[,x := rep(round(seq(-2,by=0.01),2),each= nYears*nGrP*nMod) ]
由于我的原始数据集包含数百个模型，因此此操作内存效率不高。
我需要的最重要的操作包括生成 x 的正态分布，均值 = var1 和 sd = var2，通过 group_code、Year 和模型。例如：
# key computation
dt [,norm_dist := dnorm (x,var1,var2),by= .(group_code,model )]
最后一个操作在我的桌面上非常快。但是，我还有其他操作要执行，需要对 data.table 进行子集化并且占用大量 RAM。一个例子：
dt[ x %between% c( 2,5.99),dt2 := rep_len( rev(dt [x %between% c(-2,1.99)]$value),length.out=.N),by= .(Year,model) ]
弹出以下错误：
Error: cannot allocate vector of size 1.3 Gb
我认为这个具体步骤中的问题与子集和 rev() 函数有关。
尽管如此，我用来根据 data.table dt 中的向量“x”执行一组计算的方法似乎并不合适，因为我将数据集与我需要计算的向量（“x”）。
我希望有人能教我如何有效地改进我的代码，因为我在原始数据集中有相当多的模型，大大增加了它的大小。
谢谢！
解决方法
我觉得这部分代码应该更清楚
dt[ x %between% c( 2,5.99),dt2 := rep_len( rev(dt [x %between% c(-2,1.99)]$value),length.out=.N),by= .(Year,model) ]
因为它对我来说有点像一个黑匣子。特别是因为这种双重子集是您的问题产生的地方。
这些代码位 x %between% c( 2,5.99) 和 dt[x %between% c(-2,1.99)] 应该在您的所有情况下始终产生相同的位置。您应该在代码中考虑这一点，以提高效率。
试试这样的事情让事情更清楚一点：
by_YM <- split(dt,by=c("Year","model"))
ind1 <- which(by_YM[[1]][["x"]] %between% c( 2,5.99))
ind2 <- which(by_YM[[1]][["x"]] %between% c(-2,1.99))
for(i in 1:length(by_YM)){
dt_i <- by_YM[[i]]
#val1 <- rep_len(rev(dt_i$value[ind2]),length.out=length(ind1)) #val1 is equal to val,no need for rep_len
val <- rev(dt_i$value[ind2])
by_YM[[i]] <- dt_i[ind1,dt2 := val]
}
但是我们的 dt2 列不相等，但由于我不确定最终结果应该如何，因此我无法调试它。
dt2_a <- dt[Year == 20 & model == 3,dt2]
dt2_b <- by_YM[["20.3"]][,dt2]
test <- cbind(dt2_a,dt2_b)
第二个代码也快得多。
library(microbenchmark)
microbenchmark( "new_code" = {
by_YM <- split(dt,"model"))
ind1 <- which(by_YM[[1]][["x"]] %between% c( 2,1.99))
for(i in 1:length(by_YM)){
dt_i <- by_YM[[i]]
val1 <- rep_len(rev(dt_i$value[ind2]),dt2 := val]
}},"old_code" = dt[ x %between% c( 2,model) ],times = 5)
Unit: milliseconds
expr min lq mean median uq max neval cld
new_code 155.426 156.4916 200.6587 185.0347 188.9436 317.3977 5 a
old_code 1290.909 1299.8570 1398.6866 1370.4526 1471.0569 1561.1574 5 b
试试看，祝你好运