[FIXED] Filtering Outliers - how to make median-based Hampel Function faster?

Issue

I need to use a Hampel filter on my data, stripping outliers.

I haven't been able to find an existing one in Python; only in Matlab and R.

[Matlab function description][1]

[Stats Exchange discussion of Matlab Hampel function][2]

[R pracma package vignette; contains hampel function][3]

I've written the following function, modeling it off the function in the R pracma package; however, it is far far slower than the Matlab version. This is not ideal; would appreciate input on how to speed it up.

The function is shown below-

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)

The R implementation in "pracma" package, which I am using as a model:

function (x, k, t0 = 3) 
{
    n <- length(x)
    y <- x
    ind <- c()
    L <- 1.4826
    for (i in (k + 1):(n - k)) {
        x0 <- median(x[(i - k):(i + k)])
        S0 <- L * median(abs(x[(i - k):(i + k)] - x0))
        if (abs(x[i] - x0) > t0 * S0) {
            y[i] <- x0
            ind <- c(ind, i)
        }
    }
    list(y = y, ind = ind)
}

Any help in making function more efficient, or a pointer to an existing implementation in an existing Python module would be much appreciated. Example data below; %%timeit cell magic in Jupyter indicates it currently takes 15 seconds to run:

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

Solution

A Pandas solution is several orders of magnitude faster:

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)

Timing this gives 11 ms vs 15 seconds; vast improvement.

I found a solution for a similar filter in this post.

Answered By - EHB