scipy.interpolate的splrep函数选择节点时的bug(?)

Question

我有一个关于scipy库中splrep函数的问题，我觉得这可能是个bug，所以我会把所有的代码贴出来，方便你们在自己的电脑上复现。假设我想找到一些数据的b样条表示，比如说，下面这段代码生成的数据集，它是十个高斯分布的混合，并且加了一些噪声：

import numpy as np
# First we define the number of datapoints:
ndata = 100
x = np.arange(0,1,1./np.double(ndata))
means = np.random.uniform(0,1,10)
y = 0.0
for i in range(len(means)):
    y = y+np.exp(-(x-means[i])**2./0.01)
# We add some noise to obtain the data:
data = y + np.random.normal(0,0.05,len(y))

生成的数据应该是这样的：

现在，我们用splrep和splev函数来获取这条曲线的b样条表示：

from scipy.interpolate import splrep,splev
# First define the number of knots. Let's put, say, 10 knots:
nknots = 10
# Now we crate the array of knots:
knots = np.arange(x[1],x[len(x)-1],(x[len(x)-1]-x[1])/np.double(nknots))
tck = splrep(x,data,t=knots)
fit = splev(x,tck)

如果你把到目前为止的所有内容都画出来，结果看起来是正常的：

但是，当数据点的数量和节点的数量某些组合时，会出现问题。例如，如果你用ndata = 1931和nknots = 796来运行上面的代码，我会遇到以下错误：

File "/usr/lib/python2.7/dist-packages/scipy/interpolate/fitpack.py", line 465, in
splrep raise _iermess[ier][1](_iermess[ier][0])
ValueError:     Error on input data

这给我带来了麻烦，因为上面的代码不能自动化处理。我正在处理大约19000个数据点的数据集，在这种情况下，使用while循环结合try和except会非常消耗计算资源。所以我有几个问题：

1. 你能复现这个问题吗？如果可以的话……
2. 你知道发生了什么吗？

Answer 1

你可以不直接指定节点（通过 t 参数），而是使用 s 参数。
s 参数用来控制样条曲线的光滑程度。它通过增加节点的数量来实现这一点，直到满足以下条件：

sum((w * (y - g))**2,axis=0) <= s

（g 是光滑后的样条曲线表示，w 是权重）。

Answer 2

我想出了一个解决这个问题的方法。这个问题可能是因为在两个节点之间的数据点太少了。所以我做的就是把创建节点数量的那一行代码替换成：

idx_knots = (np.arange(1,len(x)-1,(len(x)-2)/np.double(nknots))).astype('int')
knots = x[idx_knots]

这样一来，我就能确保在节点之间有足够的数据点，因为我调整了x向量的索引。