3

我有这个关于 scipy 的 splrep 函数的问题，我认为这是一个错误，所以我会发布每段代码，以便您可以在您的计算机上重现它。假设我想找到一些数据的 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)

如果您将所有内容绘制到这里，一切似乎都很好： 适合 b 样条 但是，数据点数量和结数的某些组合存在问题。例如，如果您使用 and 尝试上面的代码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可能对计算要求很高。所以我的问题是：tryexcept

  1. 你能重现这个问题吗？如果你能...
  2. 你知道发生了什么事吗？
4

2 回答 2

2

我创造了一种解决问题的方法。它可能与适合两个结之间的数据点数量较少有关，所以我所做的是将创建结数的行替换为：

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

这样，我确保节点之间有足够的数据点，因为我使用 x 向量的索引。

于 2013-01-26T16:34:28.983 回答
1

t您可以使用参数，而不是直接指定结（通过参数） s。该s参数控制样条曲线的平滑度。它通过增加结的数量来做到这一点，直到条件

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

满足（g 是平滑样条表示，w 是权重）。

于 2013-09-17T09:09:35.347 回答