使用 itertools 更快地循环

Question

我有一个函数

def getSamples():
    p = lambda x : mlab.normpdf(x,3,2) + mlab.normpdf(x,-5,1)
    q = lambda x : mlab.normpdf(x,5,14)
    k=30
    goodSamples = []
    rightCount = 0
    totalCount = 0
    while(rightCount < 100000):
        z0 = np.random.normal(5, 14)
        u0 = np.random.uniform(0,k*q(z0))
        if(p(z0) > u0):
            goodSamples.append(z0)
            rightCount += 1
        totalCount += 1
    return np.array(goodSamples)

我生成 100000 个样本的实现需要很长时间。我怎样才能用itertools 或类似的东西让它更快？

Answer 1

我想说，让这段代码更快的秘诀不在于改变循环语法。这里有几点：

1. np.random.normal 有一个额外的参数 size 可以让您一次获取多个值。我建议使用一组 1E09 元素，然后检查你的条件，看看有多少是好的。然后，您可以估计发生这种情况的可能性。
2. 要创建统一样本，为什么不使用sympy 对pdf 进行符号评估？ （我不知道这是否更快，但可能是因为您已经知道均值和方差。）
3. 同样，p 可以使用符号函数吗？

Answer 2

一般来说，性能问题是由“错误的方式”做事引起的。 Numpy 在使用时可以非常快，因为它被设计为使用，即通过利用其向量处理，将这些向量化操作传递给编译代码。来自其他编程语言/方法的两种不良做法是

1. 循环：每当您认为需要循环停止并思考时。大多数时候你不想要，事实上甚至不想要一个。在没有循环的情况下编写和运行代码要快得多。
2. 内存分配：只要知道对象的大小，就为它预先分配空间。与其他替代方案相比，内存增长速度非常慢，尤其是在 Python 列表中。

在这种情况下，很容易获得（大约）两个数量级的加速；权衡是更多的内存使用。

以下是一些代表代码，不代表盲目使用。我甚至没有验证它会产生正确的结果。它或多或少是您日常工作的直接翻译。您似乎正在使用拒绝方法从概率分布中抽取随机数。对于您的概率分布，可能有更有效的算法来执行此操作。

def getSamples2() :
    p = lambda x : mlab.normpdf(x,3,2) + mlab.normpdf(x,-5,1)
    q = lambda x : mlab.normpdf(x,5,14)
    k=30
    N = 100000 # Total number of samples we want
    Ngood = 0 # Current number of good samples
    goodSamples = np.zeros(N) # Storage for the good samples
    while Ngood < N : # Unfortunately a loop, ....
       z0 = np.random.normal(5, 14, size=N)
       u0 = np.random.uniform(size=N)*k*q(z0)
       ind, = np.where(p(z0) > u0)
       n = min(len(ind), N-Ngood)
       goodSamples[Ngood:Ngood+n] = z0[ind[:n]]
       Ngood += n
    return goodSamples

这会以块的形式生成随机数并保存好的随机数。我没有尝试优化块大小（这里我只使用N，我们想要的总数，原则上这可以/应该不同，甚至可以根据我们剩下的要生成的数量进行调整）。不幸的是，这仍然使用循环，但现在它将运行“数十”次而不是 100,000 次。这也使用了where 函数和数组切片；这些都是很好的通用工具。

在我的机器上使用%timeit 进行的一项测试中我发现

In [27]: %timeit getSamples() # Original routine
1 loops, best of 3: 49.3 s per loop
In [28]: %timeit getSamples2()
1 loops, best of 3: 505 ms per loop

Answer 3

这有点 itertools 的“魔法”，但我不确定它是否有帮助。准备一个 numpy 数组（使用零）并在不创建 python 自动增长列表的情况下填充它可能会更好。这是迭代工具和零准备。 （对于未经测试的代码，请提前原谅）

from itertools import count, ifilter, imap, takewhile
import operator


def getSamples():
    p = lambda x : mlab.normpdf(x, 3, 2) + mlab.normpdf(x, -5, 1)
    q = lambda x : mlab.normpdf(x, 5, 14)

    k = 30
    n = 100000
    samples_iter = imap(
        operator.itemgetter(1),
        takewhile(
            lambda i, s: i < n,
            enumerate(
                ifilter(lambda z: p(z) > np.random.uniform(0,k*q(z)),
                        (np.random.normal(5, 14) for _ in count()))
    )))

    goodSamples = numpy.zeros(n)
    # set values from iterator, probably there is a better way for that
    for i, sample in enumerate(samples_iter):
        goodSamples[i] = sample
    return goodSamples