1-D 插值#

分段线性插值#

如果您只需要线性（又称折线）插值，可以使用 numpy.interp 例程。它接受两个用于插值的数据数组 x 和 y，以及第三个数组 xnew，包含要在其上评估插值的点。

>>> import numpy as np
>>> x = np.linspace(0, 10, num=11)
>>> y = np.cos(-x**2 / 9.0)

构造插值

>>> xnew = np.linspace(0, 10, num=1001)
>>> ynew = np.interp(xnew, x, y)

并绘制它

>>> import matplotlib.pyplot as plt
>>> plt.plot(xnew, ynew, '-', label='linear interp')
>>> plt.plot(x, y, 'o', label='data')
>>> plt.legend(loc='best')
>>> plt.show()
../../_images/1D-1.png

numpy.interp 的一个局限性是它不允许控制外插。有关提供此类功能的替代例程，请参阅 B 样条插值部分。

三次样条#

当然，分段线性插值在数据点处会产生拐角，即线性片段连接的地方。为了生成更平滑的曲线，可以使用三次样条，其中插值曲线由具有匹配的一阶和二阶导数的立方分段组成。在代码中，这些对象通过 CubicSpline 类实例表示。使用数据数组 x 和 y 来构造实例，然后可以使用目标 xnew 值对其进行评估。

>>> from scipy.interpolate import CubicSpline
>>> spl = CubicSpline([1, 2, 3, 4, 5, 6], [1, 4, 8, 16, 25, 36])
>>> spl(2.5)
5.57

CubicSpline 对象的 __call__ 方法接受标量值和数组。它还接受第二个参数 nu，用于评估 nu 阶导数。例如，我们绘制样条的导数。

>>> from scipy.interpolate import CubicSpline
>>> x = np.linspace(0, 10, num=11)
>>> y = np.cos(-x**2 / 9.)
>>> spl = CubicSpline(x, y)

>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(4, 1, figsize=(5, 7))
>>> xnew = np.linspace(0, 10, num=1001)
>>> ax[0].plot(xnew, spl(xnew))
>>> ax[0].plot(x, y, 'o', label='data')
>>> ax[1].plot(xnew, spl(xnew, nu=1), '--', label='1st derivative')
>>> ax[2].plot(xnew, spl(xnew, nu=2), '--', label='2nd derivative')
>>> ax[3].plot(xnew, spl(xnew, nu=3), '--', label='3rd derivative')
>>> for j in range(4):
...     ax[j].legend(loc='best')
>>> plt.tight_layout()
>>> plt.show()
../../_images/1D-2.png

请注意，一阶和二阶导数根据构造是连续的，而三阶导数在数据点处会跳跃。

单调插值器#

三次样条在构造上是两次连续可微的。这可能导致样条函数在数据点之间振荡和“超调”。在这些情况下，另一种方法是使用所谓的*单调*三次插值器：这些插值器构造为仅一次连续可微，并尝试保留数据所隐含的局部形状。scipy.interpolate 提供了两种此类对象：PchipInterpolator 和 Akima1DInterpolator。为了说明，我们考虑一个包含异常值的数据。

>>> from scipy.interpolate import CubicSpline, PchipInterpolator, Akima1DInterpolator
>>> x = np.array([1., 2., 3., 4., 4.5, 5., 6., 7., 8])
>>> y = x**2
>>> y[4] += 101

>>> import matplotlib.pyplot as plt
>>> xx = np.linspace(1, 8, 51)
>>> plt.plot(xx, CubicSpline(x, y)(xx), '--', label='spline')
>>> plt.plot(xx, Akima1DInterpolator(x, y)(xx), '-', label='Akima1D')
>>> plt.plot(xx, PchipInterpolator(x, y)(xx), '-', label='pchip')
>>> plt.plot(x, y, 'o')
>>> plt.legend()
>>> plt.show()
../../_images/1D-3.png

