【问题标题】：Randomly Generating Combinations From Variable Weights从可变权重随机生成组合
【发布时间】：2015-09-10 06:07:51
【问题描述】：

非常重要的编辑：所有A_i都是独特的。 p>

问题

我有 A 个 n 唯一个对象的列表。每个对象A_i都有一个可变的百分比P_i。

我想创建一个算法来生成 k 个对象的新列表 B (k n/₂ 并且在大多数情况下 k 明显小于 ⁿ/₂。例如 n=231 , k=21)。列表 B 应该没有重复项，并且将填充来自列表 A 的对象，但具有以下限制：

对象A_i出现在B中的概率是P_i.

我尝试过的

（这些代码在 PHP 中仅用于测试目的）我首先列出了A

$list = [
    "A" => 2.5, 
    "B" => 2.5, 
    "C" => 2.5, 
    "D" => 2.5, 
    "E" => 2.5, 
    "F" => 2.5, 
    "G" => 2.5, 
    "H" => 2.5, 
    "I" => 5,   
    "J" => 5,   
    "K" => 2.5, 
    "L" => 2.5, 
    "M" => 2.5, 
    "N" => 2.5, 
    "O" => 2.5, 
    "P" => 2.5, 
    "Q" => 2.5, 
    "R" => 2.5, 
    "S" => 2.5, 
    "T" => 2.5, 
    "U" => 5,   
    "V" => 5,   
    "W" => 5,   
    "X" => 5,   
    "Y" => 5,   
    "Z" => 20   
];

一开始我尝试了以下两种算法（这些在 PHP 中只是为了测试目的）：

$result = [];

while (count($result) < 10) {
    $rnd = rand(0,10000000) / 100000;

    $sum = 0;
    foreach ($list as $key => $value) {
        $sum += $value;
        if ($rnd <= $sum) {
            if (in_array($key,$result)) {
                break;
            } else {
                $result[] = $key;
                break;
            }
        }
    }
}

与

$result = [];

while (count($result) < 10) {
    $sum = 0;
    foreach ($list as $key => $value) {
        $sum += $value;
    }

    $rnd = rand(0,$sum * 100000) / 100000;

    $sum = 0;
    foreach ($list as $key => $value) {
        $sum += $value;
        if ($rnd <= $sum) {
            $result[] = $key;
            unset($list[$key]);
            break;
        }
    }
}

这两种算法的唯一区别是遇到重复时会再次尝试，而当对象表单列表A被拾取时会删除它。事实证明，这两种算法具有相同的概率输出。

我运行第二个算法 100,000 次，并记录每个字母被选中的次数。以下数组包含基于 100,000 次测试在任何列表 B 中选择一个字母的概率。

[A] => 30.213
[B] => 29.865
[C] => 30.357
[D] => 30.198
[E] => 30.152
[F] => 30.472
[G] => 30.343
[H] => 30.011
[I] => 51.367
[J] => 51.683
[K] => 30.271
[L] => 30.197
[M] => 30.341
[N] => 30.15
[O] => 30.225
[P] => 30.135
[Q] => 30.406
[R] => 30.083
[S] => 30.251
[T] => 30.369
[U] => 51.671
[V] => 52.098
[W] => 51.772
[X] => 51.739
[Y] => 51.891
[Z] => 93.74

回顾算法时，这是有道理的。该算法错误地将原始百分比解释为在任何给定位置而不是任何列表 B 中拾取对象的概率百分比。例如，实际上，在列表 B 中选择 Z 的机会是 93%，但在索引 B_n 中选择 Z 的机会 是 20%。这不是我想要的。我希望 Z 在列表 B 中被选中的机会为 20%。

这甚至可能吗？怎么办？

编辑 1

我尝试简单地让所有 P_i = k 的总和，如果所有 P_i是相等的，但是修改了它们的值之后，它开始变得越来越错误。

初始概率

$list= [
    "A" => 8.4615,
    "B" => 68.4615,
    "C" => 13.4615,
    "D" => 63.4615,
    "E" => 18.4615,
    "F" => 58.4615,
    "G" => 23.4615,
    "H" => 53.4615,
    "I" => 28.4615,
    "J" => 48.4615,
    "K" => 33.4615,
    "L" => 43.4615,
    "M" => 38.4615,
    "N" => 38.4615,
    "O" => 38.4615,
    "P" => 38.4615,
    "Q" => 38.4615,
    "R" => 38.4615,
    "S" => 38.4615,
    "T" => 38.4615,
    "U" => 38.4615,
    "V" => 38.4615,
    "W" => 38.4615,
    "X" => 38.4615,
    "Y" =>38.4615,
    "Z" => 38.4615
];

10,000 次运行后的结果

