在R中生成不重复的随机整数对

Question

我希望能够生成不重复的随机整数对（换句话说，我不想有任何重复的整数对）。这个概念听起来很简单，但我想不出一个快速简便的解决方案。

举个例子，假设我想要使用整数序列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中绘制大量的配对时，情况就不太好了。

- Jacob H

---

如何在不重复的情况下获取 {3,3}？ - rawr

BrodieG，我正在一个项目的开头，所以我不确定这个序列有多大，可能不会是1e7。但比4要接近得多。 - Jacob H

@DavidRobinson，你的回答去哪了？ - Bryan Hanson

DavidRobinson的答案是错误的。 - Jacob H

你可以通过直接计算行和列来提高最后一行计算r.pair的性能，这是一个常数时间操作，取代了线性时间操作which：行是as.integer((x-1)%%4)+1L，列是as.integer((x-1)/4)+1L，其中x与sapply调用中的函数相同。 - Matthew Lundberg

Answer 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 个数字的 1e5 ^ 2 种可能的组合。这些 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

- BrodieG

---

你介意在你的基准测试中添加 CJ 方法吗？我很好奇它与其他方法相比如何。 - Richard Erickson

@RichardErickson，请查看更新后的内容（删除了data.table，这个奇怪的东西实际上增加了相当多的开销）。请注意，之前的答案已经在使用CJ，只是不必要地将其包装在data.table调用中（我猜这是有道理的，那会复制一个相当大的数据集）。 - BrodieG

谢谢！这是一个更好的方法。 - Jacob H

@JacobH，不用谢。另外请注意，我刚刚意识到这可以完全矢量化<headdesk>，可以再提高10-15倍的速度（请参见更新以消除sapply）。 - BrodieG

假设 size 为3，samples 为1。那么，vals 可以是9，因此你会分配 4, 0，但是4和0都不在范围1:3内。难道答案不应该是 ((vals-1) %/% size) + 1 和 (vals %% size) + 1 吗？ - Globe Theatre

Answer 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 

- Richard Erickson

---

你所链接的答案中的相关问题包含了许多有趣的方法和变化。显然这是一个具有挑战性的问题，而你做得很好！ - Bryan Hanson

谢谢！这太棒了。那个expand.grid函数非常方便。 - Jacob H

请使用CJ()代替expand.grid()。 - Arun

@Arun，CJ() 函数在哪个包中？?CJ 没有找到任何内容，而且在 Google 上搜索 "CJ" R 也没有找到任何有用的结果。 - Richard Erickson

为了以后参考，这里是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 - Richard Erickson

Answer 3

受David Robinson的初步尝试启发：

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的所有可能条目，但以随机顺序进行。这是你想要的吗？

- Bryan Hanson

---

我在写我的解决方案时尝试了你的解决方案。然而，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 - Richard Erickson

糟糕。即使对于10,000个点，combn也非常慢。这就是所谓的精简代码！ - Bryan Hanson

是的，“R” 在扩展方面可能效果不佳...浏览代码后，似乎“expand.grid”比“combn”更快。 “exapnd.grid”使用“data.frames”，而“combn”则使用矩阵。我想这就是它更快的原因。 - Richard Erickson

Answer 4

这是我的尝试。它看起来不太优雅，但对于相同的大小，它仍然比@Richard Erickson的快一点（2.0秒对2.6秒）。思路是避免创建排列，因为这可能需要很长时间并且使用大量内存。相反，我在给定范围内创建了两个随机ID样本，并检查是否存在任何重复行（对于高范围和平均样本来说，这是非常不可能的）。如果它们重复了，则创建第二列的新样本并重复所有步骤。

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

- Molx

Answer 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将返回您需要的矩阵。似乎可以很好地扩展。

- Scransom

Answer 6

这是我的解决方案。

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

但我不必处理10000000。

- skadoosh