无需替换的采样算法？
Posted
技术标签:
【中文标题】无需替换的采样算法？【英文标题】：Algorithm for sampling without replacement? 【发布时间】：2010-09-23 15:01:33 【问题描述】：我正在尝试测试特定数据聚类偶然发生的可能性。一个稳健的方法是蒙特卡罗模拟，其中数据和组之间的关联被随机重新分配大量次（例如 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 组的第二个元素，依此类推。指针重新分配很快，我会预先分配所有数据结构，但是没有替换的采样似乎是一个可能已经解决了很多次的问题。
首选逻辑或伪代码。
【问题讨论】：
【参考方案1】：这是一些基于 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。
【讨论】：
我还没有这本书，并且无法向自己证明算法的正确性。我在 java 中实现了它，并检查了人口项目是否以统一的概率进行抽样。结果令人信服。看到这个gist 在 Mathematica 中 Vitter 方法 D 的一个不加批判的实现比内置算法快几个数量级。我在这里描述它：tinyurl.com/lbldlpq @Alban - 我们可以通过考虑第一个元素来查看从 N 的总体中采样 n 个元素的问题。包含此元素的概率为 (n/N)：如果是，则问题归结为从剩余的 (N-1) 个元素中采样 (n-1) 个元素；如果不是，那么问题就归结为从剩余的 (N-1) 个元素中采样 (n) 个元素。一些变量变换会表明这是 Knuth 算法的精髓（通过增加 t）。u 处于开放、半开放还是封闭区间，(0, 1)，[0, 1) 或 [0, 1] 是否重要？ Knuth 只是说“均匀分布在零和一之间”。【参考方案2】：
基于answer by John D. Cook 的 C++ 工作代码。
#include <random>
#include <vector>
// John D. Cook, https://***.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;
【讨论】：
我编辑了您的答案，因此由于 GCC 和其他常见编译器中的线程保护，它不会太慢。但是，根据我的comment on John's answer，我不知道间隔应该是开放的、半开放的还是封闭的。这目前是半开放的。【参考方案3】：查看我对这个问题的回答Unique (non-repeating) random numbers in O(1)?。相同的逻辑应该可以完成您想要做的事情。
【讨论】：
太棒了！抱歉，当我搜索 SO（用于无替换抽样、统计、算法等）时，我没有看到该答案。也许这会像一个元问题一样引导像我这样的人找到你的原始答案。干杯！【参考方案4】：受@John D. Cook's answer 的启发，我用 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
【讨论】：
【参考方案5】：当总体规模远大于样本规模时，上述算法变得低效，因为它们具有复杂度O(n), n em> 是人口规模。
当我还是个学生的时候，我写了一些无放回均匀采样的算法，平均复杂度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)
该算法的复杂性在以下内容中进行了讨论： 鲁赞金，PS； Voytishek, A. V. 关于随机选择算法的成本。蒙特卡罗方法应用程序。 5 (1999), 没有。 1，39-54。 http://dx.doi.org/10.1515/mcma.1999.5.1.39
如果您觉得该算法有用，请提供参考。
另请参阅： P. Gupta, G. P. 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
【讨论】：
【参考方案6】：here 描述了另一种无需替换的采样算法。
这与 John D. Cook 在他的回答和 Knuth 中描述的相似，但它有不同的假设：总体规模未知，但样本可以放入内存中。这个叫做“Knuth's algorithm S”。
引用rosettacode文章：
选择可用的前 n 个项目作为样本； 对于 i > n 的第 i 个项目，有 n/i 个随机机会保留它。如果没有这个机会，样品将保持不变。如果 不是，让它随机 (1/n) 替换先前选择的 n 之一 样品的项目。 对任何后续项目重复 #2。
【讨论】：
Rosettacode 的算法名称错误：它应该是“算法 R”或“水库采样”。 “算法 S”（又名“选择抽样技术”）需要提前知道人口规模。 TAOCP - Vol 2 - §3.4.2 中描述了这两种算法【参考方案7】：我写了一个survey of algorithms for sampling without replacement。我可能有偏见，但我推荐我自己的算法，在下面用 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;
【讨论】：
与 std::sample 和 range::sample 相比如何？ 这将取决于您的特定 C++ 标准库如何实现它。在这两种情况下，文档都说“这个函数可以实现选择采样或水库采样”，所以它的执行可能与我对其中一种算法的实现类似，但您必须自己测试才能确定。以上是关于无需替换的采样算法？的主要内容，如果未能解决你的问题，请参考以下文章