我们必须拥有
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_iEnd
您需要使用定点算术而不是浮点算术来逼近涉及实数的计算。
缺少的条件是该分布具有称为“最大熵”的技术属性。像阿米特（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 python3import collectionsimport operatorimport randomdef 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 sampledef 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_distdef 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_distdef 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-10def 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 qsdef 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 qsdef 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])
