在R中的大数据集中按行生成多项式随机变量。需要加速答案

【问题标题】：Generating multinomial random variables by row in big dataset in R. Need to speed up在R中的大数据集中按行生成多项式随机变量。需要加速
【发布时间】：2018-02-08 20:46:02
【问题描述】：

我有具有以下结构的数据，我正在尝试实现一个函数，该函数使用随机多项式生成器为每一行分配一个选择，其中值<Educ_W1:Educ_W5> 作为概率向量（它们逐行相加） .因此，对于每一行，新变量的值都在 1 到 5 之间。我自己设法实现了它，但我试图找到一种更快的方法来实现它，因为在当前版本中它花费的时间太长（几天。数据包含超过 100 万个观察值）。

| IDhh|Year |Educ_W |Educ_H | Educ_W1| Educ_W2| Educ_W3| Educ_W4| Educ_W5| |----:|:----|:------|:------|---------:|---------:|---------:|---------:|---------:| | 1|1975 |2 |2 | 0.1645188| 0.3362659| 0.3940354| 0.0831637| 0.0220162| | 2|1975 |2 |2 | 0.1645188| 0.3362659| 0.3940354| 0.0831637| 0.0220162| | 5|1975 |2 |1 | 0.5103815| 0.2092249| 0.2285570| 0.0392398| 0.0125968| | 6|1975 |3 |3 | 0.0811203| 0.1535407| 0.5528233| 0.1486548| 0.0638609| | 8|1975 |1 |1 | 0.5103815| 0.2092249| 0.2285570| 0.0392398| 0.0125968| | 10|1975 |3 |2 | 0.1645188| 0.3362659| 0.3940354| 0.0831637| 0.0220162|

目前我正在通过以下方式实现该功能，但我需要很长时间。这里，变量“IDhh”唯一标识每一行。 R 包Hmisc 中的函数rMultinom 生成具有不同概率的多项随机变量。

library(dplyr)
library(tidyr)
data %>%
  select(IDhh, Year, Educ_W, Educ_H, Educ_W1 : Educ_W5) %>% 
  nest(-IDhh) %>% 
  mutate(
    wanted_W = map(data, ~ rMultinom(t(c(.x$Educ_W1, .x$Educ_W2, .x$Educ_W3, 
                                         .x$Educ_W4, .x$Educ_W5)), 1))) %>%
  unnest()

所需的输出如下所示，其中“Wanted_W”是新变量。

| IDhh| wanted_W|Year |Educ_W | Educ_W1| Educ_W2| Educ_W3| Educ_W4| Educ_W5| |-------:|--------:|:----|:------|---------:|---------:|---------:|---------:|---------:| | 18806| 3|1975 |3 | 0.1851884| 0.1577067| 0.4749609| 0.1394014| 0.0427427| | 2442099| 4|2010 |1 | 0.4436620| 0.0987973| 0.3296288| 0.1013606| 0.0265513| | 1351429| 3|1995 |3 | 0.0708855| 0.1023657| 0.5904598| 0.1784980| 0.0577910| | 250232| 3|1980 |5 | 0.0337913| 0.0347975| 0.2156134| 0.2315768| 0.4842209| | 1802868| 3|2005 |3 | 0.0371280| 0.0772428| 0.6054841| 0.2024385| 0.0777067| | 715077| 2|1985 |3 | 0.1149756| 0.1412112| 0.5458910| 0.1413975| 0.0565248|

【问题讨论】：

我建议您使用函数put 来创建可重现的示例。除此之外，我认为瓶颈是函数rMultinom。为了处理海量数据集，您需要可以矢量化的运算符。

标签： r random dplyr purrr multinomial

【解决方案1】：

而不是调用Hmisc::rMultinom 一百万次（数据中的每一行一次），您可以将概率参数作为矩阵提供给函数。然后矩阵中的每一行将定义一个不同的多项分布。

reprex::reprex_info()
#> Created by the reprex package v0.1.1.9000 on 2018-02-09

library(dplyr)
set.seed(1)

# Generate category probabilities
n <- 1e6
unifs <- replicate(5, runif(n))
probs <- sweep(unifs, 1, apply(unifs, 1, sum), "/")
colnames(probs) <- paste0("p", seq_len(ncol(probs)))

df <- as_tibble(probs)

