带替换和不带替换的加权随机选择

Question

带替换和不带替换的加权随机选择

53

最近，我需要对列表进行带权重的随机选择，包括有放回和无放回的情况。虽然有一些已知且有效的算法可以进行等权重选择，以及一些可以进行带权重无放回选择的算法（如修改后的蓄水池算法），但我找不到任何好的算法可以进行带权重有放回选择。同时，我也想避免使用蓄水池方法，因为我需要选择列表中相当大的一部分，该部分足够小而可容纳在内存中。

请问在这种情况下最佳的解决方案是什么？我自己有一些解决方案，但我希望能找到更高效、更简单或两者兼备的东西。

- Arachnid

请查看此问题 https://dev59.com/bGLVa4cB1Zd3GeqPy7j8，即使它是 C# 而不是 Python，但只有很少的代码，所以任何人都应该能理解它。 - Omu

2

对于任何需要查询的人，“蓄水池算法”可以在维基百科的“蓄水池抽样”中找到。第一篇引用的论文是Jeffrey Scott Vitter的“Random Sampling with a Reservoir”，发表于1985年3月1日的《ACM Transactions on Mathematical Software》杂志上，卷11，号1，页码为37-57。 - Kevin J. Chase

对于不带重复的加权抽样，其中权重意味着被选择的概率与权重成比例，请参见我的答案：stackoverflow.com/a/27049669/262304 请注意，有些输入没有解决方案，例如从 {'a': 3, 'b': 1, 'c: 1} 中选择2个应该使'a'的出现频率是'b'或'c'的3倍，但这是不可能的。 - ech

我甚至不确定没有替换的加权选择是否被很好地定义。对于“从{'a': 3，'b': 1，'c': 1}中选择2个”的情况，一种可能的解释是a应该比b或c出现3倍，这是不可能的。另一种可能的解释是结果ab/ba/ac/ca应该具有相对权重31，而bc/cb应该具有相对权重11。还有一种可能的解释是我们通常选择第一个值，然后根据第一次未选择的任何内容的相对权重选择第二个值-这样ab比ba更有可能。 - Karl Knechtel

出于好奇，我进行了一些研究，并发现 numpy.choice 支持这种方式 - 采用后一种解释。 - Karl Knechtel

9个回答

14

这里没有提到的一种简单方法是由Efraimidis和Spirakis提出的。在Python中，您可以从权重为正且至少有m个加权项存储在weights中的n个加权项中选择m个项目，并返回所选索引，如下所示：

import heapq
import math
import random

def WeightedSelectionWithoutReplacement(weights, m):
    elt = [(math.log(random.random()) / weights[i], i) for i in range(len(weights))]
    return [x[1] for x in heapq.nlargest(m, elt)]

这种方法在结构上与Nick Johnson提出的第一种方法非常相似。不幸的是，那种方法在选择元素时存在偏差（请参见该方法的评论）。Efraimidis和Spirakis在相关论文中证明，他们的方法等同于无重复随机抽样。

- josliber

2

可以在此处找到经过优化（2.5k gas）的log2（0..1）的Solidity版本：https://gist.github.com/k06a/af6c58fe6634e48e53929451877eb5b5 - k06a

需要对权重进行归一化吗（总和为1）？ - undefined

1

@djstrong 不是的！ - undefined

5

以下是我为不重复的加权选择所想出的方法：

def WeightedSelectionWithoutReplacement(l, n):
  """Selects without replacement n random elements from a list of (weight, item) tuples."""
  l = sorted((random.random() * x[0], x[1]) for x in l)
  return l[-n:]

这是一个O(m log m)的算法，其中m是要从中选择的列表项数。我相当确信它会正确地对项目进行加权，尽管我没有以任何正式的方式验证过。

以下是我为带有替换的加权选择所想出的方法：

def WeightedSelectionWithReplacement(l, n):
  """Selects with replacement n random elements from a list of (weight, item) tuples."""
  cuml = []
  total_weight = 0.0
  for weight, item in l:
    total_weight += weight
    cuml.append((total_weight, item))
  return [cuml[bisect.bisect(cuml, random.random()*total_weight)] for x in range(n)]

这是O(m + n log m)的时间复杂度，其中m是输入列表中的项目数，n是要选择的项目数。

- Arachnid

6

那个第一个函数很棒，但它却没有正确地给项目分配权重。考虑 WeightedSelectionWithoutReplacement([(1, 'A'), (2, 'B')], 1)。它将有1/4的概率选择A，而不是1/3。难以修复。 - Jason Orendorff

顺便提一下，更快但更复杂的算法在我的答案中：https://dev59.com/D3I95IYBdhLWcg3wzhd9#2149533 - Jason Orendorff

很好的发现@JasonOrendorff。实际上差别相当大。对于权重（1、2、3、4），你期望“1”被选择1/10的时间，但它只会被选择1/94的时间。我真的希望那能行！ - Lawrence Kesteloot

1

@JasonOrendorff：你是怎么计算1/4的？如果你有一个公式，我们能否反转它，并用能够给出正确结果的权重替换原始权重？ - Lawrence Kesteloot

@LawrenceKesteloot - 对于1/4，我的看法是这样的：(random()*1)范围在0-1之间。如果它是1，则它比(random()*2)大的概率为1/2。如果它是0，则概率为0。平均概率为：1/4。 - ech

4

在创建一个额外的O(N)大小的数据结构并用O(N)时间后，可以在O(1)时间内进行带替换的加权随机选择。该算法基于Walker和Vose开发的Alias Method，这在这里中有很好的描述。

基本思想是通过均匀的RNG以概率1/N选择直方图中的每个bin。因此，我们将遍历它，并为任何接收到过多命中的低人口bin分配多余的命中到高人口bin。对于每个bin，我们存储属于它的命中百分比和超额的伙伴bin。此版本在原地跟踪小型和大型bin，消除了需要额外堆栈的需求。它使用伙伴的索引（存储在bucket[1]中）作为已经处理它们的指示器。

以下是基于这里的C实现的最小Python实现。

def prep(weights):
    data_sz = len(weights)
    factor = data_sz/float(sum(weights))
    data = [[w*factor, i] for i,w in enumerate(weights)]
    big=0
    while big<data_sz and data[big][0]<=1.0: big+=1
    for small,bucket in enumerate(data):
        if bucket[1] is not small: continue
        excess = 1.0 - bucket[0]
        while excess > 0:
            if big==data_sz: break
            bucket[1] = big
            bucket = data[big]
            bucket[0] -= excess
            excess = 1.0 - bucket[0]
            if (excess >= 0):
                big+=1
                while big<data_sz and data[big][0]<=1: big+=1
    return data

def sample(data):
    r=random.random()*len(data)
    idx = int(r)
    return data[idx][1] if r-idx > data[idx][0] else idx

使用示例：

TRIALS=1000
weights = [20,1.5,9.8,10,15,10,15.5,10,8,.2];
samples = [0]*len(weights)
data = prep(weights)

for _ in range(int(sum(weights)*TRIALS)):
    samples[sample(data)]+=1

result = [float(s)/TRIALS for s in samples]
err = [a-b for a,b in zip(result,weights)]
print(result)
print([round(e,5) for e in err])
print(sum([e*e for e in err]))

- AShelly

4

我建议您首先查看唐纳德·克努斯的半数值算法3.4.2节。如果您的数组很大，约翰·达普纳的随机变量生成原理第3章中有更有效的算法。如果您的数组不是非常大或者您不关心尽可能地提高效率，则克努斯的简单算法可能就足够了。

- John D. Cook

6

我刚刚看了一下3.4.2节，它只涵盖了有放回和无放回的无偏选择 - 没有提到加权选择。 - Nick Johnson

§3.4.1讨论了Walker的别名方法，该方法用于带替换的加权选择。 - rob mayoff

4

以下是关于集合（或允许重复的多重集合）中元素的随机加权选择的描述，包括有和没有替换，使用O(n)空间和O(log n)时间。

它由实现一个二叉搜索树组成，按要选择的元素排序，其中每个树的节点包含：

1. 元素本身（element） 2. 元素未归一化的权重（elementweight） 3. 左子节点及其所有子节点的未归一化权重之和（leftbranchweight） 4. 右子节点及其所有子节点的未归一化权重之和（rightbranchweight）

然后我们通过下降树来随机选择BST中的一个元素。算法的粗略描述如下。算法给出了树的节点。然后将节点的leftbranchweight、rightbranchweight和elementweight的值相加，并将权重除以这个和，得到leftbranchprobability、rightbranchprobability和elementprobability的值。然后获取一个介于0和1之间的随机数（randomnumber）。

如果数字小于elementprobability，则按正常方式从BST中删除元素，更新所有必要节点的leftbranchweight和rightbranchweight，并返回元素。

否则，如果数字小于（elementprobability + leftbranchweight），则在leftchild上递归（使用leftchild作为node运行算法）。

否则，在rightchild上递归。

当我们最终使用这些权重找到要返回的元素时，我们要么简单地返回它（带替换），要么删除它并更新树中相关的权重（不带替换）。

免责声明：该算法粗略，不试图对BST的正确实现进行论述；相反，希望这个答案能够帮助那些真正需要快速加权选择而没有替代方法的人（就像我一样）。

- djhaskin987

对于加权无替换算法，这会产生错误的结果。也就是说，元素不会按其权重比例被选择。请参见https://dev59.com/D3I95IYBdhLWcg3wzhd9#27049669以获取解决方案。 - ech

1

我已经在这里实现了那个算法（或稍作改进）https://github.com/makoConstruct/jostletree.rs（不过是用 Rust 写的）。 - mako

1

这是一个旧问题，但现在numpy已经提供了简单的解决方案，因此我想提一下。当前版本的numpy是1.2版本，numpy.random.choice允许进行有或无替换的采样，并给出权重。

- Joey Jojo Jr Shabadoo

0

假设您想从列表['white','blue','black','yellow','green']中使用概率分布[0.1, 0.2, 0.4, 0.1, 0.2]无重复地随机抽取3个元素。使用numpy.random模块非常简单：

    import numpy.random as rnd

    sampling_size = 3
    domain = ['white','blue','black','yellow','green']
    probs = [.1, .2, .4, .1, .2]
    sample = rnd.choice(domain, size=sampling_size, replace=False, p=probs)
    # in short: rnd.choice(domain, sampling_size, False, probs)
    print(sample)
    # Possible output: ['white' 'black' 'blue']

将replace标志设置为True，您可以进行带替换的抽样。

更多信息请参见： http://docs.scipy.org/doc/numpy/reference/generated/numpy.random.choice.html#numpy.random.choice

- Maroxo

0

我们在编程中遇到了一个问题，即每个时期按照他们的股份比例随机选择N个候选人中的K个验证者。但这会带来以下问题：

想象一下每个候选人的概率：

0.1
0.1
0.8

在进行100万次无放回的2选3后，每个候选人的概率如下：

0.254315
0.256755
0.488930

你应该知道，在无放回的情况下，原始概率是无法实现2/3选择的。

但我们希望初始概率成为收益分配概率。否则，小型候选池就会更具利润性。因此，我们意识到采用有放回的随机选择将有所帮助——可以随机选择N中的>K个并存储每个验证器的权重以进行奖励分配。

std::vector<int> validators;
std::vector<int> weights(n);
int totalWeights = 0;

for (int j = 0; validators.size() < m; j++) {
    int value = rand() % likehoodsSum;
    for (int i = 0; i < n; i++) {
        if (value < likehoods[i]) {
            if (weights[i] == 0) {
                validators.push_back(i);
            }
            weights[i]++;
            totalWeights++;
            break;
        }

        value -= likehoods[i];
    }
}

它在数百万个样本上提供了几乎原始的奖励分布：

0.101230
0.099113
0.799657

- k06a

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- John the Statistician · Accepted Answer

使用别名方法是从一个不变列表中进行替换样本的最快途径之一。核心直觉是我们可以为加权列表创建一组等大小的存储箱，通过位操作非常有效地进行索引，以避免二分查找。如果正确完成，则每个存储箱只需要从原始列表中存储两个项目，因此可以用单个百分比表示分割。

让我们以五个等权选择的例子（a:1, b:1, c:1, d:1, e:1）为例来创建别名查找：

将权重归一化，使它们总和为1.0。(a:0.2 b:0.2 c:0.2 d:0.2 e:0.2)这是选择每个权重的概率。
找到大于或等于变量数的最小2的幂，并创建此数量的分区|p|。每个分区表示1/|p|的概率质量。在本例中，我们创建了8个分区，每个分区能够包含0.125。
取剩余权重最小的变量，并尽可能多地将其质量放入空分区。在这个例子中，我们看到a填充了第一个分区。(p1{a|null,1.0},p2,p3,p4,p5,p6,p7,p8)，其中(a:0.075, b:0.2 c:0.2 d:0.2 e:0.2)
如果分区未填满，则取剩余权重最大的变量，并用该变量填充分区。

重复步骤3和4，直到不需要将原始分区的任何权重分配给列表。

例如，如果我们运行步骤3和4的另一个迭代，我们会看到

(p1{a|null,1.0},p2{a|b,0.6},p3,p4,p5,p6,p7,p8)剩下(a:0, b:0.15 c:0.2 d:0.2 e:0.2)未被分配。

在运行时：

获得一个U(0,1)的随机数，例如二进制0.001100000
将其左移lg2(p)位，找到索引分区。因此，我们将其左移3位，得到001.1，即位置1，因此是分区2。
如果分区被拆分，则使用移位后的随机数的小数部分来决定拆分。在这种情况下，该值为0.5，且0.5 < 0.6，因此返回a。

这里有一些代码和另一个解释,但不幸的是它没有使用位移技术，我也没有验证过。