B 样条插值#

B 样条是分段多项式的另一种（如果形式上等效）表示。该基通常比幂基在计算上更稳定，并且对于包括插值、回归和曲线表示在内的各种应用很有用。详细信息请参阅 分段多项式部分，在此我们通过构造正弦函数的插值来说明其用法。

>>> x = np.linspace(0, 3/2, 7)
>>> y = np.sin(np.pi*x)

要根据数据数组 x 和 y 构造插值对象，我们使用 make_interp_spline 函数。

>>> from scipy.interpolate import make_interp_spline
>>> bspl = make_interp_spline(x, y, k=3)

此函数返回一个对象的接口与 CubicSpline 对象类似。特别是，它可以被评估到数据点并进行微分。

>>> der = bspl.derivative()      # a BSpline representing the derivative
>>> import matplotlib.pyplot as plt
>>> xx = np.linspace(0, 3/2, 51)
>>> plt.plot(xx, bspl(xx), '--', label=r'$\sin(\pi x)$ approx')
>>> plt.plot(x, y, 'o', label='data')
>>> plt.plot(xx, der(xx)/np.pi, '--', label=r'$d \sin(\pi x)/dx / \pi$ approx')
>>> plt.legend()
>>> plt.show()
../../_images/1D-4.png

请注意，通过在 make_interp_spline 调用中指定 k=3，我们要求一个三次样条（这是默认值，因此 k=3 可以省略）；三次导数的导数是二次的。

>>> bspl.k, der.k
(3, 2)

默认情况下，make_interp_spline(x, y) 的结果等同于 CubicSpline(x, y)。区别在于前者允许几种可选功能：它可以通过可选参数 k 来构造各种次数的样条，并通过可选参数 t 来构造预定义的节点。

样条插值的边界条件可以通过 make_interp_spline 函数和 CubicSpline 构造函数的 bc_type 参数进行控制。默认情况下，两者都使用“非节点”边界条件。

非三次样条#

make_interp_spline 的一个用途是构造具有线性外插的线性插值器，因为 make_interp_spline 默认会进行外插。考虑。

>>> from scipy.interpolate import make_interp_spline
>>> x = np.linspace(0, 5, 11)
>>> y = 2*x
>>> spl = make_interp_spline(x, y, k=1)  # k=1: linear
>>> spl([-1, 6])
[-2., 12.]
>>> np.interp([-1, 6], x, y)
[0., 10.]

有关更多详细信息和讨论，请参阅 外插部分。

y 的批次#

单变量插值器不仅接受一维 y 数组，还接受 y.ndim > 1。解释是 y 是 1D 数据数组的*批次*：默认情况下，y 的零维是插值轴，尾随维度是批次维度。考虑一组（一个*批次*）函数 \(f_j\) 在点 \(x_i\) 采样。通过提供一个二维数组 y，其中 y[i, j] 记录 \(f_j(x_i)\)，我们可以为所有这些函数实例化一个插值器。

为了说明

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.interpolate import make_interp_spline
>>> n = 11
>>> x = 2 * np.pi * np.arange(n) / n
>>> x.shape
(11,)
>>> y = np.stack((np.sin(x)**2, np.cos(x)), axis=1)
>>> y.shape
(11, 2)
>>> spl = make_interp_spline(x, y)
>>> xv = np.linspace(0, 2*np.pi, 51)
>>> plt.plot(x, y, 'o')
>>> plt.plot(xv, spl(xv), '-')
>>> plt.show()

有几点需要注意。最重要的一点是，这里的行为与 NumPy 的广播看起来相似，但在两方面有所不同：

  1. 即使 y 数组不是一维的，x 数组也应该是一维的：x.ndim == 1，而 y.ndim >= 1。x 与 y 之间没有广播。

  2. 默认情况下，*尾随*维度用作批次维度，这与 NumPy 使用*前导*维度作为批次维度的约定不同。

