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

Question

我正在分箱一些时间序列数据，我需要对分箱数据应用半正态过滤器。我如何在 python 中执行此操作？我在下面提供了一个玩具示例。我需要使用 std 为 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

如何实现half-gaussian过滤

Scipy 有一个名为 scipy.ndimage.gaussian_filter(). It nearly implements what we want here. Unfortunately, there's no option to use a half-gaussian instead of a gaussian. However, scipy is open-source, so we can just take the source code 的函数并将其修改为 half-gaussian.

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

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)

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

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

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

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 值。

因一个错误而关闭

此外，我想指出这里的 off-by-one 错误：

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())

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

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

Applying a half-gaussian filter to binned time series data in python

python

filtering

time-series

如何实现half-gaussian过滤

标准差的选择

因一个错误而关闭