具有权重的 C++ 随机非重复整数答案

【问题标题】：C++ random non-repeated integers with weights具有权重的 C++ 随机非重复整数
【发布时间】：2019-08-21 21:45:27
【问题描述】：

我想在（封闭）范围[0, rnd_max] 内有效地生成唯一（非重复）整数的随机样本，范围内的每个数字都可以选择，并且每个数字都与样本权重相关联（权重越大，该数字被选中的可能性就越大，如果该数字尚未包含在样本中，则概率恰好是下一个选择 weight[i] / sum(weight[not_taken])）。

我看到C++有std::discrete_distribution可以生成随机加权整数，但是如果我用它来生成随机整数并丢弃重复的整数，当要取的样本相对于可能范围的长度很大时，就会有已经采集了许多不合格的样本，导致程序效率极低。我不清楚弗洛伊德的算法是否对样本权重的情况有一些扩展 (https://math.stackexchange.com/questions/178690/whats-the-proof-of-correctness-for-robert-floyds-algorithm-for-selecting-a-sin) - 我个人想不出一个。

也可以例如使用std::discrete_distribution 将权重降至零，或执行部分加权洗牌，如此答案：C++. Weighted 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;
}

【问题讨论】：

我对 C++ 发行版不太熟悉，所以我不知道。我可以告诉你如何使用uniform_distribution 在O(n log^2 n) 总时间（每次采样的log^2 n 时间）中自己实现它。你感兴趣吗？
如果它们“不重复”，那么它们就不是随机的！
@dyukha ：是的，拜托，那也很棒。 @Adrian：是的，他们是：想象以下过程：从一个空集开始，然后使用 p[i] = {w[i] / sum(w[not taken]) if not taken, 0 otherwise} 按顺序添加元素 - 结果是随机的非重复数字。

标签： c++ random

【解决方案1】：

有一个很好的方法可以使用增强的二叉搜索树来解决这个问题。它给出了一个 O(k log n) 时间的算法来随机采样 k 个元素。

这个想法是这样的。假设您将所有元素按排序顺序存储在一个数组中，每个元素都标有其权重。然后，您可以按如下方式解决此问题（效率低下）：

在 0 和所有元素的总权重之间生成一个随机数。
遍历数组，直到找到一个元素，使得随机数在该元素跨越的“范围”内。这里，“范围”表示从该元素开始到下一个元素开始的权重窗口。
删除该元素并重复。

如果你按照上面提到的方法实现，每次选择一个随机元素都需要时间 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) 的整体解决方案运行时间。

【讨论】：

哦，不错！我有类似的总体想法，但我没有考虑平衡树。我想用二分查找+fenwick tree，也就是O(log^2 n)。
@anymous.asker，平衡树可能会很痛苦，但你可以避免它：你可以使用不平衡的 BST 并以随机顺序向树添加值（所以先洗牌，然后再添加） .结果树将以高概率平衡。另一种选择是从一开始就构建一个完美平衡的树。
@dyukha 哦，从一开始就使用完美平衡的树的想法，因为你只是在删除东西，因此不能增加高度，这是一个非常好的想法！我将编辑答案以包含它。 :-)
@anymous.asker 在不需要更新权重向量的情况下，最好将“树”存储在扁平版本中 - 作为向量。您不会删除元素，而是暂时将它们的权重设置为零（并在每次选择样本整数时更新其所有父项的权重总和；最后您应该恢复初始值）。
考虑提供伪代码来说明如何实现这个想法。另外请注意，C++ 包括std::map，它在本质上最接近标准 C++ 中的红黑树。

【解决方案2】：

将答案放入代码中：

#include <iostream>
#include <vector>
#include <random>
#include <cmath>
#define pow2(n) ( 1 << (n) ) /* https://stackoverflow.com/questions/101439/the-most-efficient-way-to-implement-an-integer-based-power-function-powint-int */



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 [when building the tree with sums of nodes weights] + k * log2(n) [when sampling the elements])——比天真的O(n * k)好）

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

EDIT2：对数值更稳健的过程进行小改动。

【讨论】：