Array
(
    [A] => 10.324
    [B] => 59.298
    [C] => 15.902
    [D] => 56.299
    [E] => 21.16
    [F] => 53.621
    [G] => 25.907
    [H] => 50.163
    [I] => 30.932
    [J] => 47.114
    [K] => 35.344
    [L] => 43.175
    [M] => 39.141
    [N] => 39.127
    [O] => 39.346
    [P] => 39.364
    [Q] => 39.501
    [R] => 39.05
    [S] => 39.555
    [T] => 39.239
    [U] => 39.283
    [V] => 39.408
    [W] => 39.317
    [X] => 39.339
    [Y] => 39.569
    [Z] => 39.522
)

【问题讨论】：

The probability that an object An appears in B is Pn. 这很棘手，我相信这不是您想要的。具体来说，如果k=n/2，至少有一半的元素应该有B_i>=1/2。
@amit 我很确定这就是我想要的，但我对我没有正确描述我的目标的可能性持开放态度。 K != n/2 bust 而不是 K < n/2，通常比 n/2 少很多，看看我上面说的示例数字。我也不明白B_i 是什么意思。
在示例中生成“A”的概率是 2.5？在这种情况下不是概率，概率一定在[0,1]范围内
@amit 这些是百分比，所以 2.5% -> .025
好的，现在我正在关注您，是的，您的术语是正确的。感谢您的澄清，将尝试提出答案。

标签： php algorithm probability

【解决方案1】：

我们必须有sum_i P_i = k，否则我们无法成功。

如前所述，问题有点简单，但你可能不喜欢这个答案，理由是它“不够随机”。

Sample a uniform random permutation Perm on the integers [0, n)
Sample X uniformly at random from [0, 1)
For i in Perm
    If X < P_i, then append A_i to B and update X := X + (1 - P_i)
    Else, update X := X - P_i
End

您需要使用定点算法而不是浮点来近似涉及实数的计算。

缺少的条件是分布具有称为“最大熵”的技术属性。像阿米特一样，我想不出一个好的方法来做到这一点。这是一个笨拙的方法。

我解决这个问题的第一个（也是错误的）直觉是将每个A_i 独立地包含在B 中，概率为P_i，然后重试直到B 的长度合适（不会有太多的重试次数，出于您可以询问 math.SE 的原因）。问题是条件作用弄乱了概率。如果P_1 = 1/3 和P_2 = 2/3 和k = 1，那么结果是

{}: probability 2/9
{A_1}: probability 1/9
{A_2}: probability 4/9
{A_1, A_2}: probability 2/9,

条件概率实际上是1/5 对应A_1 和4/5 对应A_2。

相反，我们应该替换新的概率Q_i，以产生适当的条件分布。我不知道Q_i 的封闭形式，所以我建议使用像gradient descent 这样的数值优化算法来找到它们。初始化Q_i = P_i（为什么不呢？）。使用动态规划，对于Q_i 的当前设置，给定带有l 元素的结果，可以找到A_i 是这些元素之一的概率。（我们只关心l = k 条目，但我们需要其他条目才能使循环起作用。）再做一些工作，我们可以得到整个梯度。抱歉，这太粗略了。

在 Python 3 中，使用似乎总是收敛的非线性求解方法（同时将每个 q_i 更新为其略微正确的值并归一化）：

#!/usr/bin/env python3
import collections
import operator
import random


def constrained_sample(qs):
    k = round(sum(qs))
    while True:
        sample = [i for i, q in enumerate(qs) if random.random() < q]
        if len(sample) == k:
            return sample


def size_distribution(qs):
    size_dist = [1]
    for q in qs:
        size_dist.append(0)
        for j in range(len(size_dist) - 1, 0, -1):
            size_dist[j] += size_dist[j - 1] * q
            size_dist[j - 1] *= 1 - q
    assert abs(sum(size_dist) - 1) <= 1e-10
    return size_dist


def size_distribution_without(size_dist, q):
    size_dist = size_dist[:]
    if q >= 0.5:
        for j in range(len(size_dist) - 1, 0, -1):
            size_dist[j] /= q
            size_dist[j - 1] -= size_dist[j] * (1 - q)
        del size_dist[0]
    else:
        for j in range(1, len(size_dist)):
            size_dist[j - 1] /= 1 - q
            size_dist[j] -= size_dist[j - 1] * q
        del size_dist[-1]
    assert abs(sum(size_dist) - 1) <= 1e-10
    return size_dist


def test_size_distribution(qs):
    d = size_distribution(qs)
    for i, q in enumerate(qs):
        d1a = size_distribution_without(d, q)
        d1b = size_distribution(qs[:i] + qs[i + 1 :])
        assert len(d1a) == len(d1b)
        assert max(map(abs, map(operator.sub, d1a, d1b))) <= 1e-10


def normalized(qs, k):
    sum_qs = sum(qs)
    qs = [q * k / sum_qs for q in qs]
    assert abs(sum(qs) / k - 1) <= 1e-10
    return qs


