pandas 过滤异常值 - 如何使基于中值的 Hampel 函数更快？
声明：本页面是StackOverFlow热门问题的中英对照翻译，遵循CC BY-SA 4.0协议，如果您需要使用它，必须同样遵循CC BY-SA许可，注明原文地址和作者信息，同时你必须将它归于原作者（不是我）：StackOverFlow
原文地址: http://stackoverflow.com/questions/46819260/
Warning: these are provided under cc-by-sa 4.0 license. You are free to use/share it, But you must attribute it to the original authors (not me):
StackOverFlow
Filtering Outliers - how to make median-based Hampel Function faster?
提问by EHB
I need to use a Hampel filter on my data, stripping outliers.
我需要对我的数据使用 Hampel 过滤器，去除异常值。
I haven't been able to find an existing one in Python; only in Matlab and R.
我一直无法在 Python 中找到现有的；仅在 Matlab 和 R 中。
[Matlab function description][1]
【Matlab函数说明】[1]
[Stats Exchange discussion of Matlab Hampel function][2]
[Matlab Hampel函数的Stats Exchange讨论][2]
[R pracma package vignette; contains hampel function][3]
[R pracma 包小插图；包含hampel功能][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.
我编写了以下函数，根据 R pracma 包中的函数对其进行建模；但是，它远比 Matlab 版本慢得多。这并不理想；将不胜感激有关如何加快速度的意见。
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:
“pracma”包中的 R 实现，我将其用作模型：
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:
任何使函数更高效的帮助，或指向现有 Python 模块中现有实现的指针，将不胜感激。下面的示例数据；Jupyter 中的 %%timeit cell magic 表示当前需要 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
[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
采纳答案by EHB
A Pandas solution is several orders of magnitude faster:
Pandas 解决方案要快几个数量级：
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.
计时这给出了 11 ms vs 15 秒；巨大的改进。
I found a solution for a similar filter in this post.
回答by Eduardo Osorio
Solution by @EHB above is helpful, but it is incorrect. Specifically, the rolling median calculated in median_abs_deviationis of difference, which itself is the difference between each data point and the rolling median calculated in rolling_median, but it should be the median of differences between the data in the rolling window and the median over the window. I took the code above and modified it:
上面@EHB 的解决方案很有帮助，但它是不正确的。具体来说，median_abs_deviation中计算出的滚动中位数有差异，它本身就是每个数据点与rolling_median中计算出的滚动中位数的差值，但应该是滚动窗口中的数据与窗口上的中位数差异的中位数. 我把上面的代码修改了一下：
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)
