優化循環變點測試

Question

我想在Python中編寫一個簡單的變化點查找器。下面，函數loglike（xs）返回iid正常樣本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 

如果你真的想擠進毫秒，你可以使用瓶頸的nanstd函數來代替，因爲它更快。如果你想報廢幾微秒，你可以用math.log代替np.log，因爲你只用一個數字操作，如果xs是一個數組，你可以用xs.std()代替。但在走下這條路之前，我建議你使用這個版本，並對結果進行分析以查看時間花在哪裏。

編輯

如果配置文件loglike python -m cProfile -o output yourprogram.py; runsnake output，你會看到，大部分（約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) 

我又sqeeze一個aditional的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，或者如果您不打算更改dtype，則只使用func.nanstd_1d_float64_axisNone。

我認爲這與Python中的速度一樣快。儘管如此，一半的時間都花在了數字操作上，也許Cython可以優化這一點，但隨後調用Python和Python會增加一些開銷，可以彌補這一點。

Answer 2

我會做一些基本的改變：目前，您在每次迭代時重新計算每個分區的平方和，平均值和計數，這使得這個「交叉驗證樣式」算法是二次的。

你可以做的是利用半羣結構，並在更改分區時以聯機方式計算每個元素的平方值，計數和平均值 - 基本上將np.sum中的隱式循環融合在一起。然後，您採取-0.5*n*np.log(2*np.pi*sigSq) - 0.5*n並根據更新後的值n，mean和sigSq（您將需要從平方值之和計算stddev）來計算。

與np.std相比，這將漸近地加快代碼的速度並節省一些函數調用的數值穩定性的潛在代價。你可能需要卡漢總結。

如果您不需要實際的對數似然性，您可以專注於僅最大化proflike，然後在循環外部計算lr，從而節省幾個乘數 - （您可以在函數中摺疊2 *和0.5 *一起如果你做一些基本的代數無論）

這個技術是HLearn的背後，它的交叉驗證的高性能。

編輯：一些代碼可能比其他代碼快得多。這裏最大的代價是迭代開銷和np.log調用。可能有一些fencepost錯誤，但這裏是要點：

rcnt = n 
rsum = np.sum(arr) 
rssq = np.sum(arr**2) 
lcnt = 0 
lsum = 0 
lssq = 0 

maxlike = -1e9 # or -Inf ideally 
cp = -1 

for i in arr: # arr is length n 
    lcnt += 1 
    lsum += i 
    lssq += i*i 
    lmean = lsum/lcnt 
    lvar = ((lssq - lmean**2)/(lcnt)) 
    loglike_l = lcnt*np.log(2*np.pi*lvar) - lcnt 

    rcnt -= 1 
    rsum -= i 
    rssq -= i*i 
    rmean = rsum/rcnt 
    rvar = ((rssq - rmean**2)/(rcnt)) 
    loglike_r = rcnt*np.log(2*np.pi*rvar) - rcnt 

    loglike_total = loglike_l + loglike_r 

    if maxlike < loglike_total: 
     cp = lcnt 
     maxlike = loglike_total