0

想知道是否有更快的方法来实现以下目标：

我有一个最大向量。长度long1全为零。如何在某些位置生成所有可能组合的矩阵，直到并包括最大值max1s。

以下方法有效，但当返回矩阵很大时（例如 >1e5）似乎效率很低。

### longest series of zeros
long1 <- 4
### max no. of 1s
max1s <- 2
### get combinations up to long1
f1 <- function(i) t(combinat::combn(seq.int(long1),
                                    m=i, simplify=TRUE))
### list of positions in target matrix to be made 1s
### each list element represents column positions of 1s
l1 <- sapply(1:max1s, f1)
### no. rows in return matrix
nrow1 <- sum(unlist(lapply(l1, nrow)))
### set up matrix of zeros
c2 <- matrix(0L, nrow=nrow1, ncol=long1)
### rows to start at for each 'i' in length(l1) below
nrow2 <- c(1, 1+cumsum(unlist(lapply(l1, nrow))))
for (i in 1:length(l1)){
    for (j in 1:nrow(l1[[i]])){
### now iterate over each row in that element of l1
### set relevant position in matrix to 1
        c2[nrow2[i]+(j-1), l1[[i]][j, ] ] <-  1L
    }}

在这种情况下，它1是长度为 4 的向量中 s 的所有组合，最大为 2：

> c2
      [,1] [,2] [,3] [,4]
 [1,]    1    0    0    0
 [2,]    0    1    0    0
 [3,]    0    0    1    0
 [4,]    0    0    0    1
 [5,]    1    1    0    0
 [6,]    1    0    1    0
 [7,]    1    0    0    1
 [8,]    0    1    1    0
 [9,]    0    1    0    1
[10,]    0    0    1    1

我宁愿避免使用combinat::hcubethen 消除超过一定数量的行。s ，1因为这种方法会为这样的应用程序创建不必要的大矩阵。

4

1 回答 1

1

我想您可以使用单独计算每种尺寸的组合，combn然后使用do.callwithrbind将它们组合在一起：

allcombo <- function(long1, max1s) {
  do.call(rbind, lapply(1:max1s, function(num1) {
    t(apply(combn(long1, num1), 2, function(x) {
      col = rep(0, long1)
      col[x] = 1
      col
    }))
  }))
}

我已将您发布的解决方案存储在 function 中OP。我们可以检查它们是否返回相同的值：

all.equal(OP(20, 5), allcombo(20, 5))
# [1] TRUE

现在我们可以进行基准测试（返回 21699 行）：

library(microbenchmark)
microbenchmark(OP(20, 5), allcombo(20, 5))
# Unit: milliseconds
#             expr      min       lq   median       uq      max neval
#        OP(20, 5) 242.4120 256.5791 269.7237 292.7131 556.5984   100
#  allcombo(20, 5) 150.4291 179.2588 188.4840 200.9898 448.2214   100

所以这种方法的使用combn速度要快一些（这组参数在我的电脑上是 30%）。

于 2014-05-06T18:15:19.487 回答