在R中的大数据集中按行生成多项随机变量需要加快速度

Question

我有以下结构的数据，我正在尝试实现一个函数，对于每一行，使用随机多项式生成器分配一个选项，其值为<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包rMultinom的函数Hmisc生成具有不同概率的多项式随机变量。

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|

Answer 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

Answer 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秒才能在我的电脑上运行。希望这可以帮助！

Answer 3

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

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

所有函数都是R的内置函数，probs是列式的，在某种意义上，probs的一列标识了一个多项式绘制的概率。结果是一个矩阵。每列都有结果（结果类为1，其他类为0），结果快60％