过滤异常值如何使基于媒体的Hampel函数更快？

Question

我需要对我的数据使用Hampel过滤器，剔除异常值。在

我在Python中还没有找到一个现有的；只有在Matlab和R中

【Matlab函数描述】【1】

【Matlab-Hampel函数的统计交换讨论】【2】

[R pracma package vignette；包含hampel函数][3]

我已经编写了下面的函数，用pracma包中的函数对它进行建模；但是，它比Matlab版本慢得多。这不太理想；希望能提供加快速度的建议。在

函数如下所示-

def hampel(x,k, t0=3):
    '''adapted from hampel function in R package pracma
    x= 1-d numpy array of numbers to be filtered
    k= number of items in window/2 (# forward and backward wanted to capture in median filter)
    t0= number of standard deviations to use; 3 is default
    '''
    n = len(x)
    y = x #y is the corrected series
    L = 1.4826
    for i in range((k + 1),(n - k)):
        if np.isnan(x[(i - k):(i + k+1)]).all():
            continue
        x0 = np.nanmedian(x[(i - k):(i + k+1)])
        S0 = L * np.nanmedian(np.abs(x[(i - k):(i + k+1)] - x0))
        if (np.abs(x[i] - x0) > t0 * S0):
            y[i] = x0
    return(y)

“pracma”包中的R实现，我将其用作模型：

^{pr2}$

如果有助于提高函数的效率，或者在现有Python模块中提供指向现有实现的指针，我们将不胜感激。下面的示例数据；%%timeit cell magic在Jupyter中表示当前运行需要15秒：

vals=np.random.randn(250000)
vals[3000]=100
vals[200]=-9000
vals[-300]=8922273
%%timeit
hampel(vals, k=6)

[1]：https://www.mathworks.com/help/signal/ref/hampel.html[2]：https://dsp.stackexchange.com/questions/26552/what-is-a-hampel-filter-and-how-does-it-work[3]：https://cran.r-project.org/web/packages/pracma/pracma.pdf

Answer 1

上面@EHB的解决方案是有帮助的，但它是不正确的。具体地说，在median_abs_deviation中计算的滚动中值属于difference，它本身就是每个数据点与rolling_mean中计算的滚动中值之间的差值，但它应该是滚动窗口中的数据与窗口上的中值之间差异的中值。我把上面的代码改了：

def hampel(vals_orig, k=7, t0=3):
    '''
    vals: pandas series of values from which to remove outliers
    k: size of window (including the sample; 7 is equal to 3 on either side of value)
    '''

    #Make copy so original not edited
    vals = vals_orig.copy()

    #Hampel Filter
    L = 1.4826
    rolling_median = vals.rolling(window=k, center=True).median()
    MAD = lambda x: np.median(np.abs(x - np.median(x)))
    rolling_MAD = vals.rolling(window=k, center=True).apply(MAD)
    threshold = t0 * L * rolling_MAD
    difference = np.abs(vals - rolling_median)

    '''
    Perhaps a condition should be added here in the case that the threshold value
    is 0.0; maybe do not mark as outlier. MAD may be 0.0 without the original values
    being equal. See differences between MAD vs SDV.
    '''

    outlier_idx = difference > threshold
    vals[outlier_idx] = np.nan
    return(vals)

Answer 2

熊猫解决方案的速度快了几个数量级：

def hampel(vals_orig, k=7, t0=3):
    '''
    vals: pandas series of values from which to remove outliers
    k: size of window (including the sample; 7 is equal to 3 on either side of value)
    '''
    #Make copy so original not edited
    vals=vals_orig.copy()    
    #Hampel Filter
    L= 1.4826
    rolling_median=vals.rolling(k).median()
    difference=np.abs(rolling_median-vals)
    median_abs_deviation=difference.rolling(k).median()
    threshold= t0 *L * median_abs_deviation
    outlier_idx=difference>threshold
    vals[outlier_idx]=np.nan
    return(vals)

计时这给予11毫秒对15秒；巨大的改善。在

我在this post.中找到了一个类似过滤器的解决方案