R中N个元素与q个元素的组合
Posted
技术标签:
【中文标题】R中N个元素与q个元素的组合【英文标题】：Combination of N elements with q elements in R 【发布时间】：2018-04-24 06:29:57 【问题描述】：我有N=6 元素和q=3 元素符号为0,1,2。
我想在所有可能的位置创建N=6元素的所有向量，其中2元素等于0、2元素等于1和2元素等于2 .
这些向量的数量等于combn(6,2)*combn(4,2)*combn(2,2)=90。
这是在矩阵F 中构造这些90 向量的代码：
N=6
x<-c(1:N)
#########################################
A<-combn(x,2)
B<-matrix(0,ncol(A),length(x))
for( i in 1:ncol(A) )
y<-rep(0,N)
y[A[1:nrow(A),i]]<-1
B[i,]<-y
######################################
E<-matrix(0,nrow(B),length(x)-nrow(A))
for( i in 1:nrow(B) )
q=0
for( j in 1:ncol(B) )
if( B[i,j]!=1 )
q=q+1
E[i,q]<-j
########################################
ASD<-combn(E[i,],2)
F<-matrix(0,nrow(B)*ncol(ASD),length(x))
q=0
for( i in 1:nrow(B) )
ASD<-combn(E[i,],2)
for( j in 1:ncol(ASD) )
B[i,ASD[1:nrow(ASD),j]]<-2
q=q+1
F[q,]<-B[i,]
B[i,ASD[1:nrow(ASD),j]]<-0
还有其他不那么复杂的方法吗？
【问题讨论】：
你真的不应该将变量声明为F。这令人困惑，因为F 是逻辑值FALSE 的简写。因此，如果您再次使用 F，则遵循此声明的任何需要逻辑值的函数都可能崩溃。
【参考方案1】：
这是来自包RcppAlgos的开发版本的超快单行。
devtools::install_github("jwood000/RcppAlgos")
library(RcppAlgos)
myPerms <– permuteGeneral(3,6,TRUE,"prod","==",36) - 1L
myPerms
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 1 1 2 2
[2,] 0 0 1 2 1 2
[3,] 0 0 1 2 2 1
[4,] 0 0 2 1 1 2
[5,] 0 0 2 1 2 1
[6,] 0 0 2 2 1 1
.
.
.
[,1] [,2] [,3] [,4] [,5] [,6]
[85,] 2 2 0 0 1 1
[86,] 2 2 0 1 0 1
[87,] 2 2 0 1 1 0
[88,] 2 2 1 0 0 1
[89,] 2 2 1 0 1 0
[90,] 2 2 1 1 0 0
这里有一些基准测试，其中rcppAlgo、r2eOne、r2eTwo 和 OPFun 是每个方法的代码的函数包装器。
microbenchmark(rcppAlgo(),r2eOne(),r2eTwo(),OPFun(N=6) unit = "relative")
Unit: relative
expr min lq mean median uq max neval
rcppAlgo() 1.00000 1.00000 1.00000 1.0000 1.00000 1.000000 100
r2eOne() 471.56007 473.26487 194.01669 267.9402 274.46604 8.373630 100
r2eTwo() 50.71091 48.84173 24.01617 27.8441 34.02326 2.044374 100
OPFun(N=6) 37.35899 24.38966 22.38029 19.7059 19.51935 31.18059 100
说明
由于 OP 正在寻找具有特定频率的数字的特定组合，我们可以使用 Fundamental theorem of arithmetic，它指出每个数字都可以写为质数的唯一组合的乘积。我们得到集合0, 1, 2，加1得到集合1, 2, 3。我们这样做是为了避免在取产品时出现很多零。
现在，我们的任务是找到所有组合，以使每个元素恰好出现两次。这意味着在我们将产品应用于我们的目标组合后，我们得到1*1*2*2*3*3 = 36（注意1 不是质数，但可以忽略，因为1*n = n for all n）。现在问题很简单。
我们只需找到乘积等于36 的所有组合，然后减去1 即可回到我们原来的一组数字，瞧！
一般解决方案
下面，我们有一个通用解决方案，可用于查找给定向量的所有排列，其中每个元素重复特定次数。
library(RcppAlgos) ## for primeSieve and permuteGeneral
MakePerms <- function(v, numReps, myCap = NULL)
m <- sum(numReps)
n <- length(v)
## Generate some primes using prime
## number theorem; fudging a bit to
## ensure we get n-1 prime numbers
myPs <- primeSieve(2*n*log(n))[1:(n-1)]
## Set up vector that will be tested
myV <- c(1L, myPs)
target <- prod(myV^numReps)
ps <- permuteGeneral(myV, m, TRUE, "prod", "==", target, myCap)
for (j in 1:n) ps[ps == myV[j]] <- v[j]
ps
它严重依赖于根据算术基本定理的素数分解的唯一性和一点索引（不像上面的简单情况那么简单，但仍然只有 7 行并且仍然非常快）。
我们首先创建第一个n-1 素数的向量，然后添加1 以完成myV。然后，我们将myV 的每个元素提升到numReps 给出的每个元素所需的重复次数，并获取我们的target 值。以下是一些示例：
v = c(10,13,267,1) 和 numReps = c(3,1,2,5) -->> myV = c(1,2,3,5) -->> target = 1^3 * 2^1 * 3^2 * 5^5 = 56250
v = 0:5 和 numReps = c(1,2,1,2,2,2) -->> myV = c(1,2,3,5,7,11) -->> target = 1^1 * 2^2 * 3^1 * 5^2 * 7^2 * 11^2 = 1778700
OP 示例：v = c(0,1,2) 和 numReps = c(2,2,2) -->> myV = c(1,2,3) -->> target = 1^2 * 2^2 * 3^2 = 36
在我们找到乘积等于target 值的所有排列后，我们只需使用索引将原始向量v 的内容映射到生成的矩阵。
例如，如果您在 OP 的示例中设置 N = 8，您会得到 c(0,1,2) 的所有排列，其中 0 完全重复 4 次，并且 1 和 2 重复两次。
t1 <- OPFun(N=8)
t2 <- MakePerms(0:2, c(4,2,2))
all.equal(t1[do.call(order, as.data.frame(t1)), ],
t2[do.call(order, as.data.frame(t2)), ])
[1] TRUE
microbenchmark(fun2(8), MakePerms(0:2, c(4,2,2)), unit = "relative")
Unit: relative
expr min lq mean median uq max neval
OPFun(8) 23.25099 22.56178 18.64762 19.52436 18.37387 10.90934 100
MakePerms(0:2, c(4, 2, 2)) 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 100
需要注意的是，可能的排列数量增长很快，所以像MakePerms(0:5, rep(2, 6)) 这样的尝试会失败，因为0:5 12 times 的排列总数是12^6 = 2,985,984 > 2^31 - 1（即矩阵的最大行数在Rcpp)。但是，我们并不期望所有这些都符合我们的标准，所以如果我们设置一个上限，比如10^7，我们就会成功。观察：
a <- MakePerms(0:5, rep(2, 6), 10^7)
nrow(a)
7484400
set.seed(17)
a[sample(nrow(a), 10), ]
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,] 0 5 3 3 1 2 4 4 5 1 0 2
[2,] 5 4 2 1 1 0 3 4 5 2 3 0
[3,] 2 4 5 3 5 1 3 0 1 0 4 2
[4,] 4 3 3 1 2 5 0 5 4 1 0 2
[5,] 2 2 5 3 4 1 0 3 5 1 0 4
[6,] 3 1 1 5 0 3 2 0 2 4 4 5
[7,] 1 1 4 2 0 5 4 0 3 5 3 2
[8,] 1 0 4 2 4 2 5 1 3 0 5 3
[9,] 4 3 4 1 5 0 0 2 2 1 3 5
[10,] 1 0 5 3 2 0 1 4 3 4 2 5
使用myCap也可以大大提高效率。
microbenchmark(withOutCap = MakePerms(0:5, c(1,2,1,2,1,2)),
withCap = MakePerms(0:5, c(1,2,1,2,1,2), 10^5),
times = 15)
Unit: milliseconds
expr min lq mean median uq max neval
withOutCap 219.64847 246.4718 275.04672 282.52829 299.33816 311.2031 15
withCap 22.56437 30.6904 33.30469 31.70443 37.50858 41.6095 15
identical(MakePerms(0:5, c(1,2,1,2,1,2)), MakePerms(0:5, c(1,2,1,2,1,2), 10^5))
[1] TRUE
iterpc解决方案
似乎在这一点上提供的答案是严格的学术性的，因为@StéphaneLaurent 提供的答案要好得多。超级通用，一条线，超级快！！
microbenchmark(iter = getall(iterpc(c(2,2,2), labels=c(0,1,2), ordered=TRUE)),
rcppAlg = MakePerms(0:2, c(2,2,2)))
Unit: microseconds
expr min lq mean median uq max neval
iter 428.885 453.2975 592.53164 540.154 683.9585 1165.772 100
rcppAlg 62.418 74.5205 93.44926 81.749 108.4660 216.454 100
故事随着排列数量的增加而变化。观察：
microbenchmark(iter = getall(iterpc(c(2,2,2,2), labels=c(0,1,2,3), ordered=TRUE)),
rcppAlg = MakePerms(0:3, c(2,2,2,2)),
rcppAlgCap = MakePerms(0:3, c(2,2,2,2), 5000))
Unit: microseconds
expr min lq mean median uq max neval
iter 877.246 1052.7060 1394.636 1150.0895 1265.088 8914.980 100
rcppAlg 964.446 1449.7115 2084.944 1787.9350 1906.242 10921.156 100
如果你使用myCap，MakePerms 会快一点。这并不重要，因为使用iterpc 解决方案，您甚至不必考虑将获得多少结果。很不错！！
更新
RcppAlgos（我是作者）的新版本刚刚在 CRAN 上发布。 permuteGeneral 现在有一个额外的参数，称为 freqs，它允许多重集的排列，这正是 OP 正在寻找的。p>
microbenchmark(iter = getall(iterpc(c(2,2,2,2), labels=0:3, ordered=TRUE)),
newRcppAlgos = permuteGeneral(0:3, freqs = c(2,2,2,2)))
Unit: microseconds
expr min lq mean median uq max neval
iter 457.442 482.8365 609.98678 508.6150 572.581 4037.048 100
newRcppAlgos 33.159 43.3975 56.40026 48.5665 58.194 625.691 100
microbenchmark(iter = getall(iterpc(c(5,4,3,2), labels=0:3, ordered=TRUE)),
newRcppAlgos = permuteGeneral(0:3, freqs = c(5,4,3,2)))
Unit: milliseconds
expr min lq mean median uq max neval
iter 480.25976 552.54343 567.9155 565.23066 579.0258 751.8556 100
newRcppAlgos 83.41194 87.03957 104.6279 95.67596 107.3572 181.1119 100
identical(getall(iterpc(c(5,4,3,2), labels=0:3, ordered=TRUE)),
permuteGeneral(0:3, freqs = c(5,4,3,2)))
[1] TRUE
nrow(permuteGeneral(0:3, freqs = c(5,4,3,2)))
[1] 2522520
更新 2
正如@StéphaneLaurent 所指出的，包arrangements 已作为iterpc 的替代品发布（参见@RandyLai 的cmets）。它效率更高，并且能够处理更广泛的组合问题（例如分区）。以下是较大示例的基准：
microbenchmark(arrangements = permutations(x = 0:3, freq = c(5,4,3,2)),
RcppAlgos = permuteGeneral(0:3, freqs = c(5,4,3,2)))
Unit: milliseconds
expr min lq mean median uq max neval
arrangements 97.10078 98.67154 113.5953 100.56261 131.3244 163.8912 100
RcppAlgos 92.13122 93.84818 108.1845 95.72691 101.2647 165.7248 100
...几乎相同的结果。
arrangements 的一个巨大好处是能够通过getnext 一次（或分块）获得一个排列。这允许用户生成超过2^31 - 1 的结果，并提供更大的灵活性。
有关R 中此类问题的更多信息，我在问题中写了extensive overview：R: Permutations and combinations with/without replacement and for distinct/non-distinct items/multiset。
【讨论】：
它有效！伟大的！如何判断我是否有N=8 元素q=4 元素，即0,1,2,3？同样2元素等于0，2元素等于1，2元素等于2和2元素等于3在所有可能的位置？所以我说的是N=2q时的通用解决方案。
我很高兴这有帮助。今天我有时会发布一个通用的解决方案。下班了！！！
@VassilisChasiotis，我已经更新了一般情况的答案。享受吧！
有关信息，iterpc 已被弃用，取而代之的是 arrangements。这里有一些基准：randy3k.github.io/arrangements/articles/benchmark.html
@StéphaneLaurent 是的，我是aware。自发布以来，我一直在我的答案中使用 arrangements 包。鉴于此，我将更新我的答案。【参考方案2】：
您可以使用iterpc 包：
> library(iterpc)
> I <- iterpc(c(2,2,2), labels=c(0,1,2), ordered=TRUE)
> getall(I)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 0 0 1 1 2 2
[2,] 0 0 1 2 1 2
[3,] 0 0 1 2 2 1
[4,] 0 0 2 1 1 2
[5,] 0 0 2 1 2 1
[6,] 0 0 2 2 1 1
[7,] 0 1 0 1 2 2
[8,] 0 1 0 2 1 2
[9,] 0 1 0 2 2 1
[10,] 0 1 1 0 2 2
[11,] 0 1 1 2 0 2
[12,] 0 1 1 2 2 0
[13,] 0 1 2 0 1 2
[14,] 0 1 2 0 2 1
[15,] 0 1 2 1 0 2
[16,] 0 1 2 1 2 0
[17,] 0 1 2 2 0 1
[18,] 0 1 2 2 1 0
[19,] 0 2 0 1 1 2
[20,] 0 2 0 1 2 1
[21,] 0 2 0 2 1 1
[22,] 0 2 1 0 1 2
[23,] 0 2 1 0 2 1
[24,] 0 2 1 1 0 2
[25,] 0 2 1 1 2 0
[26,] 0 2 1 2 0 1
[27,] 0 2 1 2 1 0
[28,] 0 2 2 0 1 1
[29,] 0 2 2 1 0 1
[30,] 0 2 2 1 1 0
[31,] 1 0 0 1 2 2
[32,] 1 0 0 2 1 2
[33,] 1 0 0 2 2 1
[34,] 1 0 1 0 2 2
[35,] 1 0 1 2 0 2
[36,] 1 0 1 2 2 0
[37,] 1 0 2 0 1 2
[38,] 1 0 2 0 2 1
[39,] 1 0 2 1 0 2
[40,] 1 0 2 1 2 0
[41,] 1 0 2 2 0 1
[42,] 1 0 2 2 1 0
[43,] 1 1 0 0 2 2
[44,] 1 1 0 2 0 2
[45,] 1 1 0 2 2 0
[46,] 1 1 2 0 0 2
[47,] 1 1 2 0 2 0
[48,] 1 1 2 2 0 0
[49,] 1 2 0 0 1 2
[50,] 1 2 0 0 2 1
[51,] 1 2 0 1 0 2
[52,] 1 2 0 1 2 0
[53,] 1 2 0 2 0 1
[54,] 1 2 0 2 1 0
[55,] 1 2 1 0 0 2
[56,] 1 2 1 0 2 0
[57,] 1 2 1 2 0 0
[58,] 1 2 2 0 0 1
[59,] 1 2 2 0 1 0
[60,] 1 2 2 1 0 0
[61,] 2 0 0 1 1 2
[62,] 2 0 0 1 2 1
[63,] 2 0 0 2 1 1
[64,] 2 0 1 0 1 2
[65,] 2 0 1 0 2 1
[66,] 2 0 1 1 0 2
[67,] 2 0 1 1 2 0
[68,] 2 0 1 2 0 1
[69,] 2 0 1 2 1 0
[70,] 2 0 2 0 1 1
[71,] 2 0 2 1 0 1
[72,] 2 0 2 1 1 0
[73,] 2 1 0 0 1 2
[74,] 2 1 0 0 2 1
[75,] 2 1 0 1 0 2
[76,] 2 1 0 1 2 0
[77,] 2 1 0 2 0 1
[78,] 2 1 0 2 1 0
[79,] 2 1 1 0 0 2
[80,] 2 1 1 0 2 0
[81,] 2 1 1 2 0 0
[82,] 2 1 2 0 0 1
[83,] 2 1 2 0 1 0
[84,] 2 1 2 1 0 0
[85,] 2 2 0 0 1 1
[86,] 2 2 0 1 0 1
[87,] 2 2 0 1 1 0
[88,] 2 2 1 0 0 1
[89,] 2 2 1 0 1 0
[90,] 2 2 1 1 0 0
编辑 2018-04-28
iterpc 现在已弃用，取而代之的是 arrangements。
【讨论】：
非常好！！iterpc 看起来很棒。我将在下面更新我的基准。
我认为iterpc 不会被弃用。在 R 中弃用意味着它将被删除并破坏任何依赖它的代码。他的确切措辞是“arrangements 被认为是iterpc 的替代品”。我会联系 Randy Lai 以澄清问题。【参考方案3】：
建议了两种方法：一种效率低下，一种更有效但费力。 （在这种情况下，我将“效率”等同于缩放，而不是代码量或执行时间。也就是说，只要你只创建 90 行，那么你就可以了。如果这是一个简化的问题并且您确实需要扩展到更大的矩阵，那么permutations 可能会超出内存和/或 R 的容量。）
这两种解决方案都比您的代码短一些。第一个比较清晰，只有4行代码；第二个确实有点晦涩（并且似乎进入了索引间接“inception”），但实际上仍然只需要 13 行代码。第二个可能会有所减少，但我的“游戏”时间用完了:-)
效率低下
一种方法是创建所有排列并过滤掉重复。只要您的“N”不会太大，这就会起作用。
library(gtools)
v <- rep(0:2, 2)
p <- permutations(6, 6)
p[] <- v[p]
p <- p[!duplicated(p),]
head(p)
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 0 1 2 0 1 2
# [2,] 0 1 2 0 2 1
# [3,] 0 1 2 1 0 2
# [4,] 0 1 2 1 2 0
# [5,] 0 1 2 2 0 1
# [6,] 0 1 2 2 1 0
tail(p)
# [,1] [,2] [,3] [,4] [,5] [,6]
# [85,] 2 2 0 1 0 1
# [86,] 2 2 0 1 1 0
# [87,] 2 2 0 0 1 1
# [88,] 2 2 1 0 0 1
# [89,] 2 2 1 0 1 0
# [90,] 2 2 1 1 0 0
验证每一行的每个元素都恰好有两个：
all(apply(p, 1, table) == 2)
# [1] TRUE
低效
一种效率较低（因此更加努力）的方法：使用combn(6,2) 和combn(4,2) 创建列索引矩阵，然后适当地分配“因子”。 （稍后会更有意义。）
（注意：我经常在转置矩阵方面更好地考虑这些问题；您可以轻松地做到这一点，只需调整代码以交换列/行。）
我们需要的是类似expand.grid 的东西，一次用于两列。所以我们将从较小的问题开始：
left2 <- t(combn(6, 2))
mid2 <- t(combn(4, 2))
left2
# [,1] [,2]
# [1,] 1 2
# [2,] 1 3
# [3,] 1 4
# [4,] 1 5
# [5,] 1 6
# [6,] 2 3
# [7,] 2 4
# [8,] 2 5
# [9,] 2 6
# [10,] 3 4
# [11,] 3 5
# [12,] 3 6
# [13,] 4 5
# [14,] 4 6
# [15,] 5 6
mid2
# [,1] [,2]
# [1,] 1 2
# [2,] 1 3
# [3,] 1 4
# [4,] 2 3
# [5,] 2 4
# [6,] 3 4
现在，网格将在这两个矩阵的行索引上展开。
eg <- expand.grid(a = 1:15, b = 1:6)
head(eg)
# a b
# 1 1 1
# 2 2 1
# 3 3 1
# 4 4 1
# 5 5 1
# 6 6 1
inds <- cbind(left2[eg$a,], mid2[eg$b,])
head(inds)
# [,1] [,2] [,3] [,4]
# [1,] 1 2 1 2
# [2,] 1 3 1 2
# [3,] 1 4 1 2
# [4,] 1 5 1 2
# [5,] 1 6 1 2
# [6,] 2 3 1 2
inds[25,,drop=FALSE]
# [,1] [,2] [,3] [,4]
# [1,] 3 4 1 3
这意味着，对于第 25 行，我们应该将第 3 列和第 4 列替换为第一个因子（例如 0）。然后，其余列中 (1,2,5,6)，我们应该将第 1 列和第 3 列替换为第二个因子（例如 1）。再说一遍，c(1,2,5,6)[c(1,3)] 等于将第 1 列和第 5 列替换为第二个值 (1)。 （第三个值2 将进入所有剩余的插槽。）
所以，要从上面找到c(1,2,5,6)，我们可以使用setdiff(1:6,...)：
afterleft2 <- t(apply(left2[eg$a,], 1, function(a) setdiff(1:6, a)))
head( afterleft2 )
# [,1] [,2] [,3] [,4]
# [1,] 3 4 5 6
# [2,] 2 4 5 6
# [3,] 2 3 5 6
# [4,] 2 3 4 6
# [5,] 2 3 4 5
# [6,] 1 4 5 6
afterleft2[25,,drop=FALSE]
# [,1] [,2] [,3] [,4]
# [1,] 1 2 5 6
让我们开始修复inds 第三和第四列。
inds[,3] <- afterleft2[ cbind(1:90, mid2[eg$b,1]) ]
inds[,4] <- afterleft2[ cbind(1:90, mid2[eg$b,2]) ]
head(inds)
# [,1] [,2] [,3] [,4]
# [1,] 1 2 3 4
# [2,] 1 3 2 4
# [3,] 1 4 2 3
# [4,] 1 5 2 3
# [5,] 1 6 2 3
# [6,] 2 3 1 4
inds[25,,drop=FALSE]
# [,1] [,2] [,3] [,4]
# [1,] 3 4 1 5
由此我们看到第 25 行有我们期望的“1”和“5”。
现在结束：
nr <- nrow(inds)
out <- matrix(nrow = nr, ncol = 6L)
out[cbind(1:nr,inds[,1])] <- 0L
out[cbind(1:nr,inds[,2])] <- 0L
out[cbind(1:nr,inds[,3])] <- 1L
out[cbind(1:nr,inds[,4])] <- 1L
head(out)
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 0 0 1 1 NA NA
# [2,] 0 1 0 1 NA NA
# [3,] 0 1 1 0 NA NA
# [4,] 0 1 1 NA 0 NA
# [5,] 0 1 1 NA NA 0
# [6,] 1 0 0 1 NA NA
out[25,,drop=FALSE]
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 1 NA 0 0 1 NA
我上面提到的“剩余插槽”（对于第三个值）都是NA，按设计。
out[is.na(out)] <- 2L
head(out)
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 0 0 1 1 2 2
# [2,] 0 1 0 1 2 2
# [3,] 0 1 1 0 2 2
# [4,] 0 1 1 2 0 2
# [5,] 0 1 1 2 2 0
# [6,] 1 0 0 1 2 2
out[25,,drop=FALSE]
# [,1] [,2] [,3] [,4] [,5] [,6]
# [1,] 1 2 0 0 1 2
现在进行快速健全性检查，以确保我们的 out 变量在每行中恰好有两个元素。
all(apply(out, 1, table) == 2)
# [1] TRUE
【讨论】：
在这样开头的四行代码中：out[cbind(1:nr,inds[,1])]...，nr 是未定义的。我只是添加了一个变量nr <- nrow(inds)，因为我假设您已经在您的环境中拥有。任何人，我真的很喜欢你的两种方法以及你的指导性解释。
是的，就是这样，不确定如何复制/粘贴，谢谢以上是关于R中N个元素与q个元素的组合的主要内容，如果未能解决你的问题，请参考以下文章