其次，插值轴可以通过可选的 axis 参数进行控制。上面的示例使用默认值 axis=0。对于非默认值，以下情况成立：

  • y.shape[axis] == x.size（否则会引发错误）

  • spl(xv) 的形状是 y.shape[axis:] + xv.shape + y.shape[:axis]。

虽然我们使用 make_interp_spline 演示了批处理行为，但实际上大多数单变量插值器都支持此功能：PchipInterpolator 和 Akima1DInterpolator、CubicSpline；低级多项式表示类，PPoly、BPoly 和 BSpline；以及 最小二乘拟合和样条平滑函数，make_lsq_spline 和 make_smoothing_spline。

参数化样条曲线#

到目前为止，我们考虑了样条*函数*，其中数据 y 期望显式依赖于自变量 x — 因此插值函数满足 \(f(x_j) = y_j\)。样条*曲线*将 x 和 y 数组视为平面上的点 \(\mathbf{p}_j\) 的坐标，并由附加参数（通常称为 u）参数化一条通过这些点的插值曲线。请注意，此构造可以轻松地推广到高维空间，其中 \(\mathbf{p}_j\) 是 N 维空间中的点。

如在 上一节 中所述，通过利用插值函数处理多维数据数组的事实，可以轻松地构造样条曲线。数据点对应的参数 u 的值需要用户单独提供。

参数化的选择取决于具体问题，不同的参数化可能会产生截然不同的曲线。例如，我们考虑（一个相当困难的）数据集的三个参数化，该数据集取自 Ref [1] 的第 6 章，列在 BSpline 文档字符串中。

>>> x = [0, 1, 2, 3, 4, 5, 6]
>>> y = [0, 0, 0, 9, 0, 0, 0]
>>> p = np.stack((x, y))
>>> p
array([[0, 1, 2, 3, 4, 5, 6],
       [0, 0, 0, 9, 0, 0, 0]])

我们将 p 数组的元素作为平面上七个点的坐标，其中 p[:, j] 给出了点 \(\mathbf{p}_j\) 的坐标。

首先，考虑*均匀*参数化，\(u_j = j\)。

>>> u_unif = x

其次，我们考虑所谓的*弦长*参数化，它只是连接数据点的直线段的累积长度。

\[u_j = u_{j-1} + |\mathbf{p}_j - \mathbf{p}_{j-1}|\]

对于 \(j=1, 2, \dots\) 和 \(u_0 = 0\)。此处 \(| \cdots |\) 是平面上连续点 \(p_j\) 之间的长度。

>>> dp = p[:, 1:] - p[:, :-1]      # 2-vector distances between points
>>> l = (dp**2).sum(axis=0)        # squares of lengths of 2-vectors between points
>>> u_cord = np.sqrt(l).cumsum()   # cumulative sums of 2-norms
>>> u_cord = np.r_[0, u_cord]      # the first point is parameterized at zero

最后，我们考虑有时称为*向心*参数化：\(u_j = u_{j-1} + |\mathbf{p}_j - \mathbf{p}_{j-1}|^{1/2}\)。由于额外的平方根，连续值 \(u_j - u_{j-1}\) 之间的差异将小于弦长参数化。

>>> u_c = np.r_[0, np.cumsum((dp**2).sum(axis=0)**0.25)]

现在绘制生成的曲线。

>>> from scipy.interpolate import make_interp_spline
>>> import matplotlib.pyplot as plt
>>> fig, ax = plt.subplots(1, 3, figsize=(8, 3))
>>> parametrizations = ['uniform', 'cord length', 'centripetal']
>>>
>>> for j, u in enumerate([u_unif, u_cord, u_c]):
...    spl = make_interp_spline(u, p, axis=1)    # note p is a 2D array
...
...    uu = np.linspace(u[0], u[-1], 51)
...    xx, yy = spl(uu)
...
...    ax[j].plot(xx, yy, '--')
...    ax[j].plot(p[0, :], p[1, :], 'o')
...    ax[j].set_title(parametrizations[j])
>>> plt.show()
../../_images/1D-5.png

