在 R 中生成无替换的随机整数对
Generating Random Pairs of Integers without Replacement in R
我想绘制无替换的随机整数对(换句话说,我不想要任何重复的对)。这个概念听起来很简单,但我想不出一个快速简单的解决方案。
想象一下,例如,我想使用整数序列 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 中绘制大量对时,它就不太好了。
灵感来自 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
的所有可能条目,但顺序是随机的。这是你想要的吗?
首先,我找到了如何在 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
怎么样:
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
将 return 您需要的矩阵。似乎可以缩放。
这是我的尝试。它看起来不是很优雅,但它仍然比@Richard Erickson 的快一点(2.0s 对 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
}
}
这里的关键是不要生成所有排列,因为这在内存和时间方面非常昂贵。由于您只关心两个数字,只要 (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
这是我的解决方案。
allIDX <- seq(10000000)
prtIDX <- sample(1:10000000, 10000000/2)
chlIDX <- allIDX[-prtIDX]
pairIDX <- cbind(prtIDX,chlIDX)
但我不必处理 10000000。
我想绘制无替换的随机整数对(换句话说,我不想要任何重复的对)。这个概念听起来很简单,但我想不出一个快速简单的解决方案。
想象一下,例如,我想使用整数序列 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 中绘制大量对时,它就不太好了。
灵感来自 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
的所有可能条目,但顺序是随机的。这是你想要的吗?
首先,我找到了如何在 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
怎么样:
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
将 return 您需要的矩阵。似乎可以缩放。
这是我的尝试。它看起来不是很优雅,但它仍然比@Richard Erickson 的快一点(2.0s 对 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
}
}
这里的关键是不要生成所有排列,因为这在内存和时间方面非常昂贵。由于您只关心两个数字,只要 (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
这是我的解决方案。
allIDX <- seq(10000000)
prtIDX <- sample(1:10000000, 10000000/2)
chlIDX <- allIDX[-prtIDX]
pairIDX <- cbind(prtIDX,chlIDX)
但我不必处理 10000000。