scipy.interpolate.UnivariateSpline 无论参数如何都不平滑

Question

我在使用scipy.interpolate.UnivariateSpline进行插值时，遇到了无法实现平滑的问题。根据这个函数的页面和一些之前的帖子，我认为它应该可以通过参数来提供平滑效果。

这是我的代码：

# Imports
import scipy
import pylab

# Set up and plot actual data
x = [0, 5024.2059124920379, 7933.1645067836089, 7990.4664106277542, 9879.9717114947653, 13738.60563208926, 15113.277958924193]
y = [0.0, 3072.5653360000988, 5477.2689107965398, 5851.6866463790966, 6056.3852496014106, 7895.2332350173638, 9154.2956175610598]
pylab.plot(x, y, "o", label="Actual")

# Plot estimates using splines with a range of degrees
for k in range(1, 4):
    mySpline = scipy.interpolate.UnivariateSpline(x=x, y=y, k=k, s=2)
    xi = range(0, 15100, 20)
    yi = mySpline(xi)
    pylab.plot(xi, yi, label="Predicted k=%d" % k)

# Show the plot
pylab.grid(True)
pylab.xticks(rotation=45)
pylab.legend( loc="lower right" )
pylab.show()

这是结果：

我尝试了不同的值（0.01, 0.1, 1, 2, 5, 50），还尝试了显式权重，设置为相同的值（1.0）或随机值。结果仍然没有平滑效果，而且节点的数量总是和数据点的数量一样。特别是，我希望像第4个点（7990.4664106277542, 5851.6866463790966）这样的异常值能够被平滑处理。

这是不是因为我的数据不够多？如果是这样，有没有类似的样条函数或聚类技术可以用来在这么少的数据点上实现平滑？

Answer 1

虽然我不知道有没有现成的库可以帮你做到这一点，但我建议你可以自己动手试试：首先可以在原始数据点之间创建一个样条曲线，也就是在 x 和 y 轴上都加一些节点。在你的例子中，在第4个和第5个点之间加一个节点就可以了，这样可以消除在 x=8000 附近的那个大变化。

Answer 2

简单来说，你需要更仔细地选择 s 的值。

关于 UnivariateSpline 的说明中提到：

Positive smoothing factor used to choose the number of knots. Number of 
knots will be increased until the     smoothing condition is satisfied:
sum((w[i]*(y[i]-s(x[i])))**2,axis=0) <= s

从中我们可以推测，如果你没有提供明确的权重，那么“合理”的平滑值大约是 s = m * v，其中 m 是数据点的数量，v 是数据的方差。在这种情况下，s_good ~ 5e7。

编辑：s 的合理值当然还取决于数据中的噪声水平。文档似乎建议选择 s 的范围是 (m - sqrt(2*m)) * std**2 <= s <= (m + sqrt(2*m)) * std**2，其中 std 是与想要平滑的“噪声”相关的标准差。

Answer 3

@Zhenya提到的手动在数据点之间设置节点的方法，对于噪声数据来说效果不太好，因为这种方法需要很小心地使用。不过，受到他/她建议的启发，我尝试了scikit-learn这个包里的均值漂移聚类，结果很不错。这个方法可以自动确定聚类的数量，而且在平滑数据方面表现得相当好（实际上非常平滑）。

# Imports
import numpy
import pylab
import scipy
import sklearn.cluster

# Set up original data - note that it's monotonically increasing by X value!
data = {}
data['original'] = {}
data['original']['x'] = [0, 5024.2059124920379, 7933.1645067836089, 7990.4664106277542, 9879.9717114947653, 13738.60563208926, 15113.277958924193]
data['original']['y'] = [0.0, 3072.5653360000988, 5477.2689107965398, 5851.6866463790966, 6056.3852496014106, 7895.2332350173638, 9154.2956175610598]

# Cluster data, sort it and and save
inputNumpy = numpy.array([[data['original']['x'][i], data['original']['y'][i]] for i in range(0, len(data['original']['x']))])
meanShift = sklearn.cluster.MeanShift()
meanShift.fit(inputNumpy)
clusteredData = [[pair[0], pair[1]] for pair in meanShift.cluster_centers_]
clusteredData.sort(lambda pair1, pair2: cmp(pair1[0],pair2[0]))
data['clustered'] = {}
data['clustered']['x'] = [pair[0] for pair in clusteredData]
data['clustered']['y'] = [pair[1] for pair in clusteredData]

# Build a spline using the clustered data and predict
mySpline = scipy.interpolate.UnivariateSpline(x=data['clustered']['x'], y=data['clustered']['y'], k=1)
xi = range(0, round(max(data['original']['x']), -3) + 3000, 20)
yi = mySpline(xi)

# Plot the datapoints
pylab.plot(data['clustered']['x'], data['clustered']['y'], "D", label="Datapoints (%s)" % 'clustered')
pylab.plot(xi, yi, label="Predicted (%s)" %  'clustered')
pylab.plot(data['original']['x'], data['original']['y'], "o", label="Datapoints (%s)" % 'original')

# Show the plot
pylab.grid(True)
pylab.xticks(rotation=45)
pylab.legend( loc="lower right" )
pylab.show()