从大小为N的（概率加权）集合中生成随机大小为k的子集

Question

这个问题来源于一个音乐训练游戏，我需要从12个音高类别中随机选择一个3个音符的和弦，但某些音符出现的概率比其他音符高（这样用户可以多练习那些较弱的音符）。

我原以为这个问题很简单：把每个音符的权重看作一段线段，把所有线段一个接一个地放在一起，形成一条长线段，然后在这条长线段上随机选择一个点，记录这个点落在哪个权重上，重复这个过程直到得到k个结果。

下面的Python代码演示了这种方法并没有产生正确的结果：

# Choose k items from a set of weights
# return set of winning indices
def Choose(W,k):
    import random

    cumulative = [ sum(W[:i+1]) for i in xrange(len(W)) ]
    totalWeight = cumulative[-1]

    winners = set()
    while len(winners) < k:
        rnd = random.uniform(0.0, totalWeight)

        # Returns first element of cumulative that is >= rnd
        w = next( i for i in xrange(len(cumulative)) if cumulative[i] >= rnd )
        winners.add( w )

    return winners

def Test(N):
    x = [ list(Choose( [5,3,2], 2 )) for i in xrange(int(N/2))]
    y = sum(x, [])
    z = [y.count(i) for i in (0,1,2) ]

    print z

for i in range(10):
    Test(10000)

我从3个权重[5,3,2]中生成了5000个随机组合 输出记录了每个权重出现的次数 应该是5000,3000,2000

为了确保结果的可靠性，我进行了10次实验：

python test.py 
[4173, 3331, 2496]
[4180, 3367, 2453]
[4193, 3393, 2414]
[4228, 3375, 2397]
[4207, 3388, 2405]
[4217, 3377, 2406]
[4173, 3438, 2389]
[4172, 3378, 2450]
[4174, 3371, 2455]
[4208, 3322, 2470]

结果大约是4200, 3300, 2400 并不是5000, 3000, 2000

有没有简单的方法来理解为什么这样不行呢？

有没有什么方法可以转换这些权重，比如'weight[i] -> ln(weight[i])'之类的，这样可以得到正确的结果吗？

如何才能得到正确的结果？（我更关心代码的清晰度，而不是最佳效率）

Answer 1

使用 numpy.random.choice 这个函数，并且要用到它的 p 参数：

np.random.choice(3, size=1000, p=[0.5, 0.3, 0.2])

现在再试一次，看看你会得到什么结果。

Answer 2

不放回抽样并且带有权重是个复杂的问题。

首先，想想你的直觉解决方案。你生成了5000对数据，并且你希望这5000对中有5000对包含一个1。这意味着每一对都必须包含一个1。我怀疑这并不是你想要的结果。为了得到你期望的分布，你可以先选择1，然后以0.6的概率选择2，以0.4的概率选择3。

为了实现我猜测你想要的效果，你应该做类似条件泊松抽样的操作。不过，我不知道有没有Python模块可以做到这一点，虽然几乎肯定有。R语言中的'sampling'包可以做到这一点。我在网上没有找到简单易懂的介绍。

从实际的角度来看，你可以继续你现在的做法，并调整权重，使得概率接近你想要的结果。对于你想做的事情，精确的概率似乎并不是必要的。

如果你想要一个简单的方法（虽然效率不高）来实现你的目标：

1) 先把权重归一化，使得所有权重的总和等于你想要的样本大小。以你的例子来说，0.5 + 0.3 + 0.2 = 2，所以归一化后的权重就是[1., .6, .4]。

2) 让p_i表示第i个权重，把它当作概率（它们都必须小于或等于1，否则问题就无法解决）。通过以概率p_i选择第i个元素来抽样。

3) 如果抽取的样本大小正确，就输出它；否则就重新抽取。

这里有个简单的代码示例：

import random
def sample(weights, sample_size):
    w = float(sum(weights))
    normweights = [x * sample_size / w for x in weights]
    samp = [random.random() < pi for pi in normweights]
    while sum(samp) != sample_size:
        samp = [random.random() < pi for pi in normweights]
    return [i for i,b in enumerate(samp) if b]

print(sample([.5,.3,.2], 2))

编辑：
好的，上面的算法其实不太靠谱。我会尽量记住怎么正确地做。