在R中生成随机整数对而不替换答案

【问题标题】：Generating Random Pairs of Integers without Replacement in R在R中生成随机整数对而不替换
【发布时间】：2015-04-17 01:13:22
【问题描述】：

我想在不替换的情况下绘制随机整数对（换句话说，我不想要任何重复的对）。这个概念听起来很简单，但我想不出一个快速简单的解决方案。

想象一下，例如，我想使用整数序列1:4 生成随机整数对来填充整数对的元素。还假设我想生成 5 个随机对而不进行替换。然后我希望能够生成这样的东西......

     [,1] [,2]
[1,]    1    2
[2,]    2    1
[3,]    3    3
[4,]    1    4
[5,]    4    3

在上面的示例中，没有重复的对（即行）。然而，在上述矩阵的每一列中都有重复的整数。因此，使用sample() 分别为每一列生成随机数是行不通的。

另一个看似潜在但不适用于我的上下文的解决方案是生成大量包含重复项的对，然后追溯删除这些重复项。我不能这样做，因为我需要生成特定数量的对。

我正在寻找解决此问题的有效方法。这似乎是一个简单的问题，它必须有一个简单的解决方案（即请不要嵌套 for 循环）

这是我丑陋的做法：

#This matrix maps a unique id i.e. (1:16) to a pair (i.e. the row & col of the matrix)
r.mat<-matrix(1:(4*4),4,4) 
#Drawing a random id
r.id<-sample(r.mat,5,replace=FALSE)
#Mapping the random id to a random pair
r.pair<-t(sapply(r.id, function (x) which(r.mat==x,arr.ind=TRUE)))

这适用于我的玩具示例，但是当我想从序列 1:10000000 中绘制大量对时，它就不是很好了。

【问题讨论】：

你如何得到{3,3} 没有替换
您打算从中提取每个数字的序列到底有多大？真的是1e7吗？
rawr - 我基本上从中绘制的集合是 (1,1),(2,1),(1,2),(1,3),(1,4),( 2,2) 等...所以没有替换基本上意味着我永远不会有重复的对。那有意义吗？关于如何改写问题以使其更清晰的任何建议？
BrodieG，我正处于一个项目的开始阶段，所以我不确定序列会有多大，可能不是 1e7。但比 4 更接近 1e7。
你可以通过直接计算行和列来提高最后一行计算r.pair的性能，一个常数时间的操作，代替线性时间的操作which：行是@987654327 @ 并且列为as.integer((x-1)/4)+1L，其中x 与sapply 调用中的函数相同。

标签： r random-sample

【解决方案1】：

这里的关键是不要生成所有排列，因为这在内存和时间方面都非常昂贵。由于您只关心两个数字，因此我们可以非常轻松地做到这一点，只要 (number_of_possible_values) ^ 2 小于双精度浮点中可表示的最大整数：

size <- 1e5
samples <- 100
vals <- sample.int(size ^ 2, samples)
cbind(vals %/% size + 1, vals %% size)

基本上，我们使用整数来表示每个可能的值组合。在我们的示例中，我们从直到1e5 ^ 2 的所有数字中采样，因为我们有1e5 ^ 2 可能的1e5 数字组合。这些1e10 整数中的每一个都代表其中一种组合。然后，我们通过取模作为第一个数，将整数除法作为第二个数，将该整数分解为两个分量值。

基准测试：

Unit: microseconds
                   expr        min         lq       mean
  funBrodie(10000, 100)     16.457     17.188     22.052
 funRichard(10000, 100) 542513.717 640647.919 638045.215

此外，限制应该是 ~3x1e7，并且保持相对较快：

Unit: microseconds
                  expr    min      lq     mean median      uq    max neval
 funBrodie(1e+07, 100) 18.285 20.6625 22.88209 21.211 22.4905 77.893   100

基准函数：

funRichard <- function(size, samples) {
  nums <- 1:size
  dt = CJ(nums, nums)
  dt[sample(1:dim(dt)[1], size = samples), ]
}
funBrodie <- function(size, samples) {
  vals <- sample.int(size ^ 2, samples)
  cbind(vals %/% size + 1, vals %% size)
}

并确认我们正在做类似的事情（请注意，并不是这些应该完全相同，但事实证明它们是相同的）：

set.seed(1)
resB <- funBrodie(1e4, 100)
set.seed(1)
resR <- unname(as.matrix(funRichard(1e4, 100)))
all.equal(resB, resR)
# TRUE

【讨论】：

您是否介意将CJ 方法添加到基准测试中。我很好奇它与其他方法相比如何。
@RichardErickson，查看更新（删除data.table，奇怪的是实际上增加了相当多的开销）。请注意，之前的答案已经在使用 CJ，只是不必要地将其包装在 data.table 调用中（我想这是有道理的，这会复制一个较大的数据集）。
谢谢！这是一种更好的方法。
@JacobH，不客气。另外，请注意，我刚刚意识到这可以完全矢量化以再提高 10-15 倍的速度（请参阅更新以摆脱 sapply）。
假设 size 是 3，samples 是 1。那么，vals 可能是 9，因此您将分配 4, 0，并且 4 和 0 都不在范围 1 中： 3.不应该是((vals-1) %/% size) + 1和(vals %% size) + 1才能回答原来的问题吗？

