使用 data.table 有效地模拟数据答案

【问题标题】：Simulating Data Efficiently with data.table使用 data.table 有效地模拟数据
【发布时间】：2017-06-28 02:23:17
【问题描述】：

我正在尝试从两个较小的数据集模拟一个新数据集。对我来说保持边缘很重要从最终数据集中这些较小的数据集中计数。希望这个可重现的例子能解释什么我是说。

构建虚假数据

library(data.table) # 1.10.5
set.seed(123)
meanVal <- 40

demoDat

这里我模拟了一些年龄和性别数据。每个位置总是有 2 个性别级别和 100 个年龄级别。

demoDat <- CJ(with(CJ(letters,letters[1:5]), paste0(V1,V2)), c("M","F"), 0:99)
setnames(demoDat, c("Location","Gender","Age"))
demoDat[, Val := rpois(.N, meanVal)]


       Location Gender Age Val
    1:       aa      F   0  36
    2:       aa      F   1  47
    3:       aa      F   2  29
   ---                        
25998:       ze      M  97  45
25999:       ze      M  98  38
26000:       ze      M  99  39

时间日期

此代码模拟时间数据维度。在这种情况下，日期按周间隔，但实际数据不必遵守这种统一性。数周可能会丢失。

timeDat <- with(demoDat, CJ(unique(Location), seq(from=as.Date("2016-01-01"),by=7,length.out = 52)))
setnames(timeDat, c("Location","Date"))
totals <- demoDat[, .(Val=sum(Val)), by=.(Location)]
timeDat[totals, Val := rmultinom(1:.N, i.Val, prob=rep(1,.N)), by=.EACHI,on=.(Location)]

      Location       Date Val
   1:       aa 2016-01-01 176
   2:       aa 2016-01-08 143
   3:       aa 2016-01-15 143
  ---                        
6758:       ze 2016-12-09 165
6759:       ze 2016-12-16 142
6760:       ze 2016-12-23 156

快速对账

每个位置都应该有一个 Val 列，在 demoDat 和 timeDat 数据集中总计相同。

timeDat[, sum(Val), by=.(Location)][order(-V1)][1:5]
#    Location   V1
# 1:       jb 8229
# 2:       xb 8223
# 3:       ad 8179
# 4:       nc 8176
# 5:       gd 8173
demoDat[, sum(Val), by=.(Location)][order(-V1)][1:5]
#    Location   V1
# 1:       jb 8229
# 2:       xb 8223
# 3:       ad 8179
# 4:       nc 8176
# 5:       gd 8173

所需的最终数据集

接下来，我想创建一个包含Age、Gender 和Date 变量的数据集。但我需要从demoDat 和timeDat 数据集中维护我的Val 边际总和。

我有一个策略可以完成这项任务，但它占用了相当多的 RAM。我可以采用另一种策略来一次在每个组内执行扩展吗？也许使用 .EACHI?

扩展两个数据集并合并

这是操作中昂贵的部分。数据集已扩展，因此行数等于sum(Val)。如果sum(Val) 是> 500,000,000，这可能会很昂贵。特别是因为对第二个数据集重复该操作。我希望使用.EACHI，这样只会扩展组内的数据，这应该会大大降低内存占用。

library(pryr)
memUsed <- mem_used() 
demoDatBig <- demoDat[rep(1:.N, Val), .(Location, Gender, Age, ID=rowid(Location))]
timeDatBig <- timeDat[rep(1:.N, Val), .(Location, Date, ID=rowid(Location))]
demoDatBig[timeDatBig, Date := i.Date, on=.(Location, ID)]
finalBigDat <- demoDatBig[, .(Val=.N), by=.(Location, Gender, Age, Date)]
mem_used() - memUsed
# 47 MB

所以这个操作占用了 47 MB 的 RAM，但是如果我增加 meanVal，它会显着增加。我希望在最大的Location 和ID 组上使用与此操作相同的功能所需的尽可能多的 RAM。我认为使用.EACHI 可以做到这一点，但我不确定如何。

结果数据表

       Location Gender Age       Date Val
    1:       aa      F   0 2016-01-01  36
    2:       aa      F   1 2016-01-01  47
    3:       aa      F   2 2016-01-01  29
    4:       aa      F   3 2016-01-01  40
    5:       aa      F   4 2016-01-01  24
   ---                                   