缺失数据#

我们注意到 scipy.interpolate*不支持*带有缺失数据的插值。表示缺失数据的两种流行方法是使用 numpy.ma 库的掩码数组，以及将缺失值编码为非数字 NaN。

这两种方法均未在 scipy.interpolate 中直接支持。单个例程可能提供部分支持和/或解决方法，但总的来说，该库坚决遵守 IEEE 754 语义，其中 NaN 表示*非数字*，即非法数学运算（例如除以零）的结果，而不是*缺失*。

1-D 插值的旧版接口（interp1d）#

注意

interp1d 被视为旧版 API，不建议在新代码中使用。请考虑使用 更具体的插值器。

scipy.interpolate 中的 interp1d 类是一种便捷的方法，可根据固定数据点创建函数，该函数可以使用线性插值在由给定数据定义的任何域内进行评估。通过传递构成数据的 1-D 向量来创建此类的实例。此类的一个实例定义了一个 __call__ 方法，因此可以将其视为一个函数，该函数在已知数据值之间进行插值以获得未知值。可以在实例化时指定边界行为。以下示例演示了其用于线性三次样条插值。

>>> from scipy.interpolate import interp1d

>>> x = np.linspace(0, 10, num=11, endpoint=True)
>>> y = np.cos(-x**2/9.0)
>>> f = interp1d(x, y)
>>> f2 = interp1d(x, y, kind='cubic')

>>> xnew = np.linspace(0, 10, num=41, endpoint=True)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'o', xnew, f(xnew), '-', xnew, f2(xnew), '--')
>>> plt.legend(['data', 'linear', 'cubic'], loc='best')
>>> plt.show()
"This code generates an X-Y plot of a time-series with amplitude on the Y axis and time on the X axis. The original time-series is shown as a series of blue markers roughly defining some kind of oscillation. An orange trace showing the linear interpolation is drawn atop the data forming a jagged representation of the original signal. A dotted green cubic interpolation is also drawn that appears to smoothly represent the source data."

interp1d 的“cubic”类型等同于 make_interp_spline，而“linear”类型等同于 numpy.interp，同时还允许 N 维 y 数组。

interp1d 中的另一组插值是*nearest*、*previous*和*next*，它们返回 x 轴上最近、前一个或下一个点。Nearest 和 next 可以看作是因果插值滤波器的特例。以下示例使用与上一个示例相同的数据来说明它们的用法。

>>> from scipy.interpolate import interp1d

>>> x = np.linspace(0, 10, num=11, endpoint=True)
>>> y = np.cos(-x**2/9.0)
>>> f1 = interp1d(x, y, kind='nearest')
>>> f2 = interp1d(x, y, kind='previous')
>>> f3 = interp1d(x, y, kind='next')

>>> xnew = np.linspace(0, 10, num=1001, endpoint=True)
>>> import matplotlib.pyplot as plt
>>> plt.plot(x, y, 'o')
>>> plt.plot(xnew, f1(xnew), '-', xnew, f2(xnew), '--', xnew, f3(xnew), ':')
>>> plt.legend(['data', 'nearest', 'previous', 'next'], loc='best')
>>> plt.show()
"This code generates an X-Y plot of a time-series with amplitude on the Y axis and time on the X axis. The original time-series is shown as a series of blue markers roughly defining some kind of oscillation. An orange trace showing the nearest neighbor interpolation is drawn atop the original with a stair-like appearance where the original data is right in the middle of each stair step. A green trace showing the previous neighbor interpolation looks similar to the orange trace but the original data is at the back of each stair step. Similarly a dotted red trace showing the next neighbor interpolation goes through each of the previous points, but it is centered at the front edge of each stair."