从数组中加权随机选择

Question

我想从一个数组中随机选择一个元素，但每个元素都有已知的选择概率。

所有概率（在数组内）之和为1。

您会建议哪个算法最快、最适合进行大量计算？

例如：

id => chance
array[
    0 => 0.8
    1 => 0.2
]

对于这个伪代码，所讨论的算法在多次调用时应该会以 0 id 上的四个元素对应一个 1 id 的元素统计返回。

- Mikulas Dite

Answer 1

计算离散累积概率密度函数（CDF）的列表 - 或者简单地说是权重的累积和数组。然后在0到所有权重之和（如果你的情况下可能为1）的范围内生成一个随机数，使用二分搜索在你的离散CDF数组中找到这个随机数并获取对应于该条目的值 - 这就是您的加权随机数。

- Sven Marnach

---

这个二分查找每次查询需要进行9步操作，计算公式为log2(500) = 9。 - thejh

生成0到权重总和之间的随机数，谁能保证生成的随机数将在累积分布函数数组中？假设我们有[0.1 0.2 0.4 0.3]作为权重数组。累积分布函数数组将是[0.1 0.3 0.7 1.0]。rand值必须在0到1.0之间生成。然后可以是例如0.62，但该值不在累积分布函数数组中。 - Mazzy

@Mazzy：您正在寻找包含您生成的随机数的区间 - 在这种情况下，该区间为0.3到0.7。当然，您不能期望出现精确值，但是使用二分查找来找到区间仍然有效。 - Sven Marnach

@Mazzy：二分查找可以轻松地用于查找您要查找的值所在的区间，这就是您所需要的全部。大多数编程语言标准库中的二分查找实现不需要找到确切的值，例如C++中的lower_bound()或Python中的bisect_left()。 - Sven Marnach

我无法完全理解Sven Marnach的方法。@SvenMarnach - 以A = [1,2,3]，Pa = [0.2，0.6，0.2]为例=意味着我想在长期内看到2三倍于1或3。在这种情况下，cdf（Pa）= [0.2，0.8，1.0]现在，在0和1之间生成一个随机数，假设我得到0.9。 我如何进行二进制搜索？ 我选择2还是3？请澄清一下。 您的解决方案目前看起来最好。谢谢 - Karan Kapoor

Answer 2

这个算法很直接

rand_no = rand(0,1)
for each element in array 
     if(rand_num < element.probablity)
          select and break
     rand_num = rand_num - element.probability

- Rohit J

---

这不会起作用，因为我有机会，而不是领域。| 即使有人对这个答案进行了负面评价，它给了我一个可行的想法。限制非常简单地计算，不应影响性能。 - Mikulas Dite

@Mikulas 假设您有离散的机会和在0到1之间均匀分布的随机数，它将给出与其权重相等的概率。对于您的情况，有80％的机会随机数小于.8，因此将选择第一个元素，20％的机会大于.8，在这种情况下将选择第二个元素。 - Rohit J

不需要排序也可以工作，如果您想在选择元素后删除它，则比二分查找更快。 - Rohit J

抱歉，我的问题是如果我有两个重量相同的元素怎么办？在这种情况下，我只会得到数组中两个元素中的第一个吗？还是我错了？ 如果有两个重量相同的元素，你只会获得这两个元素中的第一个。 - arpho

@arpho 我在 JavaScript 中测试了你的假设（https://jsfiddle.net/mwbrympx/1/）。看起来你是错的。 - 4castle

Answer 3

我发现这篇文章在理解这个问题方面最有用。 Stackoverflow上的这个问题也许是你要找的东西。

---

我相信最优的解决方案是使用Alias Method（维基百科）。 它需要O(n)时间进行初始化，O(1)时间进行选择，并且需要O(n)内存。

以下是生成投掷加权n-面骰子结果的算法（从此处选择长度为n的数组元素非常简单），摘自这篇文章。作者假设您有一个用于投掷公平骰子的函数（floor(random() * n)）和翻转偏倚硬币的函数（random() < p）。

算法：Vose's Alias Method

初始化：

1. 创建大小为n的数组Alias和Prob。
2. 创建两个工作列表Small和Large。
3. 将每个概率乘以n。
4. 对于每个缩放后的概率p_i：
   1. 如果p_i<1，则将i添加到Small。
   2. 否则（p_i ≥1 ），将i添加到Large。
5. 当Small和Large都不为空时（可能先清空Large）：
   
   1. 从Small中移除第一个元素，并将其称为l。
   2. 从Large中移除第一个元素，并将其称为g。
   3. 设置Prob[l] = p_l。
   4. 设置Alias[l] = g。
   5. 设置p_g := (p_g+p_l)−1。（这是一种更稳定的数值计算方式）
   6. 如果p_g<1，则将g添加到Small中。
   7. 否则（p_g ≥ 1），将g添加到Large中。
6. 当Large不为空时：
   
   1. 从Large中移除第一个元素，并将其称为g。
   2. 设置Prob[g] = 1。
