我有一个由1和0组成的矩阵，行是个人，列是事件，1表示事件发生在个人身上，0表示没有。
我想找出哪一组（在示例中）5列/事件覆盖了最多的行/个人。

测试数据

#Make test data
set.seed(123)
d <- sapply(1:300, function(x) sample(c(0,1), 30, T, c(0.9,0.1)))
colnames(d) <- 1:300
rownames(d) <- 1:30

字符串

我的努力

我最初的尝试只是基于将5列的集合与最高的colMeans相结合：

#Get top 5 columns with highest row coverage
col_set <- head(sort(colMeans(d), decreasing = T), 5)

#Have a look the set
col_set

>
      197       199        59        80        76 
0.2666667 0.2666667 0.2333333 0.2333333 0.2000000

#Check row coverage of the column set
sum(apply(d[,colnames(d) %in% names(col_set)], 1, sum) > 0) / 30 #top 5

>
[1] 0.7

然而，这个集合并没有覆盖大多数行。我通过伪随机抽样10.000个不同的5列集合来测试，然后找到覆盖率最高的集合：

#Get 5 random columns using colMeans as prob in sample
##Random sample 10.000 times
set.seed(123)
result <- lapply(1:10000, function(x){
  col_set2 <- sample(colMeans(d), 5, F, colMeans(d))
  cover <- sum(apply(d[,colnames(d) %in% names(col_set2)], 1, sum) > 0) / 30 #random 5
  list(set = col_set2, cover = cover)
})

##Have a look at the best set
result[which.max(sapply(result, function(x) x[["cover"]]))]

>
[[1]]
[[1]]$set
        59        169        262         68        197 
0.23333333 0.10000000 0.06666667 0.16666667 0.26666667 

[[1]]$cover
[1] 0.7666667

型
提供colMeans到sample的原因是具有最高覆盖率的列是我最感兴趣的列。
因此，使用伪随机抽样，我可以收集一组覆盖率比只使用前5列更高的列。然而，由于我的实际数据集比示例更大，我正在寻找一种更有效和合理的方法来找到具有最高覆盖率的列集。

编辑

对于感兴趣的，我决定microbenchmark提供的3个解决方案：

#Defining G. Grothendieck's coverage funciton outside his solutions
coverage <- function(ix) sum(rowSums(d[, ix]) > 0) / 30

#G. Grothendieck top solution
solution1 <- function(d){
  cols <- tail(as.numeric(names(sort(colSums(d)))), 20)
  co <- combn(cols, 5)
  itop <- which.max(apply(co, 2, coverage))
  co[, itop]
}

#G. Grothendieck "Older solution"
solution2 <- function(d){
  require(lpSolve)
  ones <- rep(1, 300)
  res <- lp("max", colSums(d), t(ones), "<=", 5, all.bin = TRUE, num.bin.solns = 10)
  m <- matrix(res$solution[1:3000] == 1, 300)
  cols <- which(rowSums(m) > 0)
  co <- combn(cols, 5)
  itop <- which.max(apply(co, 2, coverage))
  co[, itop]
}

#user2554330 solution
bestCols <- function(d, n = 5) {
  result <- numeric(n)
  for (i in seq_len(n)) {
    result[i] <- which.max(colMeans(d))
    d <- d[d[,result[i]] != 1,, drop = FALSE]
  }
  result
}

#Benchmarking...
microbenchmark::microbenchmark(solution1 = solution1(d),
                               solution2 = solution2(d),
                               solution3 = bestCols(d), times = 10)

>
Unit: microseconds
      expr        min         lq       mean      median         uq       max neval
 solution1 390811.850 497155.887 549314.385 578686.3475 607291.286 651093.16    10
 solution2  55252.890  71492.781  84613.301  84811.7210  93916.544 117451.35    10
 solution3    425.922    517.843   3087.758    589.3145    641.551  25742.11    10

型