我正在尝试测试偶然发生特定数据聚类的可能性。一种稳健的方法是蒙特卡罗模拟，其中数据和组之间的关联被随机重新分配大量次（例如 10,000 次），并且使用聚类度量来比较实际数据与模拟以确定 ap价值。
我已经完成了大部分工作，将分组映射到数据元素的指针，所以我计划随机重新分配指向数据的指针。问题：什么是无需替换的快速采样方法，以便在复制数据集中随机重新分配每个指针？
例如（这些数据只是一个简化的例子）：
数据（n=12 个值） - A 组：0.1、0.2、0.4 / B 组：0.5、0.6、0.8 / C 组：0.4、0.5 / D 组：0.2、0.2、0.3、0.5
对于每个复制数据集，我将拥有相同的集群大小（A=3、B=3、C=2、D=4）和数据值，但会将这些值重新分配给集群。
为此，我可以生成 1-12 范围内的随机数，分配 A 组的第一个元素，然后生成 1-11 范围内的随机数并分配 A 组的第二个元素，依此类推。指针重新分配很快，并且我将预先分配所有数据结构，但是没有替换的采样似乎是一个以前可能已经解决过很多次的问题。
逻辑或伪代码优先。
这是一些基于 Knuth 的书 Seminumeric Algorithms 的算法 3.4.2S 的无替换采样代码。
void SampleWithoutReplacement
(
int populationSize, // size of set sampling from
int sampleSize, // size of each sample
vector<int> & samples // output, zero-offset indicies to selected items
)
{
// Use Knuth's variable names
int& n = sampleSize;
int& N = populationSize;
int t = 0; // total input records dealt with
int m = 0; // number of items selected so far
double u;
while (m < n)
{
u = GetUniform(); // call a uniform(0,1) random number generator
if ( (N - t)*u >= n - m )
{
t++;
}
else
{
samples[m] = t;
t++; m++;
}
}
}
Jeffrey Scott Vitter 在“An Efficient Algorithm for Sequential Random Sampling”中提出了一种更有效但更复杂的方法，ACM Transactions on Mathematical Software, 13(1), March 1987, 58-67。
基于John D. Cook 回答的C++ 工作代码。
#include <random>
#include <vector>
// John D. Cook, https://stackoverflow.com/a/311716/15485
void SampleWithoutReplacement
(
int populationSize, // size of set sampling from
int sampleSize, // size of each sample
std::vector<int> & samples // output, zero-offset indicies to selected items
)
{
// Use Knuth's variable names
int& n = sampleSize;
int& N = populationSize;
int t = 0; // total input records dealt with
int m = 0; // number of items selected so far
std::default_random_engine re;
std::uniform_real_distribution<double> dist(0,1);
while (m < n)
{
double u = dist(re); // call a uniform(0,1) random number generator
if ( (N - t)*u >= n - m )
{
t++;
}
else
{
samples[m] = t;
t++; m++;
}
}
}
#include <iostream>
int main(int,char**)
{
const size_t sz = 10;
std::vector< int > samples(sz);
SampleWithoutReplacement(10*sz,sz,samples);
for (size_t i = 0; i < sz; i++ ) {
std::cout << samples[i] << "\t";
}
return 0;
}
请参阅我对这个问题的回答O(1) 中的唯一（非重复）随机数？. 相同的逻辑应该完成您想要做的事情。
受@John D. Cook 的回答启发，我在 Nim 中编写了一个实现。起初我很难理解它是如何工作的，所以我进行了广泛的评论，还包括一个例子。也许它有助于理解这个想法。另外，我稍微更改了变量名称。
iterator uniqueRandomValuesBelow*(N, M: int) =
## Returns a total of M unique random values i with 0 <= i < N
## These indices can be used to construct e.g. a random sample without replacement
assert(M <= N)
var t = 0 # total input records dealt with
var m = 0 # number of items selected so far
while (m < M):
let u = random(1.0) # call a uniform(0,1) random number generator
# meaning of the following terms:
# (N - t) is the total number of remaining draws left (initially just N)
# (M - m) is the number how many of these remaining draw must be positive (initially just M)
# => Probability for next draw = (M-m) / (N-t)
# i.e.: (required positive draws left) / (total draw left)
#
# This is implemented by the inequality expression below:
# - the larger (M-m), the larger the probability of a positive draw
# - for (N-t) == (M-m), the term on the left is always smaller => we will draw 100%
# - for (N-t) >> (M-m), we must get a very small u
#
# example: (N-t) = 7, (M-m) = 5
# => we draw the next with prob 5/7
# lets assume the draw fails
# => t += 1 => (N-t) = 6
# => we draw the next with prob 5/6
# lets assume the draw succeeds
# => t += 1, m += 1 => (N-t) = 5, (M-m) = 4
# => we draw the next with prob 4/5
# lets assume the draw fails
# => t += 1 => (N-t) = 4
# => we draw the next with prob 4/4, i.e.,
# we will draw with certainty from now on
# (in the next steps we get prob 3/3, 2/2, ...)
if (N - t)*u >= (M - m).toFloat: # this is essentially a draw with P = (M-m) / (N-t)
# no draw -- happens mainly for (N-t) >> (M-m) and/or high u
t += 1
else:
# draw t -- happens when (M-m) gets large and/or low u
yield t # this is where we output an index, can be used to sample
t += 1
m += 1
# example use
for i in uniqueRandomValuesBelow(100, 5):
echo i
当总体规模远大于样本规模时，上述算法变得低效，因为它们的复杂度为O ( n )，n是总体规模。
当我还是学生的时候，我写了一些无需替换的均匀采样算法，它们的平均复杂度为O ( s log s )，其中s是样本大小。这是二叉树算法的代码，平均复杂度为O ( s log s )，在 R 中：
# The Tree growing algorithm for uniform sampling without replacement
# by Pavel Ruzankin
quicksample = function (n,size)
# n - the number of items to choose from
# size - the sample size
{
s=as.integer(size)
if (s>n) {
stop("Sample size is greater than the number of items to choose from")
}
# upv=integer(s) #level up edge is pointing to
leftv=integer(s) #left edge is poiting to; must be filled with zeros
rightv=integer(s) #right edge is pointig to; must be filled with zeros
samp=integer(s) #the sample
ordn=integer(s) #relative ordinal number
ordn[1L]=1L #initial value for the root vertex
samp[1L]=sample(n,1L)
if (s > 1L) for (j in 2L:s) {
curn=sample(n-j+1L,1L) #current number sampled
curordn=0L #currend ordinal number
v=1L #current vertice
from=1L #how have come here: 0 - by left edge, 1 - by right edge
repeat {
curordn=curordn+ordn[v]
if (curn+curordn>samp[v]) { #going down by the right edge
if (from == 0L) {
ordn[v]=ordn[v]-1L
}
if (rightv[v]!=0L) {
v=rightv[v]
from=1L
} else { #creating a new vertex
samp[j]=curn+curordn
ordn[j]=1L
# upv[j]=v
rightv[v]=j
break
}
} else { #going down by the left edge
if (from==1L) {
ordn[v]=ordn[v]+1L
}
if (leftv[v]!=0L) {
v=leftv[v]
from=0L
} else { #creating a new vertex
samp[j]=curn+curordn-1L
ordn[j]=-1L
# upv[j]=v
leftv[v]=j
break
}
}
}
}
return(samp)
}
该算法的复杂性在：Rouzankin, PS；Voytishek, AV 关于随机选择算法的成本。蒙特卡罗方法应用程序。5 (1999), 没有。1，39-54。 http://dx.doi.org/10.1515/mcma.1999.5.1.39
如果您发现该算法有用，请提供参考。
另见：P. Gupta，GP Bhattacharjee。(1984) 一种有效的无替换随机抽样算法。国际计算机数学杂志 16:4，第 201-209 页。DOI: 10.1080/00207168408803438
Teuhola, J. 和 Nevalainen, O. 1982。两种有效的无替换随机抽样算法。/IJCM/, 11(2): 127–140。DOI: 10.1080/00207168208803304
在上一篇论文中，作者使用哈希表并声称他们的算法具有O ( s ) 复杂度。还有一种更快的哈希表算法，很快就会在 pqR（pretty quick R）中实现： https ://stat.ethz.ch/pipermail/r-devel/2017-October/075012.html
此处描述了另一种无需替换的采样算法。
它类似于 John D. Cook 在他的回答和 Knuth 中描述的那个，但它有不同的假设：人口规模未知，但样本可以适合记忆。这称为“Knuth 算法 S”。
引用 Rosettacode 文章：
- 在可用时选择前 n 个项目作为样本；
- 对于 i > n 的第 i 个项目，有 n/i 个随机机会保留它。如果没有这个机会，样品将保持不变。如果没有，让它随机 (1/n) 替换先前选择的 n 个样本项目之一。
- 对任何后续项目重复 #2。
我写了一份关于无放回抽样算法的调查。我可能有偏见，但我推荐我自己的算法，在下面用 C++ 实现，为许多 k、n 值提供最佳性能，并为其他值提供可接受的性能。randbelow(i)假定返回一个公平选择的小于 的随机非负整数i。
void cardchoose(uint32_t n, uint32_t k, uint32_t* result) {
auto t = n - k + 1;
for (uint32_t i = 0; i < k; i++) {
uint32_t r = randbelow(t + i);
if (r < t) {
result[i] = r;
} else {
result[i] = result[r - t];
}
}
std::sort(result, result + k);
for (uint32_t i = 0; i < k; i++) {
result[i] += i;
}
}