我们必须拥有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
您需要使用定点算术而不是浮点算术来逼近涉及实数的计算。
缺少的条件是该分布具有称为“最大熵”的技术属性。像阿米特(Amit)一样,我想不出一个好方法。这是一种笨拙的方式。
解决这个问题的第一个(也是错误的)本能是将每个事件独立地包含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
,所以我建议用像梯度下降这样的数值优化算法找到它们。初始化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])