在R中的大数据集中按行生成多项式随机变量。需要加快速度
Generating multinomial random variables by row in big dataset in R. Need to speed up
我有具有以下结构的数据，我正在尝试实现一个函数，该函数为每一行分配一个选择，使用随机多项式生成器将值 <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|
您可以在没有 rMultinom 函数的情况下执行此操作，方法是生成一个随机统一变量并检查它位于哪个区间，如 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 秒。希望这对您有所帮助！
与其调用 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
我在这里寻找同样问题的解决方案。老实说我没有找到它，但也许我可以给你一个更好的解决方案。
z = mapply(rmultinom, n = 1, size = 1, prob = split(probs, c(col(probs))))
所有函数都是 R 的内置函数，probs 是按列排列的，在某种意义上，probs 的一列标识一次抽取多项式的概率。结果是一个矩阵。每列都有结果（结果的 class 为 1，其他为 0）并且快 60%
我有具有以下结构的数据，我正在尝试实现一个函数，该函数为每一行分配一个选择，使用随机多项式生成器将值 <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|
您可以在没有 rMultinom 函数的情况下执行此操作，方法是生成一个随机统一变量并检查它位于哪个区间，如 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 秒。希望这对您有所帮助！
与其调用 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
我在这里寻找同样问题的解决方案。老实说我没有找到它，但也许我可以给你一个更好的解决方案。
z = mapply(rmultinom, n = 1, size = 1, prob = split(probs, c(col(probs))))
所有函数都是 R 的内置函数，probs 是按列排列的，在某种意义上，probs 的一列标识一次抽取多项式的概率。结果是一个矩阵。每列都有结果（结果的 class 为 1，其他为 0）并且快 60%