将半高斯滤波器应用于python中的分箱时间序列数据

Question

我正在对一些时间序列数据进行分箱，我需要对分箱数据应用半正态过滤器。我怎样才能在python中做到这一点？我在下面提供了一个玩具示例。我需要使用标准为 0.25（或其他任何值）的半高斯滤波器对 Xbinned 进行平滑处理。我很确定半高斯应该面向正向时间方向。

import numpy as np

X = np.random.randint(2, size=100) #example random process

bin_size =  5

Xbinned = []

for i in range(0, len(X)+1, bin_size):
    Xbinned.append(sum(X[i:i+(bin_size-1)])/bin_size)

Answer 1

如何实现半高斯滤波

Scipy 有一个名为scipy.ndimage.gaussian_filter() 的函数。它几乎实现了我们在这里想要的。不幸的是，没有选择使用半高斯而不是高斯。但是 scipy 是开源的，所以我们可以把source code 修改为半高斯。

我使用了这个源代码，并删除了这个特殊情况不需要的所有部分。最后，我得到了这个：

import scipy.ndimage

def halfgaussian_kernel1d(sigma, radius):
    """
    Computes a 1-D Half-Gaussian convolution kernel.
    """
    sigma2 = sigma * sigma
    x = np.arange(0, radius+1)
    phi_x = np.exp(-0.5 / sigma2 * x ** 2)
    phi_x = phi_x / phi_x.sum()

    return phi_x

def halfgaussian_filter1d(input, sigma, axis=-1, output=None,
                      mode="constant", cval=0.0, truncate=4.0):
    """
    Convolves a 1-D Half-Gaussian convolution kernel.
    """
    sd = float(sigma)
    # make the radius of the filter equal to truncate standard deviations
    lw = int(truncate * sd + 0.5)
    weights = halfgaussian_kernel1d(sigma, lw)
    origin = -lw // 2
    return scipy.ndimage.convolve1d(input, weights, axis, output, mode, cval, origin)

这是如何工作的简短摘要：

1. 首先，它生成一个卷积核。它使用公式e^(-1/2 * (x/sigma)^2) 生成高斯分布。它会一直持续到距离中心 4 个标准差为止。
2. 接下来，它将内核与您的信号进行卷积。它将内核调整为从当前时间步开始，而不是以当前时间步为中心。

在你的信号上尝试这个，我得到这样的结果：

array([0.59979879, 0.6       , 0.40006707, 0.59993293, 0.79993293,
       0.40013414, 0.20006707, 0.59986586, 0.40006707, 0.4       ,
       0.99979879, 0.00033535, 0.59979879, 0.40006707, 0.00013414,
       0.59979879, 0.20013414, 0.00006707, 0.19993293, 0.59986586])

标准差的选择

如果您选择 0.25 的标准差，这对您的信号几乎没有影响。以下是它使用的卷积权重：[0.99966465 0.00033535]。换句话说，这对信号的影响不到 0.1%。

我建议使用更大的 sigma 值。

因一个错误而关闭

另外，我想在这里指出一个错误：

for i in range(0, len(X)+1, bin_size):
    Xbinned.append(sum(X[i:i+(bin_size-1)])/bin_size)

Numpy 范围不包含在内，因此 i 到 i+(bin_size-1) 的范围实际上捕获 4 个元素，而不是 5 个。

要解决此问题，您可以将其更改为：

for i in range(0, len(X), bin_size):
    Xbinned.append(X[i:i+bin_size].mean())

（另外，我修复了循环规范中的一个错误，并使用了一个 numpy 快捷方式来查找平均值。）