system.time({
  probs <- as.matrix(df %>% select(p1:p5))
  res <- df %>% 
    mutate(rcat = Hmisc::rMultinom(probs, 1))
})
#>    user  system elapsed 
#>    9.25    0.15    9.50

res
#> # A tibble: 1,000,000 x 6
#>            p1         p2         p3         p4         p5  rcat
#>         <dbl>      <dbl>      <dbl>      <dbl>      <dbl> <chr>
#>  1 0.14607852 0.07709049 0.33798110 0.22154639 0.21730349    p4
#>  2 0.12813691 0.23952958 0.11025717 0.31642808 0.20564827    p4
#>  3 0.19137423 0.24349984 0.06855848 0.23421041 0.26235703    p3
#>  4 0.30227095 0.03050219 0.27667295 0.28389810 0.10665580    p3
#>  5 0.10096040 0.03334545 0.07350112 0.38768513 0.40450791    p4
#>  6 0.32430441 0.22123172 0.13317669 0.08001760 0.24126959    p2
#>  7 0.32710720 0.14134942 0.25371663 0.20344497 0.07438178    p1
#>  8 0.21841291 0.23480314 0.25563400 0.06838794 0.22276200    p3
#>  9 0.21164692 0.19809418 0.15415735 0.15095640 0.28514514    p1
#> 10 0.02220492 0.23105648 0.35661756 0.08688459 0.30323645    p3
#> # ... with 999,990 more rows

【讨论】：

【解决方案2】：

您可以在不使用 rMultinom 函数的情况下执行此操作，方法是生成一个随机 Uniform 变量并检查它所在的区间，如 here 所述：

set.seed(1234)
n=1000000
library(data.table)

# Sample data -----------------------------------------------------------
create_probs <- function(x)
{
  y = sample(1:10,x)
  y = as.list(y/sum(y))
  return(y)
}
p_dt = data.table(id=1:n)
p_dt =p_dt[,c("Educ_w1","Educ_w2","Educ_w3","Educ_w4","Educ_w5"):=create_probs(5),by=1:nrow(p_dt)]


# Function --------------------------------------------------------------
p_dt[,U:=runif(1,0,1),1:nrow(p_dt)]
p_dt = p_dt[,Educ_w:=min(which(cumsum(unlist(.SD))>U)),1:nrow(p_dt),
            .SDcols=c("Educ_w1","Educ_w2","Educ_w3","Educ_w4","Educ_w5")]

head(p_dt)

样本输出：

   id    Educ_w1    Educ_w2   Educ_w3    Educ_w4    Educ_w5          U Educ_w
1:  1 0.06666667 0.20000000 0.1666667 0.26666667 0.30000000 0.49320836      4
2:  2 0.36842105 0.05263158 0.1052632 0.26315789 0.21052632 0.54415445      4
3:  3 0.25925926 0.18518519 0.1111111 0.37037037 0.07407407 0.65840751      4
4:  4 0.29032258 0.09677419 0.3225806 0.06451613 0.22580645 0.26604797      1
5:  5 0.22222222 0.16666667 0.1111111 0.05555556 0.44444444 0.05887458      1
6:  6 0.31034483 0.17241379 0.2758621 0.20689655 0.03448276 0.98659704      5

这个功能部分在我的电脑上运行大约需要 8 秒。希望这会有所帮助！

【讨论】：

您好弗洛里安，非常感谢您的回答。我没有这样做的原因是许多人共享相同的概率向量。使用您的方法，它们最终都会得到相同的选择，而使用“rMultinom”，对于 Educ_W1:Educ_W5 中具有相同值的行，结果仍然是随机的。最好的！
糟糕，对不起！我误读了这个问题，并假设概率也是随机确定的。我已经用我希望更接近您正在寻找的解决方案更新了我的答案，尽管我看到还添加了另一个看起来非常简单和有希望的答案；）

【解决方案3】：

我在这里寻找相同问题的解决方案。老实说我没有找到它，但也许我可以给你一个更好的解决方案。

z = mapply(rmultinom, n = 1, size = 1, prob = split(probs, c(col(probs))))

所有函数都是 R 的内置函数，probs 是按列排列的，从某种意义上说，probs 的一列标识了一次多项式抽签的概率。结果是一个矩阵。每列都有结果（结果类为 1，其他类为 0），速度快 60%

【讨论】：