无需替换的采样算法？

Question

我正在尝试测试特定数据聚类偶然发生的可能性。一个稳健的方法是蒙特卡罗模拟，其中数据和组之间的关联被随机重新分配大量次（例如 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 组的第二个元素，依此类推。指针重新分配很快，我会预先分配所有数据结构，但是没有替换的采样似乎是一个可能已经解决了很多次的问题。

首选逻辑或伪代码。

Answer 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。

Answer 2

基于answer by 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;
}

Answer 3

查看我对这个问题的回答Unique (non-repeating) random numbers in O(1)?。相同的逻辑应该可以完成您想要做的事情。

Answer 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

Answer 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

Answer 6

here 描述了另一种无需替换的采样算法。

这与 John D. Cook 在他的回答和 Knuth 中描述的相似，但它有不同的假设：总体规模未知，但样本可以放入内存中。这个叫做“Knuth's algorithm S”。

引用rosettacode文章：

1. 选择可用的前 n 个项目作为样本；
2. 对于 i > n 的第 i 个项目，有 n/i 个随机机会保留它。如果没有这个机会，样品将保持不变。如果 不是，让它随机 (1/n) 替换先前选择的 n 之一 样品的项目。
3. 对任何后续项目重复 #2。

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