Efficiently create random sample from expand.grid output without using expand.grid-CodePudding

Let's suppose I have 3 vectors (just as an example). Now I want to obtain a random sample over all possible combinations of those 3. Usually, I'd do this:

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, ]

That's fine for small vectors. But if those vectors grow in size, my_grid will get very big very fast. So the question is, how to create the result by only using idx and the three vectors?

CodePudding user response：

To avoid the expand.grid, you can use a Cantor expansion of the integers. This is a generalization of the binary expansion. See the above link for details. In short, each integer in 1:(n1*n2*n3) has a Cantor expansion (x1, x2, x3) with x1 in 1:n1, x2 in 1:n2, x3 in 1:n3. The binary expansion is the case n1 = n2 = n3 = 2.

Here is the code for your example:

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]])
}

The above link provides a Rcpp function to replace intToCantor.

CodePudding user response：

This should work:

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])]
)

Basically, the columns of the output of expand.grid for each vector are formed by blocks of values and the length of each block is a product of the lengths of the "preceding" vectors. So you just divide the index by this product and than take modulo from this number to find which value would be on that position.