7. 当Small不为空时：这只可能是由于数值不稳定引起的。
   
   1. 从Small中移除第一个元素，并将其称为l。
   2. 设置Prob[l] = 1。

生成：

1. 从一个有n个面的公正骰子上生成一个随机数，称其为第i面。
2. 投掷一个有偏的硬币，其正面朝上的概率为 Prob[i]。
3. 如果硬币正面朝上，则返回 i。
4. 否则，返回 Alias[i]。

- Simon Baumgardt-Wellander

Answer 4

这是一个Ruby的实现例子：

def weighted_rand(weights = {})
  raise 'Probabilities must sum up to 1' unless weights.values.inject(&:+) == 1.0
  raise 'Probabilities must not be negative' unless weights.values.all? { |p| p >= 0 }
  # Do more sanity checks depending on the amount of trust in the software component using this method,
  # e.g. don't allow duplicates, don't allow non-numeric values, etc.
  
  # Ignore elements with probability 0
  weights = weights.reject { |k, v| v == 0.0 }   # e.g. => {"a"=>0.4, "b"=>0.4, "c"=>0.2}

  # Accumulate probabilities and map them to a value
  u = 0.0
  ranges = weights.map { |v, p| [u += p, v] }   # e.g. => [[0.4, "a"], [0.8, "b"], [1.0, "c"]]

  # Generate a (pseudo-)random floating point number between 0.0(included) and 1.0(excluded)
  u = rand   # e.g. => 0.4651073966724186
  
  # Find the first value that has an accumulated probability greater than the random number u
  ranges.find { |p, v| p > u }.last   # e.g. => "b"
end

使用方法：

weights = {'a' => 0.4, 'b' => 0.4, 'c' => 0.2, 'd' => 0.0}

weighted_rand weights

大致预期内容如下：

sample = 1000.times.map { weighted_rand weights }
sample.count('a') # 396
sample.count('b') # 406
sample.count('c') # 198
sample.count('d') # 0

- knugie

---

刚刚用了这个，发现我认识这个名字！谢谢@wolfgang-teuber！ - Abe Petrillo

使用这种方法需要注意的一点是，如果您的权重为1.0，其余为0.0，则此方法将不能按预期工作。我们的权重是以环境变量的形式存在的，当我们将其中一个权重设置为1.0（即始终为真）时，它产生了相反的影响。这只是给其他人提个醒！ - Abe Petrillo

@AbePetrillo 我更新了 weighted_rand 方法以解决你所描述的问题。 - knugie

Answer 5

一个ruby的例子

#each element is associated with its probability
a = {1 => 0.25 ,2 => 0.5 ,3 => 0.2, 4 => 0.05}

#at some point, convert to ccumulative probability
acc = 0
a.each { |e,w| a[e] = acc+=w }

#to select an element, pick a random between 0 and 1 and find the first   
#cummulative probability that's greater than the random number
r = rand
selected = a.find{ |e,w| w>r }

p selected[0]

- krusty.ar

---

在这个算法中，最后一个元素永远不会被选中，因为它的概率是1.0，而rand将始终在0和1之间。 - Matt Darby

Answer 6

按如下方式，每个样本的期望时间为O(1)。

计算每个元素i的累积分布函数CDF F(i)，即小于或等于i的概率之和。

将元素i的范围r(i)定义为区间[F(i-1),F(i)]。

对于每个区间[(i-1)/n,i/n]，创建一个桶，其中包含范围与该区间重叠的元素列表。只要你足够小心，这需要总共O(n)的时间来处理完整的数组。

当你随机抽取数组时，只需计算随机数所在的桶，并与列表中的每个元素进行比较，直到找到包含它的区间。

一个样本的成本是O(随机选择的列表的期望长度) <= 2。

- jonderry

---

如果权重相差悬殊，该算法的最坏时间复杂度为O(n)。可能会出现所有区间都属于同一个桶的情况。如果没有对权重进行额外限制，这绝对不是O(1)，甚至不是O(log n)。 - Sven Marnach

最坏情况很少发生。如果所有的n个区间都重叠在一个桶中，那么几乎所有的查询都只需要与一个区间进行比较。在实践中，这将比二分搜索快得多。如果您坚持优化最坏情况，可以在每个桶内执行二分搜索，使得每个查询的成本在最坏情况下为O(lg(最大桶长度))，在期望中为O(随机选择列表的lg长度的期望值)，仍然是O(1)。 - jonderry

谢谢，看起来很不错。我需要运行一些试验才能确定它是否比我的解决方案中的CDF方法更快。 - Mikulas Dite

平均而言，它是O(1)，这是真的。+1 - Sven Marnach

@Mikulas Dite，值得强调的是，这也是一种CDF数组解决方案，与纯二分查找的区别有点像在数组中使用二分查找和哈希搜索元素之间的区别。另一种看待它的方式是计算CDF数组，而不是在其上进行二分查找，而是将随机数哈希到对应于桶起始位置的数组索引。然后，您可以使用任何搜索策略（例如，暴力线性搜索或二分搜索）进一步缩小范围以找到正确的采样元素。 - jonderry

