在 R 中有效地采样彩票号码答案

【问题标题】：Efficiently Sampling Lottery Numbers in R在 R 中有效地采样彩票号码
【发布时间】：2020-10-27 17:42:15
【问题描述】：

我想编写一个函数，对 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() 的情况下，我有点迷失了。如果有人能解决这个子集问题或指出一个完全不同的方向来引导我找到解决方案，我将不胜感激。

【问题讨论】：

您能否将所有可能的组合生成为矩阵的行，然后仅对行号进行采样（在这种情况下为 1-~810 万）？生成完整集需要几秒钟，但索引子集很快。
@mrhellmann 好主意，我刚试过这个，如果我创建矩阵大约需要 9 秒，采样一百万行只需几毫秒。如果你想写这个作为答案，我很乐意接受。

标签： r function sampling

【解决方案1】：

replicate 在这个规模上实际上并没有做得那么好。使用即时编译（现在在 R 中使用了几年），for 循环可以更快，尤其是当我们可以精确地预分配数据结构时。我们也可以避开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 中，然后需要检查每个维度的维度，然后再决定是否可以将其简化为矩阵，并且必须进行转换。所有这些开销都被一个预先分配的整数矩阵跳过了，速度大约提高了 2 倍。

我们还可以与 vapply 进行比较（快速测试显示 vapply 只比循环慢一点），但我认为要从中获得更快的速度，您需要并行运行 -这将是一个不错的选择，并且可能可以让您获得几乎等于您使用的核心数量的加速。

sample.int 几乎只是对 C 代码的调用，因此使用 Rcpp 可能不会做得更好 - 我认为并行化是提高速度的最佳选择。

【讨论】：

【解决方案2】：

由于生成此大小的集合的所有组合只需要几秒钟，因此可能值得这样做，然后将其子集用于“彩票”。下面我使用sample() 生成了 100 万行索引（包括替换和不替换），并在整个集合上用括号括起来以生成可能的票证。

如果您需要经常执行此操作，或者在不同时间执行此操作，则可能值得保存完整的组合集，而不是每次都重新生成它。几乎所有的处理都在生成全套组合。之后选择“票”很快。

时间显示创建所有组合需要大约 6 秒，100 万个索引大约需要 0.2 秒，100 万行的括号子集大约需要 0.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 日创建}

【讨论】：

@dww 上面的版本有生成唯一票证和重复票证的方法，不是吗？注意最后 5 行代码。
你可以使用 comb_n 包中的函数 Rfast，这比 R 的基础快。

【解决方案3】：

由于效率是主要关注点，因此有几个软件包，arrangements 和 RcppAlgos^* 考虑到了这一特定任务。 p>

在开始之前，我们首先说明当我们使用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

现在，使用 RcppAlgos 包，如果请求的结果数量少于结果总数（在我们的例子中超过 50 亿），我们可以保证结果是唯一的：

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的作者

【讨论】：