无需使用expand.grid即可从expand.grid输出高效创建随机样本

Question

假设我有 3 个向量（仅作为示例）。现在我想获得这 3 种所有可能组合的随机样本。通常，我会这样做：

x <- 1:3
y <- 10:12
z <- 15:18

N <- length(x) * length(y) * length(z) # Length of the resulting grid
idx <- sample(1:N, 10, replace = FALSE)
my_grid <- expand.grid(x, y, z)
result <- my_grid[idx, ]

这对于小向量来说没问题。但如果这些向量的大小变大，

my_grid

就会变得非常快。那么问题是，如何仅使用

result

和三个向量来创建

idx

？

Answer 1

这应该有效：

X <- list(x, y, z)
X.len <- sapply(X, length)

# modify '%%' to return values 1..n instead of 0..(n-1)
mod2 <- function(x, y) (x-1) %% y + 1

result <- sapply(seq_along(X), function(i) 
  X[[i]][mod2(ceiling(idx / prod(X.len[seq_len(i-1)])),  
              X.len[i])]
)

基本上，每个向量的

expand.grid

输出的列由值块构成，每个块的长度是“前面”向量的长度的乘积。因此，您只需将索引除以该乘积，然后对该数字取模即可找到该位置上的值。

Answer 2

为了避免

expand.grid

，您可以使用整数的 Cantor 展开。这是二进制展开式的推广。详情请参阅上面的链接。简而言之，

1:(n1*n2*n3)

中的每个整数都有一个康托展开式

(x1, x2, x3)

，其中

x1

中的

1:n1

，

x2

中的

1:n2

，

x3

中的

1:n3

。二进制展开就是这种情况

n1 = n2 = n3 = 2

.

这是您的示例的代码：

intToCantor <- function(n, sizes) {
  l <- c(1, cumprod(sizes))
  epsilon <- numeric(length(sizes))
  while(n>0){
    k <- which.min(l<=n)
    e <- floor(n / l[k-1])
    epsilon[k-1] <- e
    n <- n - e*l[k-1]
  }
  epsilon
}

CantorToInt <- function(epsilon, sizes) {
  sum(epsilon * c(1, cumprod(sizes[1:(length(epsilon)-1)])))
}

x <- 1:3
y <- 10:12
z <- 15:18

sizes <- c(length(x), length(y), length(z))
N <- prod(sizes)
n <- 10
idx <- sample(1:N, n, replace = FALSE)

result <- matrix(NA_real_, nrow = n, ncol = length(sizes))
for(i in 1:n) {
  indices <- 1 + intToCantor(idx[i] - 1, sizes = sizes)
  result[i, ] <- c(x[indices[1]], y[indices[2]], z[indices[3]])
}

上面的链接提供了一个Rcpp函数来替代

intToCantor

。

---

编辑

@Robert Hacken 的回答提供了一种更有效的方法。事实上，他的答案使用了隐藏的康托展开式，比我的更快：

mod2 <- function(x, y) (x-1) %% y + 1

CantorExpansion <- function(n, sizes) {
  p <- cumprod(c(1, head(sizes, -1)))
  mod2(ceiling(n / p), sizes)
}

---

编辑

Rcpp功能更快：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
IntegerVector CantorRcpp(int n, std::vector<int> sizes) {
  IntegerVector epsilon(sizes.size(), 1);
  std::vector<int>::iterator it = sizes.begin();
  it = sizes.insert(it, 1);
  int G[sizes.size()];
  std::partial_sum(sizes.begin(), sizes.end(), G, std::multiplies<int>());
  n--;
  int k;
  while(n > 0) {
    k = 1;
    while(G[k] <= n) {
      k += 1;
    }
    int d = G[k-1];
    epsilon(k-1) = 1 + n / d;
    n = n % d;
  }
  return epsilon;
}

/*** R
library(microbenchmark)

CantorExpansion <- function(n, sizes) {
  p <- cumprod(c(1L, head(sizes, -1L)))
  1L + ((ceiling(n / p) - 1L) %% sizes)
}

sizes <- 2L:9L
Robert <- function() {
  L <- vector("list", length = prod(sizes))
  for(n in seq_len(prod(sizes))) {
    L[[n]] <- CantorExpansion(n, sizes)
  }
}
Rcpp <- function() {
  L <- vector("list", length = prod(sizes))
  for(n in seq_len(prod(sizes))) {
    L[[n]] <- CantorRcpp(n, sizes)
  }
}
microbenchmark(
  Robert = Robert(),
  Rcpp = Rcpp(),
  times = 10L
)
# Unit: milliseconds
#    expr       min        lq      mean    median        uq      max neval cld
#  Robert 1658.3666 1690.2473 1743.9026 1728.0874 1765.9536 1910.122    10   b
#    Rcpp  693.8287  764.4371  841.9733  801.9504  947.8049 1050.848    10  a 
*/

一个人可以改进这两个功能，因为累积乘积只能计算一次。

Answer 3

如果您想要一个随机子集而不生成所有组合的完整表，即

expand.grid

，您可以尝试，例如

lst <- list(x, y, z)
n <- 10
res <- list()
repeat {
    if (length(res) == n) break
    v <- list(sapply(lst, sample, 1))
    if (!v %in% res) {
        res <- append(res, v)
    }
}
as.data.frame(do.call(rbind, res))

你会得到类似的东西

   V1 V2 V3
1   2 12 15
2   2 11 18
3   3 10 18
4   3 11 18
5   3 12 16
6   2 11 17
7   2 12 18
8   1 12 15
9   1 12 17
10  2 10 16