请注意，您在这里的保证比通常的“最坏情况”评估要好，因为通过构造，您的访问是已知为随机的... - comingstorm

Answer 7

这是我在生产环境中使用的PHP代码：

/**
 * @return \App\Models\CdnServer
*/
protected function selectWeightedServer(Collection $servers)
{
    if ($servers->count() == 1) {
        return $servers->first();
    }

    $totalWeight = 0;

    foreach ($servers as $server) {
        $totalWeight += $server->getWeight();
    }

    // Select a random server using weighted choice
    $randWeight = mt_rand(1, $totalWeight);
    $accWeight = 0;

    foreach ($servers as $server) {
        $accWeight += $server->getWeight();

        if ($accWeight >= $randWeight) {
            return $server;
        }
    }
}

- Agent Coop

Answer 8

使用pickup gem的Ruby解决方案：

require 'pickup'

chances = {0=>80, 1=>20}
picker = Pickup.new(chances)

例子：

5.times.collect {
  picker.pick(5)
}

输出如下：

[[0, 0, 0, 0, 0], 
 [0, 0, 0, 0, 0], 
 [0, 0, 0, 1, 1], 
 [0, 0, 0, 0, 0], 
 [0, 0, 0, 0, 1]]

- devstopfix

Answer 9

在我看来，最有效的算法是为数组的每个元素生成一个随机数，该随机数从具有由该元素权重给定的参数的指数分布中绘制。当你遍历数组时，保留具有最低“排序数字”的元素。在这种情况下，特定元素具有数组最低排序号的概率与数组元素的权重成比例。

详细信息和代码如下。

该算法是O(n)，不涉及排序或额外存储，并且选择可以在一次遍历数组中完成。权重必须大于零，但不必总和为任何特定值。

附加功能：如果您将每个数组元素的排序号码存储在其中，则可以选择按增加的排序号对数组进行排序，以获得数组的随机排序，在其中具有更高权重的元素很有可能较早出现（我发现这在决定要选择哪个DNS SRV记录以决定查询哪台计算机时非常有用）。

其他算法：重复随机抽样需要每次通过数组进行新的遍历；对于无重复随机选择，可以按照递增排序号的顺序对数组进行排序，并以此顺序读取k个元素。

请参阅维基百科关于指数分布的页面（特别是有关这种变量集合最小值分布的注释），证明上述结论的正确性，同时也为生成这种变量的技术提供了指针：如果T在[0,1)内具有均匀随机分布，则Z=-log(1-T)/w（其中w是分布的参数；这里是相关元素的权重）具有指数分布。

也就是说：

对于数组中的每个元素i，计算zi = -log(1-T)/wi，其中T是从[0,1)均匀分布中抽取的值，wi是第i个元素的权重。

在遍历数组时，保留到目前为止具有最低zi值的元素的引用。

元素i将以概率wi/(w1+w2+...+wn)被选中。

以下是Python中的示例，它对权重数组进行一次遍历，进行10000次试验。

import math, random

random.seed()

weights = [10, 20, 50, 20]
nw = len(weights)
results = [0 for i in range(nw)]

n = 10000
while n > 0: # do n trials
    smallest_i = 0
    smallest_z = -math.log(1-random.random())/weights[0]
    for i in range(1, nw):
        z = -math.log(1-random.random())/weights[i]
        if z < smallest_z:
            smallest_i = i
            smallest_z = z

    # we have selected element 'smallest_i'
    results[smallest_i] += 1 # accumulate our choices

    n -= 1

for i in range(nw):
    print("{} -> {}".format(weights[i], results[i]))

---

编辑（历史记录）：发布后，我确信自己不可能是第一个想到这个解决方案的人，再次搜索后发现确实如此。

• 在类似问题的回答中，Joe K建议了这个算法（并且还指出可能有其他人先想到了这个算法）。
• 与此同时，那个问题的另一个回答则指向了Efraimidis和Spirakis（预印本），该文描述了一种类似的方法。
• 我非常确定，看起来Efraimidis和Spirakis实际上是同一个指数分布算法的伪装形式，并且这得到了维基百科关于蓄水池抽样的页面的一句话评论的证实：“等价地，这个算法的一个更加数值稳定的表述就是上面的指数分布算法”。该引用来自Richard Arratia的一系列讲义；指数分布的相关属性在第1.3节中提到（其中提到类似的内容在某些领域是“熟知的事实”），但没有提到它与Efraimidis和Spirakis算法的关系。

- Norman Gray

Answer 10

“幸运轮”O(n)算法，仅适用于小数组：

function pickRandomWeighted(array, weights) {
    var sum = 0;
    for (var i=0; i<weights.length; i++) sum += weights[i];
    for (var i=0, pick=Math.random()*sum; i<weights.length; i++, pick-=weights[i])
        if (pick-weights[i]<0) return array[i];
}

- Sarsaparilla