Python的interp1d与UnivariateSpline比较

Question

我正在尝试把一些MatLab的代码移植到Scipy上，试了两个来自scipy.interpolate的函数，一个是interp1d，另一个是UnivariateSpline。使用interp1d的结果和MatLab中的interp1d函数一致，但使用UnivariateSpline得到的结果却不一样，有时候差别还挺大的。

f = interp1d(row1,row2,kind='cubic',bounds_error=False,fill_value=numpy.max(row2))
return  f(interp)

f = UnivariateSpline(row1,row2,k=3,s=0)
return f(interp)

有没有人能给点建议？我的x值不是均匀分布的，不过我不太明白这会有什么影响。

Answer 1

对我来说是有效的，

from scipy import allclose, linspace
from scipy.interpolate import interp1d, UnivariateSpline

from numpy.random import normal

from pylab import plot, show

n = 2**5

x = linspace(0,3,n)
y = (2*x**2 + 3*x + 1) + normal(0.0,2.0,n)

i = interp1d(x,y,kind=3)
u = UnivariateSpline(x,y,k=3,s=0)

m = 2**4

t = linspace(1,2,m)

plot(x,y,'r,')
plot(t,i(t),'b')
plot(t,u(t),'g')

print allclose(i(t),u(t)) # evaluates to True

show()

这让我得到了，

Answer 2

结果之所以不同（但都可能是正确的），是因为UnivariateSpline和interp1d使用的插值方法不一样。

• interp1d会用你提供的x点作为节点，构建一个平滑的B样条曲线。

• UnivariateSpline是基于FITPACK的，它同样构建平滑的B样条曲线。不过，FITPACK会尝试选择一些新的节点，以便更好地拟合数据（可能是为了最小化某种误差，或者考虑曲率的惩罚，类似的意思）。你可以通过g.get_knots()查看它使用了哪些节点。

所以，你得到不同结果的原因就是插值算法不同。如果你想要在数据点上有节点的B样条，可以使用interp1d或splmake。如果你想要FITPACK的效果，就用UnivariateSpline。在数据密集的情况下，两种方法的结果是一样的，但当数据稀疏时，结果可能会不同。

（我怎么知道这些的呢：我看了代码 :-)

Answer 3

我刚遇到同样的问题。

简短回答

可以使用 InterpolatedUnivariateSpline 来解决：

f = InterpolatedUnivariateSpline(row1, row2)
return f(interp)

详细回答

UnivariateSpline 是一种“针对给定数据点的一维平滑样条曲线”，而 InterpolatedUnivariateSpline 则是“针对给定数据点的一维插值样条曲线”。前者是用来平滑数据的，而后者是一种更常见的插值方法，能够产生与 interp1d 相同的结果。下面的图展示了这两者的区别。

下面的代码可以用来重现这张图。

import scipy.interpolate as ip

#Define independent variable
sparse = linspace(0, 2 * pi, num = 20)
dense = linspace(0, 2 * pi, num = 200)

#Define function and calculate dependent variable
f = lambda x: sin(x) + 2
fsparse = f(sparse)
fdense = f(dense)

ax = subplot(2, 1, 1)

#Plot the sparse samples and the true function
plot(sparse, fsparse, label = 'Sparse samples', linestyle = 'None', marker = 'o')
plot(dense, fdense, label = 'True function')

#Plot the different interpolation results
interpolate = ip.InterpolatedUnivariateSpline(sparse, fsparse)
plot(dense, interpolate(dense), label = 'InterpolatedUnivariateSpline', linewidth = 2)

smoothing = ip.UnivariateSpline(sparse, fsparse)
plot(dense, smoothing(dense), label = 'UnivariateSpline', color = 'k', linewidth = 2)

ip1d = ip.interp1d(sparse, fsparse, kind = 'cubic')
plot(dense, ip1d(dense), label = 'interp1d')

ylim(.9, 3.3)

legend(loc = 'upper right', frameon = False)

ylabel('f(x)')

#Plot the fractional error
subplot(2, 1, 2, sharex = ax)

plot(dense, smoothing(dense) / fdense - 1, label = 'UnivariateSpline')
plot(dense, interpolate(dense) / fdense - 1, label = 'InterpolatedUnivariateSpline')
plot(dense, ip1d(dense) / fdense - 1, label = 'interp1d')

ylabel('Fractional error')
xlabel('x')
ylim(-.1,.15)

legend(loc = 'upper left', frameon = False)

tight_layout()