【解决方案2】：

首先，我找到了如何在SO 上生成对。但是，这并没有扩展，所以我查看了 ?combn 并找到了 expand.grid 函数。

接下来，我使用data.table 包，因为它可以很好地处理大数据（请参阅它的文档了解原因）。

## the data.table library does well with large data sets
library(data.table)

## Small dummy dataset
pairOne = 1:10
pairTwo = 1:2
nSamples = 3

system.time({
dt = data.table(expand.grid(pairOne, pairTwo))
dt2 = dt[sample(1:dim(dt)[1], size = nSamples), ]
})
#   user  system elapsed 
#  0.002   0.001   0.001 

## Large dummy dataset
pairOne = 1:10000
pairTwo = 1:10000
length(pairOne) * length(pairTwo)
nSamples = 1e5
system.time({
dt = data.table(expand.grid(pairOne, pairTwo))
dt2 = dt[sample(1:dim(dt)[1], size = nSamples), ]
})
#   user  system elapsed 
#  2.576   1.276   3.862

【讨论】：

链接答案中的链接问题包含许多有趣的方法和变体。显然是一个具有挑战性的问题，而且你做得很好！
谢谢！这很棒。那个expand.grid函数很方便。
使用CJ() 而不是expand.grid()。
@Arun，CJ() 在什么包中？ ?CJ 没有找到任何东西，谷歌搜索“CJ”R 也没有找到任何有用的结果。
为了将来参考，这里是CJ方法dt = CJ(pairOne, pairTwo); system.time({ dt = CJ(pairOne, pairTwo) dt2 = dt[sample(1:dim(dt)[1], size = nSamples), ] })。这个方法总共0.42秒。哇，data.table 从未停止让我惊讶！ @Arun

【解决方案3】：

灵感来自大卫罗宾逊最初的尝试：

set.seed(1)
np <- 1000 # number of elements desired
M1 <- t(combn(1:np, 2))
sam <- sample(1:nrow(M1), np, replace = FALSE)
M2 <- M1[sam,]
anyDuplicated(M2) # returns FALSE

这将使用M1 的所有可能条目，但顺序是随机的。这是你想要的吗？

【讨论】：

我在编写我的解决方案时尝试了您的解决方案。但是，combn 给出一个大数字错误（例如，1e7）：test = combn(1e7,2) 给出这个错误：Error in matrix(r, nrow = len.r, ncol = count) : invalid 'ncol' value (too large or NA) In addition: Warning message: In combn(1e+07, 2) : NAs introduced by coercion
无赖。加上combn 即使是 10,000 点也很慢。精益代码就这么多！
是的，R 无法扩大规模...查看代码，expand.grid 似乎比 combn 快。 exapnd.grid 正在使用 data.frames 而 combn 正在使用矩阵。我想知道这是否是它更快的原因。

【解决方案4】：

这是我的尝试。它看起来不是很优雅，但它仍然比@Richard Erickson 的快一点（2.0s vs 2.6s，对于相同的尺寸）。这个想法是避免创建排列，因为这可能会花费大量时间并使用大量内存。相反，我在给定范围内创建了两个随机 ID 样本，并检查是否有任何行发生重复（这对于高范围和平均样本来说不太可能）。如果它们重复，则为第 2 列创建一个新样本并重复所有内容。

range <- 1e8
n <- 1e5
ids1 <- sample(range, n)
ids2 <- sample(range, n)
mat1 <- cbind(ids1, ids2)
found = FALSE
while(!found) {
  if (any(duplicated(rbind(mat1, mat1[,2:1])))) {
    ids2 <- sample(range, n)
    mat1 <- cbind(ids1, ids2)
  } else {
    found=TRUE
  }
}

【讨论】：

【解决方案5】：

怎么样：

no.pairs.needed <- 4 # or however many you want
npairs<-0
pairs <- NULL
top.sample.range <- 10000  # or whatever

while (npairs < no.pairs.needed){
  newpair <- matrix(data=sample(1:top.sample.range,2), nrow=1, ncol=2)
 if(!anyDuplicated(rbind(pairs, newpair))){
    pairs <- rbind(pairs, newpair)
    npairs <- npairs+1
  }
}

然后对象pairs 将返回您需要的矩阵。似乎可以扩展。

【讨论】：

【解决方案6】：

这是我的解决方案。

allIDX <- seq(10000000)
prtIDX <- sample(1:10000000, 10000000/2)
chlIDX <- allIDX[-prtIDX]
pairIDX <- cbind(prtIDX,chlIDX)

但我不必处理 10000000。

【讨论】：