32430:       ze      M  96 2016-12-16   7
32431:       ze      M  96 2016-12-23  34
32432:       ze      M  97 2016-12-23  45
32433:       ze      M  98 2016-12-23  38
32434:       ze      M  99 2016-12-23  39

该解决方案有望通过这些测试

#### Test 1
test1 <- finalBigDat[, .(Val = sum(Val)), by=.(Location, Gender, Age)]
test1[demoDat, ValCheck := i.Val, on=.(Location, Gender, Age)]
test1[Val != ValCheck]
#Empty data.table (0 rows) of 5 cols: Location,Gender,Age,Val,ValCheck

#### Test 2
test2 <- finalBigDat[, .(Val = sum(Val)), by=.(Location, Date)]
test2[timeDat, ValCheck := i.Val, on=.(Location, Date)]
test2[Val != ValCheck]
#Empty data.table (0 rows) of 4 cols: Location,Date,Val,ValCheck

结果

我浏览了这两种解决方案并跟踪了两者的内存和系统时序。这两种解决方案都很棒，并且是对我目前拥有的解决方案的巨大升级。 @swihart 的解决方案可以难以置信地扩展到大 meanVal，所以我选择了这个作为接受的答案。当meanVal 没有那么大时，Heather 的回答会有所帮助。大大小小的meanVal 都经常出现，所以我都需要。

   meanVal            Ans            Time      Mem    Rows
1:      40     Mike.Gahan  0.6245470 secs 44.54293   32434
2:      40 Heather Turner  0.6391492 secs 38.65355 1352000
3:      40        swihart 11.1602619 secs 66.97550 1352000
4:     400     Mike.Gahan  2.593275 secs 437.23832   32611
5:     400 Heather Turner  1.303993 secs  38.79871 1352000
6:     400        swihart 11.736836 secs  66.97550 1352000
7:    4000     Mike.Gahan 30.390986 secs 4364.51501   32629
8:    4000 Heather Turner  6.279249 secs   38.79871 1352000
9:    4000        swihart 11.427965 secs   66.97550 1352000
10:   20000     Mike.Gahan -------did not finish----------
11:   20000 Heather Turner 23.78948 secs 36.30617 1352000
12:   20000        swihart 11.53811 secs 66.97550 1352000
13:   30000     Mike.Gahan -------did not finish----------
14:   30000 Heather Turner 537.6459  secs 57.15375 1352000
15:   30000        swihart 11.970013 secs 66.97474 1352000

【问题讨论】：

请在“占用相当多的内存”上填写数字，您希望提高多少？ 10%？ 90%？特别是请隔离哪一行中的哪个表达式正在浪费 RAM。请参阅有关内存分析的现有问题。如果需要，您可以随时使用 HDFS 等磁盘支持的存储。
meanVal 应该是 40，而不是 400，对吗？
是的..每个位置总是有 2 个性别级别和 100 个年龄级别。
没有。时间点不会遵守网格。
您对上面的语言不清楚 - 占用 RAM 的不是进程（本身） - 最终对象正在占用 RAM。 ?mem_used 显示 R 当前正在使用的内存，而不是给定进程使用的最大内存量。这也可以通过在您的“大”data.tables 上使用object.size 来验证 - 您会看到它们正在占用您的内存。那么 - 您是想让创建行为更高效，还是希望生成的对象占用更少的 RAM？

标签： r performance memory data.table

【解决方案1】：

我针对不同大小的meanVal 运行了您的方法，并看到了生成demoDatBig 和timeDatBig 的方法的缩放问题。我有一种方法（附在这篇文章的底部），它生成cartDat——日期和性别年龄组的笛卡尔交叉，对于meanVal和sum(Val)的增加是稳健的，如下表所示列出了正在讨论的 data.tables 的 object.size() 的结果：

| meanVal  | sum(Val) | demoDatBig (MB)  | timeDatBig (MB)  | cartDat (MB)  |
|----------|----------|------------------|------------------|---------------|
|      40  |     1e6  |            27.8  |            15.9  |          67.1 |
|     400  |     1e7  |           277.6  |           158.7  |          67.1 |
|   4,000  |     1e8  |         2,776.8  |         1,586.8  |          67.1 |
|  40,000  |     1e9  |        27,770.3  |        15,868.7  |          67.1 |

