python - 优化变点测试循环

Question

我正在尝试用 Python 编写一个简单的变化点查找器。下面，函数 loglike(xs) 返回独立同分布正态样本 xs 的最大化对数似然。函数 most_probable_cp(xs) 循环遍历 xs 中间 ~75% 的每个点，并使用似然比找到 xs 中最可能的变化点。

我正在使用二元分割，并且正在自举以获取似然比的临界值，因此我需要调用 most_probable_cp() 数千次。有什么办法可以加快速度吗？ Cython 会有帮助吗？我从未使用过它。

import numpy as np

def loglike(xs):
    n = len(xs)
    mean = np.sum(xs)/n
    sigSq = np.sum((xs - mean)**2)/n
    return -0.5*n*np.log(2*np.pi*sigSq) - 0.5*n


def most_probable_cp(xs, left=None, right=None):
    """
    Finds the most probable changepoint location and corresponding likelihood for xs[left:right]
    """
    if left is None:
        left = 0

    if right is None:
        right = len(xs)

    OFFSETPCT = 0.125
    MINNOBS = 12

    ys = xs[left:right]
    offset = min(int(len(ys)*OFFSETPCT), MINNOBS)
    tLeft, tRight = left + offset, right - offset
    if tRight <= tLeft:
        raise ValueError("left and right are too close together.")

    maxLike = -1e9
    cp = None
    dataLike = loglike(ys)
    # Bottleneck is below.
    for t in xrange(tLeft, tRight):
        profLike = loglike(xs[left:t]) + loglike(xs[t:right])
        lr = 2*(profLike - dataLike)
        if lr > maxLike:
            cp = t
            maxLike = lr

    return cp, maxLike

Answer 1

首先，使用 Numpy 的标准差实现。这不仅会更快，而且会更稳定。

def loglike(xs):
    n = len(xs)
    return -0.5 * n * np.log(2 * np.pi * np.std(xs)) - 0.5 * n

如果你真的想压缩毫秒数，你可以使用 bottleneck 的 nanstd 函数，因为它更快。如果你想舍弃微秒，你可以用 math.log 替换 np.log，因为你只对一个数字进行操作，如果 xs 是一个数组，您可以使用 xs.std() 代替。但在走这条路之前，我建议您使用这个版本，并对结果进行剖析以查看时间花在了哪里。

编辑

如果你像那样分析日志，python -m cProfile -o output yourprogram.py; runsnake 输出，您会看到大部分(大约 80%)的时间都花在了计算 np.std 上。这是我们的第一个目标。正如我之前所说，最好的调用是使用 bottleneck.nanstd。

import bottleneck as bn

def loglike(xs):
    n = len(xs)
    return -0.5 * n * np.log(2 * np.pi * bn.nanstd(xs)) - 0.5 * n

在我的基准测试中，它实现了 8 倍的加速，而且只有 30% 的时间。 len 是 5%，所以没有必要进一步研究它。将 np.log 和 np.pi 替换为它们的数学对应项，并取公因数，我可以再次将时间减半。

return -0.5 * n * (math.log(2 * math.pi * bn.nanstd(xs)) - 1)

我还可以压缩额外的 10% 来稍微损害可读性:

factor = math.log(2*math.pi)

def loglike(xs):
    n = len(xs)
    return -0.5 * n * (factor + math.log(bn.nanstd(xs)) - 1)

编辑2

如果你想真正推送它，你可以将 bn.nanstd 替换为专门的功能。在循环之前，定义 std, _ = bn.func.nansum_selector(xs, axis=0) 并使用它代替 bn.nanstd，或者只是 func.nanstd_1d_float64_axisNone 如果你不会改变数据类型。

而且我认为这与在 Python 中获得的速度一样快。尽管如此，仍有一半时间花在了数字运算上，也许 Cython 能够对此进行优化，但随后调用和调用 Python 会增加开销，可以弥补这一点。