从可变权重随机生成组合
Randomly Generating Combinations From Variable Weights
非常重要的编辑： 所有 Ai 都是 独一无二.
问题
我有一个列表 A of n unique 对象。每个对象 Ai 都有一个可变百分比 Pi .
我想创建一个算法来生成一个包含 k 个对象的新列表 B (k < n/2 在大多数情况下 k 明显小于 n/2。例如 n=231 , k =21)。列表 B 不应有重复项，并将填充来自列表 A 的对象，但具有以下限制：
The probability that an object Ai appears in B is Pi.
我试过的
（这些片段在 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[=77 中选择 Z 的机会=]n 是 20%。这不是我想要的。我希望 Z 在列表 B 中被选中的几率为 20%。
这可能吗？怎么做到的？
编辑 1
我尝试简单地计算所有 Pi = k 的总和，如果所有 Pi是相等的，但是修改了它们的值后，就开始越来越错了。
初始概率
$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
)
让我们分析一下。
替换：（不是你想要的，但更容易分析）。
给定大小为 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。
计算每个元素的概率后，检查它是否确实是概率space（总和为1，所有概率都在[0,1]）。如果不是 - 对于给定的概率和 k.
则无法完成
无替换：这比较棘手，因为现在 q_i,j != q_i,t（选择不是 i.i.d）。这里的概率计算会更加棘手，目前我不确定如何计算它们，需要在 运行 时间内完成，我想是在创建列表的过程中。
（删除了一个我几乎可以肯定是有偏见的解决方案）。
除非我的数学技能比我认为的要弱得多，否则在列表 B 中找到示例中列表 A 中的元素的平均机会应该是 10/26 = 0.38。
如果你降低任何物体的这个机会，一定有其他人有更高的机会。
此外，列表 A 中的概率无法计算：它们太低：您无法填写列表/您没有足够的元素可供选择。
假设以上是正确的（或足够正确），这意味着在您的列表 A 中，您的平均体重必须是随机选择的平均机会。反过来，这意味着您在列表 a 中的概率总和不等于 100。
除非我完全错了，那就是...
必须有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,
条件概率实际上是 A_1 的 1/5 和 A_2 的 4/5。
相反，我们应该替换产生适当条件分布的新概率 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])
非常重要的编辑： 所有 Ai 都是 独一无二.
问题
我有一个列表 A of n unique 对象。每个对象 Ai 都有一个可变百分比 Pi .
我想创建一个算法来生成一个包含 k 个对象的新列表 B (k < n/2 在大多数情况下 k 明显小于 n/2。例如 n=231 , k =21)。列表 B 不应有重复项，并将填充来自列表 A 的对象，但具有以下限制：
The probability that an object Ai appears in B is Pi.
我试过的
（这些片段在 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[=77 中选择 Z 的机会=]n 是 20%。这不是我想要的。我希望 Z 在列表 B 中被选中的几率为 20%。
这可能吗？怎么做到的？
编辑 1
我尝试简单地计算所有 Pi = k 的总和，如果所有 Pi是相等的，但是修改了它们的值后，就开始越来越错了。
初始概率
$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
)
让我们分析一下。 替换：（不是你想要的，但更容易分析）。
给定大小为 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。
计算每个元素的概率后，检查它是否确实是概率space（总和为1，所有概率都在[0,1]）。如果不是 - 对于给定的概率和 k.
无替换：这比较棘手，因为现在 q_i,j != q_i,t（选择不是 i.i.d）。这里的概率计算会更加棘手，目前我不确定如何计算它们，需要在 运行 时间内完成，我想是在创建列表的过程中。
（删除了一个我几乎可以肯定是有偏见的解决方案）。
除非我的数学技能比我认为的要弱得多，否则在列表 B 中找到示例中列表 A 中的元素的平均机会应该是 10/26 = 0.38。
如果你降低任何物体的这个机会，一定有其他人有更高的机会。
此外，列表 A 中的概率无法计算：它们太低：您无法填写列表/您没有足够的元素可供选择。
假设以上是正确的（或足够正确），这意味着在您的列表 A 中，您的平均体重必须是随机选择的平均机会。反过来，这意味着您在列表 a 中的概率总和不等于 100。
除非我完全错了，那就是...
必须有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,
条件概率实际上是 A_1 的 1/5 和 A_2 的 4/5。
相反，我们应该替换产生适当条件分布的新概率 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])