我的方法的关键是在未展开的源数据表demoDat 和timeDat 之间生成笛卡尔交叉，然后使用“迭代多元超几何采样”(IMHS) 方案来保留两个源数据表。为了让IMHS 具有 R 功能，我从 CRAN 中取出 R 包 BiasedUrn 并重新编译它，以便它可以处理 52 种颜色（在我们的应用程序中，日期）。如果需要调整给定位置的最大日期数，请告诉我，我将重新编译。因此，R package BiasedUrn52 is on github.

我的解决方案通过了test1 和test2 并保留了边缘。然而，与 OP 程序相比，它似乎在更多的日期中分配了性别年龄边缘。请允许我详细说明：

如果我们取timeDat的前5行：

> head(demoDat,5)
   Location Gender Age Val
1:       aa      F   0  36
2:       aa      F   1  47
3:       aa      F   2  29
4:       aa      F   3  40
5:       aa      F   4  50

还有finalBigDat的前6位：

> head(finalBigDat,6)
   Location Gender Age       Date Val
1:       aa      F   0 2016-01-01  36
2:       aa      F   1 2016-01-01  47
3:       aa      F   2 2016-01-01  29
4:       aa      F   3 2016-01-01  40
5:       aa      F   4 2016-01-01  24
6:       aa      F   4 2016-01-08  26

我们看到 F-0 性别年龄组的全部 36 个被归于 2016-01-01，而 F-4 组的 50 个分布在 2016-01-01 (24) 和 2016-01 -08 (26)，但没有其他日期 (50=24+26)。

IMHS 方法在更多日期之间分配边际（我不确定这是否需要 - 请告诉我）。例如，IMHS 占据了 F-0 组的 36 个，而不是像 finalBigDat 那样将所有 36 个放在 2016-01-01 上，而是将它们分散到更多的日期（查看 seq.Draws）：

> cartDat[Location=='aa' & Gender=="F" & Age==0,
+         c("Location", "Gender", "Age", "Date", "seq.Draws"),
+         with=FALSE]
    Location Gender Age       Date seq.Draws
 1:       aa      F   0 2016-01-01         1
 2:       aa      F   0 2016-01-08         0
 3:       aa      F   0 2016-01-15         1
 4:       aa      F   0 2016-01-22         1
 5:       aa      F   0 2016-01-29         0
 6:       aa      F   0 2016-02-05         0
 7:       aa      F   0 2016-02-12         0
 8:       aa      F   0 2016-02-19         0
 9:       aa      F   0 2016-02-26         0
10:       aa      F   0 2016-03-04         0
11:       aa      F   0 2016-03-11         0
12:       aa      F   0 2016-03-18         0
13:       aa      F   0 2016-03-25         3
14:       aa      F   0 2016-04-01         1
15:       aa      F   0 2016-04-08         0
16:       aa      F   0 2016-04-15         0
17:       aa      F   0 2016-04-22         1
18:       aa      F   0 2016-04-29         1
19:       aa      F   0 2016-05-06         0
20:       aa      F   0 2016-05-13         2
21:       aa      F   0 2016-05-20         0
22:       aa      F   0 2016-05-27         0
23:       aa      F   0 2016-06-03         0
24:       aa      F   0 2016-06-10         0
25:       aa      F   0 2016-06-17         1
26:       aa      F   0 2016-06-24         2
27:       aa      F   0 2016-07-01         0
28:       aa      F   0 2016-07-08         0
29:       aa      F   0 2016-07-15         0
30:       aa      F   0 2016-07-22         1
31:       aa      F   0 2016-07-29         0
32:       aa      F   0 2016-08-05         1
33:       aa      F   0 2016-08-12         1
34:       aa      F   0 2016-08-19         1
35:       aa      F   0 2016-08-26         1
36:       aa      F   0 2016-09-02         1
37:       aa      F   0 2016-09-09         2
38:       aa      F   0 2016-09-16         0
39:       aa      F   0 2016-09-23         1
40:       aa      F   0 2016-09-30         0
41:       aa      F   0 2016-10-07         2
42:       aa      F   0 2016-10-14         3
43:       aa      F   0 2016-10-21         0
44:       aa      F   0 2016-10-28         1
45:       aa      F   0 2016-11-04         1
46:       aa      F   0 2016-11-11         1
47:       aa      F   0 2016-11-18         0
48:       aa      F   0 2016-11-25         0
49:       aa      F   0 2016-12-02         2
50:       aa      F   0 2016-12-09         1
51:       aa      F   0 2016-12-16         1
52:       aa      F   0 2016-12-23         1

