如何获取加权高斯滤波器

Question

我有一组带权重的 x,y 点，如下所示（完整数据集在这里）：

#  x       y     w
-0.038  2.0127  0.71
0.058   1.9557  1
0.067   2.0016  0.9
0.072   2.0316  0.83
...

我需要找到一条平滑的线，根据分配给每个点的重要性进行调整，即更高的权重意味着数据点应该具有更高的相关性。

这是我目前的代码，基本上将一个高斯滤波器应用于数据（我从这个问题中得到了灵感：Python中的线平滑算法？）：

import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import gaussian_filter1d

# Read data from file.
data = np.loadtxt('data_file', unpack=True)
x, y, w = data[0], data[1], data[2]

# Return evenly spaced numbers over a specified interval.
t = np.linspace(0, 1, len(x))
t2 = np.linspace(0, 1, 100)    
# One-dimensional linear interpolation.
x2 = np.interp(t2, t, x)
y2 = np.interp(t2, t, y)

# Obtain Gaussian filter with fixed sigma value.
sigma = 7
x3 = gaussian_filter1d(x2, sigma)
y3 = gaussian_filter1d(y2, sigma)

# Make plot.
cm = plt.cm.get_cmap('RdYlBu')
plt.scatter(x, y, marker="o", c=w, s=40, cmap=cm, lw=0.5, vmin=0, vmax=1)
plt.plot(x3, y3, "r", lw=2)
plt.show()

这段代码生成以下绘图（蓝色点具有更高的权重值）：

问题在于此拟合并未考虑到每个点分配的权重。如何将该信息引入高斯滤波器中？

- Gabriel

---

许多scipy插值允许指定权重。http://docs.scipy.org/doc/scipy/reference/interpolate.html - tom10

@tom10，您方便指向一个具体的方案吗？据我所知，UnivariateSpline函数有权重，但需要x数组递增，而我的数据不符合该条件。 - Gabriel

你是否知道numpy.interp文档中所说的：“xp:数据点的x坐标，必须单调递增。”和“不会检查x坐标序列xp是否单调递增。如果xp不是单调递增的，则结果是无意义的。”所以，对于任何一个，最好从排序开始（使用numpy.argsort以便您可以在所有轴上都与x一致）。 - tom10

@tom10，我不确定我理解你的意思。如果我对数据进行排序，那么这将影响我试图拟合/插值的点云的形状。您能否详细说明一下您的答案是什么意思？谢谢。 - Gabriel

很不幸，我没有在任何我能够访问的计算机上安装scipy，所以我无法提供完整的示例。但是，简单来说：排序不应该改变云的形状，它只是意味着x值是有序的。也就是说，如果你有一堆(x,y)对，每个对中的值创建了云的形状，你只需要按照它们的x值对这些值进行排序。因此，对于分离的x、y、w数组，你需要按照x对所有的x、y和w进行排序：例如，ysorted = y[np.argsor(x)]，即，这些对是相同的，但是重新排序了。 - tom10

@tom10 哦，我现在明白你的意思了！一旦我弄清楚如何按照你所提到的方式对数组进行排序，我就会尽快修复问题。非常感谢你指出这一点，我完全忘记了它！ - Gabriel

Answer 1

请注意，以下思路是一种解决方法而不是精确的解决方案，但值得一试。

思路是使用“w”权重参数来重复相应的“x”和“y”值。因此，如果您将“w”缩放到例如范围“[1,10]”，则“x”中的所有相应值以及“y”中的所有相应值都将被重复10次，当“w”等于10时。也就是说，新的“x”和“y”将被创建。通过这种方式，我们将重量作为值在“x”和“y”中的频率纳入其中。完成这些操作后，将新值提供给算法，希望能够得到所需的结果，如下面的示例所示。

对于第一个图，蓝色到红色的光谱对应于较低到较高的权重。标题中的数字是上述“复制因子”。

对于第二个图，我们没有更改您的颜色格式。

- Developer

---

虽然这不完全是我所要求的，但似乎符合我的需求，因此我将其标记为已接受。非常感谢@开发者！ - Gabriel