選擇組合的子集

Question

假設我有一個20×5的矩陣，我想選擇矩陣的子集並且用它們做一些計算。進一步假設每個子矩陣是7×5。我當然可以做

ncomb <- combn(20, 7) 

這給我的7個索引所有可能的組合，我可以用這些來獲取子矩陣。但是用一個小的20×5矩陣，已經有77520個可能的組合。所以我想隨機抽樣一些組合，例如5000個組合。

一種可能性是以下內容：

ncomb <- combn(20, 7) 
ncombsub <- ncomb[, sample(77520, 5000)] 

換句話說，我得到所有可能的組合，然後隨機選擇只5000組合。但我想如果我有一個更大的矩陣計算所有可能的組合，將會有問題 - 比方說100 X 7.

所以我想知道是否有辦法在沒有首先獲得所有可能的組合的情況下獲得組合的子集。

Answer 1

最後我做什麼@Roland建議，通過修改combn()和字節編譯代碼：

combn_sub <- function (x, m, nset = 5000, seed=123, simplify = TRUE, ...) { 
    stopifnot(length(m) == 1L) 
    if (m < 0) 
     stop("m < 0", domain = NA) 
    if (is.numeric(x) && length(x) == 1L && x > 0 && trunc(x) == 
     x) 
     x <- seq_len(x) 
    n <- length(x) 
    if (n < m) 
     stop("n < m", domain = NA) 
    m <- as.integer(m) 
    e <- 0 
    h <- m 
    a <- seq_len(m) 
    len.r <- length(r <- x[a]) 
    count <- as.integer(round(choose(n, m))) 
    if(count < nset) nset <- count 
    dim.use <- c(m, nset)  

    ##-----MOD 1: Change the output matrix size-------------- 
    out <- matrix(r, nrow = len.r, ncol = nset) 

    if (m > 0) { 
     i <- 2L 
     nmmp1 <- n - m + 1L 

     ##----MOD 2: Select a subset of indices 
     set.seed(seed) 
     samp <- sort(c(1, sample(2:count, nset - 1))) 

     ##----MOD 3: Start a counter. 
     counter <- 2L  

     while (a[1L] != nmmp1) { 
      if (e < n - h) { 
       h <- 1L 
       e <- a[m] 
       j <- 1L 
      } 
      else { 
       e <- a[m - h] 
       h <- h + 1L 
       j <- 1L:h 
      } 
      a[m - h + j] <- e + j 

      #-----MOD 4: Whenever the counter matches an index in samp, 
      #a combination of row indices is produced and stored in the matrix `out` 
      if(samp[i] == counter){ 
       out[, i] <- x[a] 
       if(i == nset) break 
       i <- i + 1L 
      } 
      #-----Increase the counter by 1 for each iteration of the while-loop 
      counter <- counter + 1L 
     } 
    } 
    array(out, dim.use) 
} 

library("compiler") 
comb_sub <- cmpfun(comb_sub) 

Answer 2

你的方法：

op <- function(){ 
    ncomb <- combn(20, 7) 
    ncombsub <- ncomb[, sample(choose(20,7), 5000)] 
    return(ncombsub) 
} 

了不同的策略，簡單的樣品從原來的矩陣七行5000次（直到5000個唯一的行組合找到了一個新的樣本來替換任何重複的樣品）：

me <- function(){ 
    rowsample <- replicate(5000,sort(sample(1:20,7,FALSE)),simplify=FALSE) 
    while(length(unique(rowsample))<5000){ 
    rowsample <- unique(rowsample) 
    rowsample <- c(rowsample, 
        replicate(5000-length(rowsample), 
           sort(sample(1:20,7,FALSE)),simplify=FALSE)) 
    } 
    return(do.call(cbind,rowsample)) 
} 

這應該是更有效率，因爲它可以防止您必須首先計算所有組合，隨着矩陣變大，這會變得昂貴。

然而，一些基準測試顯示情況並非如此。至少在這個矩陣：

library(microbenchmark) 
microbenchmark(op(),me()) 

Unit: milliseconds 
expr  min  lq median  uq  max neval 
op() 184.5998 201.9861 206.3408 241.430 299.9245 100 
me() 411.7213 422.9740 429.4767 474.047 490.3177 100