OP 方法和IMHS cartDat 方法之间的分布差异只是一个旁白。边缘被保留，如下所示。

timeDat 的边缘被保留：

> cartDat[, sum(seq.Draws), by=.(Location, Date)]
      Location       Date  V1
   1:       aa 2016-01-01 176
   2:       aa 2016-01-08 143
   3:       aa 2016-01-15 143
   4:       aa 2016-01-22 154
   5:       aa 2016-01-29 174
  ---                        
6756:       ze 2016-11-25 169
6757:       ze 2016-12-02 148
6758:       ze 2016-12-09 165
6759:       ze 2016-12-16 142
6760:       ze 2016-12-23 156
> timeDat
      Location       Date Val
   1:       aa 2016-01-01 176
   2:       aa 2016-01-08 143
   3:       aa 2016-01-15 143
   4:       aa 2016-01-22 154
   5:       aa 2016-01-29 174
  ---                        
6756:       ze 2016-11-25 169
6757:       ze 2016-12-02 148
6758:       ze 2016-12-09 165
6759:       ze 2016-12-16 142
6760:       ze 2016-12-23 156

demoDat 的边缘也是如此：

> cartDat[, sum(seq.Draws), by=.(Location, Gender, Age)]
       Location Gender Age V1
    1:       aa      F   0 36
    2:       aa      F   1 47
    3:       aa      F   2 29
    4:       aa      F   3 40
    5:       aa      F   4 50
   ---                       
25996:       ze      M  95 48
25997:       ze      M  96 41
25998:       ze      M  97 45
25999:       ze      M  98 38
26000:       ze      M  99 39
> demoDat
       Location Gender Age Val
    1:       aa      F   0  36
    2:       aa      F   1  47
    3:       aa      F   2  29
    4:       aa      F   3  40
    5:       aa      F   4  50
   ---                        
25996:       ze      M  95  48
25997:       ze      M  96  41
25998:       ze      M  97  45
25999:       ze      M  98  38
26000:       ze      M  99  39

这里是IMHS cartDat 方法和一些测试：

#Cartesian cross of demoDat and timeDat
devtools::install_github("swihart/BiasedUrn52")
library(BiasedUrn52)
setkey(timeDat, Location)
setkey(demoDat, Location, Gender, Age)
cartDat <- demoDat[timeDat, allow.cartesian=TRUE]
setkeyv(cartDat, key(demoDat))
cartDat
cartDat[,group:=.GRP,by=c("Gender", "Age") ]
cartDat[,demoDat.Val:=Val]
cartDat[,timeDat.Val:=i.Val]
setcolorder(cartDat, c("Location", 
                       "group",
                       "Gender",
                       "Age",
                       "Val",
                       "demoDat.Val",
                       "Date",
                       "timeDat.Val",
                       "i.Val"))

#Define Iterative Multivariate Hypergeometric Sampling function
imhs <- function(.N, Val, i.Val, group){

  grp.ind <- unique(group)
  num.grp <- max(group)
  grp.size <- as.numeric(table(group))

  draws <- rep(NA, length(group))
  for(grp in grp.ind){

    if(grp==1){
      draws[group==1] = rMFNCHypergeo(1, 
                                      i.Val[group==1], 
                                      Val[group==1][1], 
                                      rep(1/grp.size[grp.ind==1],grp.size[grp.ind==1])
      )
      i.Val[group==2]= i.Val[group==1]-draws[group==1]
    }else{
      draws[group==grp] = rMFNCHypergeo(1, 
                                        i.Val[group==grp], 
                                        Val[group==grp][1], 
                                        rep(1/grp.size[grp.ind==grp],grp.size[grp.ind==grp])
      )
      if(grp<=num.grp){
        i.Val[group==(grp+1)]= i.Val[group==grp]-draws[group==grp]
      }
    }

  }

  list(i.Val, draws)
}


# run it the data.table way:
cartDat[,
        c("seq.Val", "seq.Draws") := imhs(.N, demoDat.Val, timeDat.Val, group),        
        by=c("Location") ]

# take a look:
cartDat

# reconciliation
demoDat[, sum(Val), by=.(Location)][order(-V1)]
cartDat[, sum(seq.Draws), by=.(Location)][order(-V1)]

