在 R 中有效抽样彩票号码
Efficiently Sampling Lottery Numbers in R
我想编写一个函数，对 n 张彩票进行采样，每个彩票的 6 个号码从 1 到 45没有替换。但是，我需要高效地执行此操作，这意味着没有循环或类似循环的函数。 （我想 Rcpp 也可以，但我更喜欢基于 R 的矢量化解决方案）
无限制求解：
lottery_inef <- function(n){
t(replicate(n,
sample(1:45, 6)))
}
所以这里我得到一个矩阵，其中每一行对应一张彩票。现在，如果我想模拟数百万张彩票，这会变得非常慢，因此我对矢量化解决方案很感兴趣。
我的想法是：
lottery_ef <- function(n){
m <- matrix(sample(1:45, n*6, replace = TRUE), ncol = 6)
# somehow subset the matrix without a loop to remove all the
# rows that have non-unique values as in the lottery we can only draw each number once
}
对于高效版本，我在没有循环或 apply() 的子集化点上有点迷茫。如果有人能解决这个子集化问题或指出一个完全不同的方向来引导我找到解决方案，我将不胜感激。
replicate 实际上并没有在这个规模上做得那么好。使用即时编译（现在在 R 中使用了几年），for 循环可以更快，尤其是当我们可以精确地 pre-allocate 数据结构时。我们也可以避免 t():
lottery_inef <- function(n){
t(replicate(n,
sample(1:45, 6)))
}
lottery_preall <- function(n){
m = matrix(NA_integer_, nrow = n, ncol = 6)
for(i in 1:n) {
m[i, ] = sample.int(45L, size = 6)
}
m
}
nn = 1e6
microbenchmark::microbenchmark(
lottery_inef(nn),
lottery_preall(nn),
times = 2
)
# Unit: seconds
# expr min lq mean median uq max neval
# lottery_inef(nn) 9.400862 9.400862 9.571756 9.571756 9.742649 9.742649 2
# lottery_preall(nn) 4.948216 4.948216 5.454482 5.454482 5.960749 5.960749 2
replicate 将结果累加到 list 中，然后需要检查每个结果的维度，然后再决定是否可以将其简化为矩阵，并且必须进行转换。所有这些开销都被一个 pre-allocated 整数矩阵跳过了大约 2x speed-up.
我们也可以比较，比如 vapply（快速测试显示 vapply 比循环慢一点），但我认为要从中获得更快的速度，你会需要并行 运行 - 这在这里是一个不错的选择，可能会让你 speed-up 几乎等于你使用的内核数量。
sample.int 几乎只是对 C 代码的调用，所以使用 Rcpp 可能不会做得更好——我认为并行化是提高速度的最佳选择。
由于为这个大小的集合生成所有组合只需要几秒钟，因此可能值得这样做，然后将其子集用于 'lottery tickets'。下面我使用 sample() 生成了 100 万行索引（包括替换和不替换），并对整个集合进行了括号子集化以生成可能的票证。
如果您需要经常这样做，或者在不同的时间这样做，可能值得保存完整的组合集而不是每次 re-generate。几乎所有的处理都是生成完整的组合集。之后选择 'tickets' 很快。
时间显示创建所有组合花费了约 6 秒，创建 100 万个索引花费了约 .2 秒，对 100 万行的括号子集创建了约 .1 秒。
set.seed(2)
tictoc::tic() #included for timing
# All possible lotto combinations as matrix, 1 per row
lotto_all <- t(combn(1:45, 6))
tictoc::toc() #included for timing
#> 5.899 sec elapsed
# A look at the data:
head(lotto_all)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1 2 3 4 5 6
#> [2,] 1 2 3 4 5 7
#> [3,] 1 2 3 4 5 8
#> [4,] 1 2 3 4 5 9
#> [5,] 1 2 3 4 5 10
#> [6,] 1 2 3 4 5 11
# Getting index (row) numbers for our 'tickts' with & without replacement
tictoc::tic()
sample_indices_no_replacement <- sample(1:nrow(lotto_all), size = 1e6, replace = F)
tictoc::toc()
#> 0.178 sec elapsed
sample_indices_w_replacement <- sample(1:nrow(lotto_all), size = 1e6, replace = T)
# The number combinations of our 'tickets'
tictoc::tic()
sample_tickets_no_rep <- lotto_all[sample_indices_no_replacement,]
tictoc::toc()
#> 0.097 sec elapsed
sample_tickets_rep <- lotto_all[sample_indices_w_replacement,]
# A look at the sample tickets:
head(sample_tickets_no_rep)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 8 12 14 31 34 44
#> [2,] 6 10 16 26 32 36
#> [3,] 3 4 10 15 41 43
#> [4,] 2 3 5 17 33 36
#> [5,] 7 17 24 25 35 40
#> [6,] 32 33 34 36 39 43
# See that there are some duplicates using replacement = T
length(unique(sample_indices_no_replacement))
#> [1] 1000000
length(unique(sample_indices_w_replacement))
#> [1] 941309
由 reprex package (v0.3.0)
于 2020-10-27 创建
由于效率是主要关注点，因此有几个包，arrangements 和 RcppAlgos*, 记住了这个特定任务。
在开始之前，我们首先说明当我们使用sample时，我们无法控制结果的唯一性。每次绘制都来自均匀分布，因此我们可以多次重复绘制相同的排列。使用@Gregor 的函数，我们有：
set.seed(42)
system.time(a <- lottery_inef(1e6))
user system elapsed
7.640 0.345 7.984
sum(duplicated(a))
[1] 86
set.seed(42)
system.time(b <- lottery_preall(1e6))
user system elapsed
3.673 0.256 3.929
sum(duplicated(b))
[1] 86
虽然使用包 arrangements 速度更快，但我们仍然看到相同的行为：
set.seed(42)
system.time(arng <- arrangements::permutations(45, 6, nsample = 1e6))
user system elapsed
0.761 0.021 0.785
sum(duplicated(arng))
[1] 108
现在，如果请求的结果数量少于结果总数（在我们的例子中超过 50 亿），使用包 RcppAlgos 我们保证结果是唯一的：
RcppAlgos::permuteCount(45, 6)
[1] 5864443200
system.time(algosSer <- RcppAlgos::permuteSample(45, 6,
n = 1e6,
seed = 42))
user system elapsed
0.560 0.001 0.561
sum(duplicated(algosSer))
[1] 0
此外，我们可以通过 nThreads 参数利用多个线程来实现更快的速度。
system.time(algosPar <- RcppAlgos::permuteSample(45, 6,
n = 1e6,
seed = 42,
nThreads = 4))
user system elapsed
0.574 0.001 0.280
## Results are the same as the serial version
identical(algosPar, algosSer)
[1] TRUE
* 我是 RcppAlgos
的作者
我想编写一个函数，对 n 张彩票进行采样，每个彩票的 6 个号码从 1 到 45没有替换。但是，我需要高效地执行此操作，这意味着没有循环或类似循环的函数。 （我想 Rcpp 也可以，但我更喜欢基于 R 的矢量化解决方案）
无限制求解：
lottery_inef <- function(n){
t(replicate(n,
sample(1:45, 6)))
}
所以这里我得到一个矩阵，其中每一行对应一张彩票。现在，如果我想模拟数百万张彩票，这会变得非常慢，因此我对矢量化解决方案很感兴趣。
我的想法是：
lottery_ef <- function(n){
m <- matrix(sample(1:45, n*6, replace = TRUE), ncol = 6)
# somehow subset the matrix without a loop to remove all the
# rows that have non-unique values as in the lottery we can only draw each number once
}
对于高效版本，我在没有循环或 apply() 的子集化点上有点迷茫。如果有人能解决这个子集化问题或指出一个完全不同的方向来引导我找到解决方案，我将不胜感激。
replicate 实际上并没有在这个规模上做得那么好。使用即时编译（现在在 R 中使用了几年），for 循环可以更快，尤其是当我们可以精确地 pre-allocate 数据结构时。我们也可以避免 t():
lottery_inef <- function(n){
t(replicate(n,
sample(1:45, 6)))
}
lottery_preall <- function(n){
m = matrix(NA_integer_, nrow = n, ncol = 6)
for(i in 1:n) {
m[i, ] = sample.int(45L, size = 6)
}
m
}
nn = 1e6
microbenchmark::microbenchmark(
lottery_inef(nn),
lottery_preall(nn),
times = 2
)
# Unit: seconds
# expr min lq mean median uq max neval
# lottery_inef(nn) 9.400862 9.400862 9.571756 9.571756 9.742649 9.742649 2
# lottery_preall(nn) 4.948216 4.948216 5.454482 5.454482 5.960749 5.960749 2
replicate 将结果累加到 list 中，然后需要检查每个结果的维度，然后再决定是否可以将其简化为矩阵，并且必须进行转换。所有这些开销都被一个 pre-allocated 整数矩阵跳过了大约 2x speed-up.
我们也可以比较，比如 vapply（快速测试显示 vapply 比循环慢一点），但我认为要从中获得更快的速度，你会需要并行 运行 - 这在这里是一个不错的选择，可能会让你 speed-up 几乎等于你使用的内核数量。
sample.int 几乎只是对 C 代码的调用，所以使用 Rcpp 可能不会做得更好——我认为并行化是提高速度的最佳选择。
由于为这个大小的集合生成所有组合只需要几秒钟，因此可能值得这样做，然后将其子集用于 'lottery tickets'。下面我使用 sample() 生成了 100 万行索引（包括替换和不替换），并对整个集合进行了括号子集化以生成可能的票证。
如果您需要经常这样做，或者在不同的时间这样做，可能值得保存完整的组合集而不是每次 re-generate。几乎所有的处理都是生成完整的组合集。之后选择 'tickets' 很快。
时间显示创建所有组合花费了约 6 秒，创建 100 万个索引花费了约 .2 秒，对 100 万行的括号子集创建了约 .1 秒。
set.seed(2)
tictoc::tic() #included for timing
# All possible lotto combinations as matrix, 1 per row
lotto_all <- t(combn(1:45, 6))
tictoc::toc() #included for timing
#> 5.899 sec elapsed
# A look at the data:
head(lotto_all)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 1 2 3 4 5 6
#> [2,] 1 2 3 4 5 7
#> [3,] 1 2 3 4 5 8
#> [4,] 1 2 3 4 5 9
#> [5,] 1 2 3 4 5 10
#> [6,] 1 2 3 4 5 11
# Getting index (row) numbers for our 'tickts' with & without replacement
tictoc::tic()
sample_indices_no_replacement <- sample(1:nrow(lotto_all), size = 1e6, replace = F)
tictoc::toc()
#> 0.178 sec elapsed
sample_indices_w_replacement <- sample(1:nrow(lotto_all), size = 1e6, replace = T)
# The number combinations of our 'tickets'
tictoc::tic()
sample_tickets_no_rep <- lotto_all[sample_indices_no_replacement,]
tictoc::toc()
#> 0.097 sec elapsed
sample_tickets_rep <- lotto_all[sample_indices_w_replacement,]
# A look at the sample tickets:
head(sample_tickets_no_rep)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 8 12 14 31 34 44
#> [2,] 6 10 16 26 32 36
#> [3,] 3 4 10 15 41 43
#> [4,] 2 3 5 17 33 36
#> [5,] 7 17 24 25 35 40
#> [6,] 32 33 34 36 39 43
# See that there are some duplicates using replacement = T
length(unique(sample_indices_no_replacement))
#> [1] 1000000
length(unique(sample_indices_w_replacement))
#> [1] 941309
由 reprex package (v0.3.0)
于 2020-10-27 创建由于效率是主要关注点，因此有几个包，arrangements 和 RcppAlgos*, 记住了这个特定任务。
在开始之前，我们首先说明当我们使用sample时，我们无法控制结果的唯一性。每次绘制都来自均匀分布，因此我们可以多次重复绘制相同的排列。使用@Gregor 的函数，我们有：
set.seed(42)
system.time(a <- lottery_inef(1e6))
user system elapsed
7.640 0.345 7.984
sum(duplicated(a))
[1] 86
set.seed(42)
system.time(b <- lottery_preall(1e6))
user system elapsed
3.673 0.256 3.929
sum(duplicated(b))
[1] 86
虽然使用包 arrangements 速度更快，但我们仍然看到相同的行为：
set.seed(42)
system.time(arng <- arrangements::permutations(45, 6, nsample = 1e6))
user system elapsed
0.761 0.021 0.785
sum(duplicated(arng))
[1] 108
现在，如果请求的结果数量少于结果总数（在我们的例子中超过 50 亿），使用包 RcppAlgos 我们保证结果是唯一的：
RcppAlgos::permuteCount(45, 6)
[1] 5864443200
system.time(algosSer <- RcppAlgos::permuteSample(45, 6,
n = 1e6,
seed = 42))
user system elapsed
0.560 0.001 0.561
sum(duplicated(algosSer))
[1] 0
此外，我们可以通过 nThreads 参数利用多个线程来实现更快的速度。
system.time(algosPar <- RcppAlgos::permuteSample(45, 6,
n = 1e6,
seed = 42,
nThreads = 4))
user system elapsed
0.574 0.001 0.280
## Results are the same as the serial version
identical(algosPar, algosSer)
[1] TRUE
* 我是 RcppAlgos