def approximate_qs(ps, reps=100):
    k = round(sum(ps))
    qs = ps[:]
    for j in range(reps):
        size_dist = size_distribution(qs)
        for i, p in enumerate(ps):
            d = size_distribution_without(size_dist, qs[i])
            d.append(0)
            qs[i] = p * d[k] / ((1 - p) * d[k - 1] + p * d[k])
        qs = normalized(qs, k)
    return qs


def test(ps, reps=100000):
    print(ps)
    qs = approximate_qs(ps)
    print(qs)
    counter = collections.Counter()
    for j in range(reps):
        counter.update(constrained_sample(qs))
    test_size_distribution(qs)
    print("p", "Actual", sep="\t")
    for i, p in enumerate(ps):
        print(p, counter[i] / reps, sep="\t")


if __name__ == "__main__":
    test([2 / 3, 1 / 2, 1 / 2, 1 / 3])

【讨论】：

我无法理解这种方法，介意通过它进入聊天室吗？ chat.stackexchange.com/rooms/25165/…
@HurricaneDevelopment 我因为一些愚蠢的原因无法登录，但我现在真的没有时间聊天。也许我稍后会发布一个实现。
@HurricaneDevelopment 发布了一个似乎可行的 Python 实现。
运行这个实现会产生全 0 的实际值。预期的最终产品是什么？
@HurricaneDevelopment 实际分布。你没有在 Python 2 中运行它，是吗？

【解决方案2】：

让我们分析一下。替换：（不是您想要的，但更易于分析）。

给定一个大小为k 的列表L 和元素a_i，a_i 在列表中的概率由您的值p_i 表示。

让我们检查a_i 位于列表中某个索引j 的概率。让我们将该概率表示为q_i,j。请注意，对于列表中的任何索引t，q_i,j = q_i,t - 所以我们可以简单地说q_i_1=q_i_2=...=q_i_k=q_i。

a_i 在列表中任意位置的概率表示为：

1-(1-q_i)^k

但它也是p_i - 所以我们需要解方程

1-(1-q_i)^k = pi
1 - (1-q_i)^k -pi = 0

一种方法是newton-raphson method。

计算每个元素的概率后，检查它是否确实是一个概率空间（总和为 1，所有概率都在 [0,1] 中）。如果不是 - 对于给定的概率和k，它无法完成。

没有替换：这更棘手，因为现在q_i,j != q_i,t（选择不是 i.i.d）。这里的概率计算会更棘手，我目前不确定如何计算它们，我想在创建列表期间需要在运行时完成。

（删除了一个我几乎可以肯定是有偏见的解决方案）。

【讨论】：

是的，所以您对 With Replacements 的逻辑是完美的，我对此进行了测试并且它有效，但正如您所说，它不适用于 With Replacements .我有一种感觉，每一个单品被挑选出来后，所有q_i都必须重新计算。我知道问题的一个条件是sum_i p_i = k 但是这对每个步骤都是必要的吗？例如，在我选择第一个项目后，所有 q_i 将重新计算，然后应该 sum_i p_i still = k 其中 k 现在是 k - 1？
@HurricaneDevelopment 是的，这是必要的。似乎没有计算q_is 的好方法。当然，简单的归一化不是答案。
@DavidEisenstat 我正在尝试复制算法使用替换，但结果似乎与输入百分比不匹配，介意看看吗？它太多不适合评论所以这里是一个pastebin。 pastebin.com/S5jfaJyF
@HurricaneDevelopment 我没有看到您检查q_is 总和为1 的部分。我猜他们没有。
@DavidEisenstat 是的，你是对的，所以 q_i 的总和必须为 1 并且 p_i 的总和也必须为 k?

【解决方案3】：

除非我的数学技能比我认为的列表 A 中的元素在列表 B 中找到的平均几率低很多，否则应该是 10/26 = 0.38。
如果您降低任何对象的此机会，则必须有其他对象具有更高的机会。此外，您来自列表 A 的概率无法计算：它们太低：您无法填写您的列表/您没有足够的元素可供选择。

假设上述内容正确（或足够正确），这意味着在您的列表 A 中，您的平均体重必须是随机选择的平均机会。反过来，这意味着您在列表 a 中的概率总和不等于 100。

除非我完全错了，否则就是……

【讨论】：

您并没有完全错，但不幸的是，我尝试了其中的一些方法，但没有成功。您的 10/26 非常接近，但是我相信因为没有替代品，实际的 calc 是 25_C_9 / 26_C_10，结果是 .3846。无论如何，想法是一样的。我在某处读到所有p_i 的总和必须等于K.26 * .3846 = 10。我最初用所有p_i = .3846 尝试了 10,000 次，它成功了！但是在更改保持约束sum_i p_i = k 的数字后，它停止工作。请查看我的编辑以查看结果。