algorithm - 如何每次不放样重新称量(条件抽样)？

Question

考虑一个包含 N 行权重的数据集。这是基本算法:

1. 将权重归一化，使它们的总和为 1。
2. 将权重备份到另一列以记录样本概率
3. 随机选择 1 行(不放回)，给定样本概率，并将其添加到样本数据集中
4. 从原始数据集中移除抽取的权重，并通过归一化剩余行的权重重新计算样本概率
5. 重复步骤3和4，直到样本中的权重之和达到或超过阈值(假设0.6)

这是一个玩具示例:

import pandas as pd
import numpy as np

def sampler(n):
    df = pd.DataFrame(np.random.rand(n), columns=['weight'])
    df['weight'] = df['weight']/df['weight'].sum()
    df['samp_prob'] = df['weight']

    samps = pd.DataFrame(columns=['weight'])

    while True:
        choice = np.random.choice(df.index, 1, replace=False, p=df['samp_prob'])[0]
        samps.loc[choice, 'weight'] = df.loc[choice, 'weight']
        df.drop(choice, axis=0, inplace=True)
        df['samp_prob'] = df['weight']/df['weight'].sum()
        if samps['weight'].sum() >= 0.6:
            break
    return samps

玩具示例的问题是随着 n 大小的增加，运行时间呈指数增长:

Answer 1

出发进近

一些观察:

• 导致创建新数据帧的每次迭代丢弃的行对性能没有帮助。
• 看起来不太容易矢量化，但应该很容易处理底层数组数据以提高性能。这个想法是使用掩码并避免重新创建数据帧或数组。一开始，我们将使用两列数组，对应于名为:'weights' 的列和 'samp_prob' .

所以，考虑到这些，开始的方法应该是这样的——

def sampler2(n):
    a = np.random.rand(n,2)
    a[:,0] /= a[:,0].sum()
    a[:,1] = a[:,0]
    N = len(a)

    idx = np.arange(N)
    mask = np.ones(N,dtype=bool)
    while True:
        choice = np.random.choice(idx[mask], 1, replace=False, p=a[mask,1])[0]
        mask[choice] = 0
        a_masked = a[mask,0]
        a[mask,1] = a_masked/a_masked.sum()
        if a[~mask,0].sum() >= 0.6:
            break
    out = a[~mask,0]
    return out

改进 #1

后来的观察表明数组的第一列在迭代中没有改变。因此，我们可以通过预先计算总和然后在每次迭代中优化第一列的屏蔽总和，a[~mask,0].sum()将只是总和减去 a_masked.sum() . Thsi 引导我们进行第一项改进，如下所列 -

def sampler3(n):
    a = np.random.rand(n,2)
    a[:,0] /= a[:,0].sum()
    a[:,1] = a[:,0]
    N = len(a)

    idx = np.arange(N)
    mask = np.ones(N,dtype=bool)
    a0_sum = a[:,0].sum()
    while True:
        choice = np.random.choice(idx[mask], 1, replace=False, p=a[mask,1])[0]
        mask[choice] = 0
        a_masked = a[mask,0]
        a_masked_sum = a_masked.sum()
        a[mask,1] = a_masked/a_masked_sum
        if a0_sum - a_masked_sum >= 0.6:
            break
    out = a[~mask,0]
    return out

改进 #2

现在，假设第一列在迭代之间没有变化，则可以通过使用两个单独的数组来改进对二维数组列的切片和 mask 。这给了我们一个修改后的版本，就像这样 -

def sampler4(n):
    a = np.random.rand(n)
    a /= a.sum()
    b = a.copy()
    N = len(a)

    idx = np.arange(N)
    mask = np.ones(N,dtype=bool)
    a_sum = a.sum()
    while True:
        choice = np.random.choice(idx[mask], 1, replace=False, p=b[mask])[0]
        mask[choice] = 0
        a_masked = a[mask]
        a_masked_sum = a_masked.sum()
        b[mask] = a_masked/a_masked_sum
        if a_sum - a_masked_sum >= 0.6:
            break
    out = a[~mask]
    return out

运行时测试 -

In [250]: n = 1000

In [251]: %timeit sampler(n) # original app
     ...: %timeit sampler2(n)
     ...: %timeit sampler3(n)
     ...: %timeit sampler4(n)
1 loop, best of 3: 655 ms per loop
10 loops, best of 3: 50 ms per loop
10 loops, best of 3: 44.9 ms per loop
10 loops, best of 3: 38.4 ms per loop

In [252]: n = 2000

In [253]: %timeit sampler(n) # original app
     ...: %timeit sampler2(n)
     ...: %timeit sampler3(n)
     ...: %timeit sampler4(n)
1 loop, best of 3: 1.32 s per loop
10 loops, best of 3: 134 ms per loop
10 loops, best of 3: 119 ms per loop
10 loops, best of 3: 100 ms per loop

因此，我们得到 17x+ 和 13x+ 最终版本比 n=1000 的原始方法提速和 n=2000尺寸!