不规则间隔点的高斯和滤波器

Question

我有一组点，分别用 (x,y) 表示，这两个是向量。

from pylab import *
x = sorted(random(30))
y = random(30)
plot(x,y, 'o-')

现在我想用高斯平滑这些数据，并且只在 x 轴上的某些特定（均匀间隔的）点上进行评估。比如说：

x_eval = linspace(0,1,11)

我听说这个方法叫做“高斯和滤波器”，但到目前为止我还没有找到在 numpy/scipy 中的实现，尽管乍一看这似乎是个标准问题。

因为 x 值不是均匀分布的，所以我不能使用 scipy.ndimage.gaussian_filter1d。

通常，这种平滑处理是通过傅里叶空间进行的，然后与核相乘，但我不太确定这是否适用于不规则间隔的数据。

谢谢大家的想法。

Answer 1

根据@Jaime的回答，我写了一个函数，里面有一些额外的说明，并且可以丢弃那些离数据点很远的估算值。

我觉得可以通过自助法（bootstrapping）来获取这个估算值的置信区间，但我还没有尝试过。

def gaussian_sum_smooth(xdata, ydata, xeval, sigma, null_thresh=0.6):
    """Apply gaussian sum filter to data.
    
    xdata, ydata : array
        Arrays of x- and y-coordinates of data. 
        Must be 1d and have the same length.
    
    xeval : array
        Array of x-coordinates at which to evaluate the smoothed result
    
    sigma : float
        Standard deviation of the Gaussian to apply to each data point
        Larger values yield a smoother curve.
    
    null_thresh : float
        For evaluation points far from data points, the estimate will be
        based on very little data. If the total weight is below this threshold,
        return np.nan at this location. Zero means always return an estimate.
        The default of 0.6 corresponds to approximately one sigma away 
        from the nearest datapoint.
    """
    # Distance between every combination of xdata and xeval
    # each row corresponds to a value in xeval
    # each col corresponds to a value in xdata
    delta_x = xeval[:, None] - xdata

    # Calculate weight of every value in delta_x using Gaussian
    # Maximum weight is 1.0 where delta_x is 0
    weights = np.exp(-0.5 * ((delta_x / sigma) ** 2))

    # Multiply each weight by every data point, and sum over data points
    smoothed = np.dot(weights, ydata)

    # Nullify the result when the total weight is below threshold
    # This happens at evaluation points far from any data
    # 1-sigma away from a data point has a weight of ~0.6
    nan_mask = weights.sum(1) < null_thresh
    smoothed[nan_mask] = np.nan
    
    # Normalize by dividing by the total weight at each evaluation point
    # Nullification above avoids divide by zero warning shere
    smoothed = smoothed / weights.sum(1)

    
    return smoothed

Answer 2

我先说明一下，这个问题更像是数字信号处理（DSP）的问题，而不是编程问题……

不过，针对你的问题，有一个简单的两步解决方案。

第一步：重新采样数据

为了说明这个过程，我们可以创建一个随机的数据集，数据的采样间隔不一样：

import numpy as np
x = np.cumsum(np.random.randint(0,100,100))
y = np.random.normal(0,1,size=100)

这样就会得到类似于下面的结果：

我们可以使用简单的线性插值来重新采样这些数据：

nx = np.arange(x.max()) # choose new x axis sampling
ny = np.interp(nx,x,y) # generate y values for each x

这样就把我们的数据转换成了：

第二步：应用滤波器

在这个阶段，你可以使用一些通过 scipy 提供的工具，给数据应用一个高斯滤波器，并设置一个特定的sigma值：

import scipy.ndimage.filters as filters
fx = filters.gaussian_filter1d(ny,sigma=100)

将这个结果与原始数据进行对比，我们得到：

这里的 sigma 值决定了滤波器的宽度。

Answer 3

对于非常大的数据集，这种方法可能会出现问题，但你想要的正确计算方法可以这样进行：

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(0) # for repeatability
x = np.random.rand(30)
x.sort()
y = np.random.rand(30)

x_eval = np.linspace(0, 1, 11)
sigma = 0.1

delta_x = x_eval[:, None] - x
weights = np.exp(-delta_x*delta_x / (2*sigma*sigma)) / (np.sqrt(2*np.pi) * sigma)
weights /= np.sum(weights, axis=1, keepdims=True)
y_eval = np.dot(weights, y)

plt.plot(x, y, 'bo-')
plt.plot(x_eval, y_eval, 'ro-')
plt.show()