使用C++生成具有权重的不重复随机整数

Question

我希望能够高效地生成一个非重复整数的随机样本，其范围为闭区间[0，rnd_max]，每个数字都可以选择，并且每个数字都与一个样本权重相关（权重越大，被选择的概率就越大，如果它在样本中尚未被选择，则概率恰好为weight[i]/sum(weight[not_taken])），

我了解到C ++有std::discrete_distribution可以生成随机加权整数，但是如果我使用它来生成随机整数并丢弃重复的整数，当要选取的样本相对于可能范围的长度很大时，将会有许多已经被选择的失败样本，导致高度低效的过程。不清楚Floyd算法是否有一些扩展到具有样本权重的情况（https://math.stackexchange.com/questions/178690/whats-the-proof-of-correctness-for-robert-floyds-algorithm-for-selecting-a-sin）- 我个人想不出来。

也可以例如使用将权重降为零的std::discrete_distribution，或执行部分加权洗牌，如在此答案中所示：C ++。加权std::shuffle - 但在该答案中，每次迭代都需要重新生成std::discrete_distribution，因此运行时间变为二次（它需要每次循环通过传递给它的权重）。

我想知道C++中什么是高效的加权随机样本以获得唯一整数，对于不同的样本大小（例如从可用范围中取1％到90％的样本）。

#include <vector>
#include <random>
#include <algorithm>

int main()
{
    size_t rnd_max = 1e5;
    size_t ntake = 1e3;

    unsigned int seed = 12345;
    std::mt19937 rng(seed);
    std::gamma_distribution<double> rgamma(1.0, 1.0);
    std::vector<double> weights(rnd_max);
    for (double &w : weights) w = rgamma(rng);

    std::vector<int> chosen_sample(ntake);
    // sampler goes here...

    return 0;
}

- anymous.asker

---

我对C++的分布不是很熟悉，所以我不知道有哪些。但我可以告诉你如何使用uniform_distribution自己实现，总时间复杂度为O(n log^2 n)（每个采样需要log^2 n的时间）。这对你有兴趣吗？ - user2956272

如果它们“不重复”，那么它们就不是随机的！ - Adrian Mole

@dyukha：好的，如果可以的话，麻烦您也这样做一下。@Adrian：是的，它们是这样的：想象一下以下过程：从一个空集合开始，然后使用“p[i] = {w[i] / sum(w[not taken]) if not taken, 0 otherwise}”连续添加元素——结果是随机不重复数字。 - anymous.asker

Answer 1

使用增广二叉搜索树可以优雅地解决这个问题。它提供了一个在 O(k log n) 时间内随机抽取 k 个元素的算法。

思路如下：假设你将所有元素按照排序顺序标记权重后存储在数组中，那么你可以通过以下方式（效率低下）解决此问题：

1. 生成介于 0 和所有元素总权重之间的随机数。 2. 迭代数组，直到找到一个元素，使得随机数处于该元素跨越的“范围”内。这里，“范围”表示从该元素开始到下一个元素开始的权重窗口。 3. 删除该元素并重复以上步骤。

如果按照上述实现，每次随机选择元素的时间复杂度为 O(n)：你必须迭代整个数组，然后在选择元素后删除某个单独的元素。这不是很好；总运行时间为 O(kn)。

我们可以稍微改进这个想法。在存储数组中的所有元素时，让每个元素存储其实际重量和其前面所有元素的组合重量。现在，要找到你要抽样的元素，你不需要使用线性搜索。相反，你可以在数组上使用二进制搜索，在时间复杂度为O(log n)的情况下找到你的元素。然而，这种方法的总运行时间仍然是每次迭代O(n)，因为这是选择的元素的成本，所以我们仍然处于O(kn)的领域。

然而，如果你将元素存储在一个平衡的二叉搜索树中，而不是在一个排序的数组中，其中每个元素都存储了其左子树中所有元素的重量，则可以模拟上述算法（二进制搜索被替换为对树的遍历）。此外，这具有一个优点，即从树中删除一个元素可以在O(log n)的时间内完成，因为它是一棵平衡的BST。

（如果您想知道如何进行遍历以找到所需的元素，请快速搜索“order statistics tree”。这里的想法本质上是这个想法的概括。）

