r - 计算 36 个元素的向量的 N 个随机排列

Question

我有一个包含 36 个元素的向量 V，18 个是“0”，18 个是“1”。 我想计算这个向量的 N 个随机(不是第一个 N)排列。

我可以这样做:

library(combinat)
N <- 100 # or 200, 300, 500... max 1000
V <- c(rep(0, 18), rep(1, 18))
n <- factorial(36) # total number of unique possible permutations
p <- unique(permn(V))[sample(1:n, N)]

但我很快遇到了组合爆炸问题，因为 sample(1:n, N) 返回 Error in 1:n : result would be too long a vector

和

permn(V) 返回 向量错误("list", gamma(n + 1)) : vector size specified is too large

还有其他(更好的)方法可以做到这一点吗？

Answer 1

首先，没有 factorial(36) 结果，因为您有重复的元素。如果我们这样做了，要获取总数，我们可以使用 gmp 包来获取:

gmp::factorialZ(36)
Big Integer ('bigz') :
[1] 371993326789901217467999448150835200000000

我们实际处理的称为 multisets 的排列(正如@JakubBucek 在评论中指出的那样)。使用 RcppAlgos 包(我编写的)或 arrangements 包，我们可以轻松正确地计算结果总数，更重要的是生成所需的结果。

首先，实际结果数:

arrangements::npermutations(0:1, freq = c(18, 18), bigz = TRUE)
Big Integer ('bigz') :
[1] 9075135300

RcppAlgos::permuteCount(0:1, freqs = c(18, 18))
[1] 9075135300

这是组合学的结果。也就是说，我们必须除以相似元素的排列数的乘积:

gmp::factorialZ(36) / gmp::pow.bigz(gmp::factorialZ(18), 2)
Big Rational ('bigq') :
[1] 9075135300

现在，生成随机排列。对于包 arrangements，我们使用 nsample 参数。此外，我们可以设置可重复性的种子:

set.seed(123)
r1 <- arrangements::permutations(0:1, freq = c(18, 18), nsample = 10)

set.seed(123)
r2 <- arrangements::permutations(0:1, freq = c(18, 18), nsample = 10)

dim(r1)
[1] 10 36

identical(r1, r2)
[1] TRUE

## only showing 10 columns
head(r1[,1:10])
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    0    0    0    0    1    1    0    1    1     1
[2,]    1    0    1    1    1    1    1    1    1     0
[3,]    0    0    0    0    0    1    1    0    0     0
[4,]    1    1    1    0    0    1    0    1    0     0
[5,]    0    1    1    0    0    1    1    1    0     1
[6,]    0    0    0    1    1    1    0    1    1     1

对于 RcppAlgos，我们使用类似的参数 n 和 seed 调用 permuteSample:

r3 <- RcppAlgos::permuteSample(0:1, freqs = c(18, 18), n = 10, seed = 42)
r4 <- RcppAlgos::permuteSample(0:1, freqs = c(18, 18), n = 10, seed = 42)

identical(r3, r4)
[1] TRUE

dim(r3)
[1] 10 36

这两个软件包也非常有效。生成 1000 个随机排列只需不到一秒的时间:

system.time(RcppAlgos::permuteSample(0:1, freqs = c(18, 18), n = 1000))
 user  system elapsed 
0.051   0.000   0.052 

system.time(arrangements::permutations(0:1, freq = c(18, 18), nsample = 1000))
 user  system elapsed 
0.249   0.000   0.249