# do the checks for the margins:
cartDat[, sum(seq.Draws), by=.(Location, Date)]
timeDat
cartDat[, sum(seq.Draws), by=.(Location, Gender, Age)]
demoDat


# such different sizes due to distributing across more dates:
nrow(demoDat)
nrow(cartDat)
nrow(cartDat[seq.Draws != 0])
nrow(finalBigDat)
nrow(cartDat[seq.Draws != 0])/nrow(finalBigDat)

# attain and print object sizes for cartDat
print(object.size(cartDat), units = "Mb")
print(object.size(cartDat[seq.Draws!=0]), units="Mb")

# attain and print object sizes for demoDatBig, timeDatBig, finalBigData
print(object.size(demoDatBig), units = "Mb")
print(object.size(timeDatBig), units = "Mb")
print(object.size(finalBigDat), units = "Mb")



## (OP) The solution would pass these tests:
finalBigDat2 <- cartDat

#### Test 1 (change to sum(seq.Draws))
test1 <- finalBigDat2[, .(Val = sum(seq.Draws)), by=.(Location, Gender, Age)]
test1[demoDat, ValCheck := i.Val, on=.(Location, Gender, Age)]
test1[Val != ValCheck]
#Empty data.table (0 rows) of 5 cols: Location,Gender,Age,Val,ValCheck

#### Test 2 (change to sum(seq.Draws))
test2 <- finalBigDat2[, .(Val = sum(seq.Draws)), by=.(Location, Date)]
test2[timeDat, ValCheck := i.Val, on=.(Location, Date)]
test2[Val != ValCheck]
#Empty data.table (0 rows) of 4 cols: Location,Date,Val,ValCheck

【讨论】：

是否有可能让它在 104 周而不是 52 周内工作？您是否将 Makevars 文件更改为 52，因为它会更快？
我将 makevars 从 30 增加到 52。我可以增加到 104。我会在这个周末的某个时候做。
不用担心。我已经更改并重新编译。我只是想知道那个最大值是否存在是有原因的？如果我将该值增加太多会不会很危险？我不确定。
我也想知道。可能会在描述中找到 Agner Fog 的电子邮件。
刚刚给 Agner 发了邮件。最大颜色设置高一点没问题，但是如果颜色数量多，有些功能会花费很长时间。

【解决方案2】：

通过加入位置和 ID，您当前的方法可以系统地将日期与性别：年龄类别中的观察结果相匹配，例如将第一个日期分配给前 176 个观察值（涵盖前四个性别：年龄类别和第五个类别的一部分）等等。

相反，您可以通过使用 stats 包中的 r2dtable 函数模拟双向 (gender:age):date 列联表来模拟每个性别：年龄类别中每个日期的计数这样边际总数是固定的。

首先在每个位置创建性别、年龄和日期的组合

setkey(timeDat, Location)
setkey(demoDat, Location)
finalBigDat <- demoDat[timeDat, .(Location, Gender, Age, Date),
                       allow.cartesian=TRUE]

然后在每个Location内应用r2dtable，将行总计设置为每个性别：年龄类别的计数，将列总计设置为每个日期的计数：

setkey(finalBigDat, Location)
finalBigDat[, Val := c(r2dtable(1, demoDat[.BY, Val], timeDat[.BY, Val])[[1]]),
            by = Location]
head(finalBigDat)
#    Location Gender Age       Date Val
# 1:       aa      F   0 2016-01-01   0
# 2:       aa      F   1 2016-01-01   3
# 3:       aa      F   2 2016-01-01   2
# 4:       aa      F   3 2016-01-01   2
# 5:       aa      F   4 2016-01-01   1
# 6:       aa      F   5 2016-01-01   0

这将根据需要通过测试 1 和 2。无论meanVal 的值是多少，finalBigDat 的尺寸总是相同的，但是rd2table 将需要更长的时间来运行更高的meanVal。以下是一些说明性的时间安排：

| meanVal  | sum(Val) | time (s) |
|----------|----------|----------|
|      40  |     1e6  |     0.36 | 
|    4000  |     1e8  |     7.72 |
|   20000  |     5e8  |    35.09 |

我确实尝试了 meanVal = 40000，但在大约 300 秒时停止了代码运行，所以如果你想设置一个高的 meanVal，你可能会遇到困难。

【讨论】：