按照@dyukha的建议，您可以通过在O(n)时间内从项目中构建完美平衡树（实际上这些项目不需要排序即可使用此技术-您明白原因吗？），然后每次需要删除某些内容时使用标准树删除算法来获得每个操作的O(log n)时间。这将给出总体解决方案运行时间为O(k log n)。

- templatetypedef

---

哦，不错！我有一个类似的想法，但我没有考虑到平衡树。我想使用二分查找+树状数组，这是O(log^2 n)。 - user2956272

@anymous.asker，平衡树可能会让人头疼，但你可以避免这种情况：你可以使用不平衡的二叉搜索树，并以随机顺序添加值到树中（因此首先进行洗牌，然后再添加）。结果树将具有高概率的平衡性。另一种选择是从一开始就构建一个完美平衡的树。 - user2956272

@dyukha 哦，从一开始就使用完全平衡的树，因为你只删除元素，所以无法增加高度的想法真的很好！我会编辑答案并包含这个建议。 :-) - templatetypedef

在不需要更新权重向量的情况下，最好将“树”以压平的形式存储 - 作为一个向量。您不会删除元素，而是将它们的权重暂时设置为零（每次选择样本整数时更新所有父级的权重总和; 最后您应该恢复初始值）。 - ALX23z

考虑提供伪代码来说明如何实现这个想法。另外，请注意C++包括 std::map，这是标准C++中最接近红黑树的东西。 - Peter O.

Answer 2

将答案转化为代码：

#include <iostream>
#include <vector>
#include <random>
#include <cmath>
#define pow2(n) ( 1 << (n) ) /* https://dev59.com/hnVD5IYBdhLWcg3wE3No */



int main()
{
    /* random and very biased set of weights */
    std::vector<double> weights{1, 1, 10000, 1, 30000, 1, 1, 500000};
    int rnd_max = weights.size();
    int ntake = 3;

    /* initialize random sampler */
    unsigned int seed = 12345;
    std::mt19937 rng(seed);

    /* determine smallest power of two that is larger than N */
    int tree_levels = ceil(log2((double) rnd_max));

    /* initialize vector with place-holders for perfectly-balanced tree */
    std::vector<double> tree_weights(pow2(tree_levels + 1));

    /* compute sums for the tree leaves at each node */
    int offset = pow2(tree_levels) - 1;
    for (int ix = 0; ix < rnd_max; ix++) {
        tree_weights[ix + offset] = weights[ix];
    }
    for (int ix = pow2(tree_levels+1) - 1; ix > 0; ix--) {
        tree_weights[(ix - 1) / 2] += tree_weights[ix];
    }

    /* sample according to uniform distribution */
    double rnd_subrange, w_left;
    double curr_subrange;
    int curr_ix;
    std::vector<int> sampled(ntake);
    for (int el = 0; el < ntake; el++) {

        /* go down the tree by drawing a random number and
           checking if it falls in the left or right sub-ranges */
        curr_ix = 0;
        curr_subrange = tree_weights[0];
        for (int lev = 0; lev < tree_levels; lev++) {
            rnd_subrange = std::uniform_real_distribution<double>(0, curr_subrange)(rng);
            w_left = tree_weights[2 * curr_ix + 1];
            curr_ix = 2 * curr_ix + 1 + (rnd_subrange >= w_left);
            curr_subrange = tree_weights[curr_ix];
        }

        /* finally, add element from this iteration */
        sampled[el] = curr_ix - offset;

        /* now remove the weight of the chosen element */
        tree_weights[curr_ix] = 0;
        for (int lev = 0; lev < tree_levels; lev++) {
            curr_ix = (curr_ix - 1) / 2;
            tree_weights[curr_ix] =   tree_weights[2 * curr_ix + 1]
                                    + tree_weights[2 * curr_ix + 2];
        }
    }

    std::cout << "sampled integers: [ ";
    for (int a : sampled) std::cout << a << " ";
    std::cout << "]" << std::endl;
    return 0;
}

从有偏差的权重中得到预期的输出：

sampled integers: [ 7 4 2 ]

（请注意，时间复杂度为O(n [使用节点权重构建树] + k * log2(n) [对元素进行采样]) - 比朴素的O(n * k)更好）

编辑：更新答案以适用于潜在的非唯一权重。

编辑2：对于更具数值稳健性的过程进行了小改动。

- anymous.asker