四联光电智能照明论坛

标题: Python插值函数学习笔记 [打印本页]

作者: admin 时间: 2016-7-21 11:49
标题: Python插值函数学习笔记
插值-interpolate：http://hyry.dip.jp/tech/book/pag ... 10-interpolate.html

插值是通过已知的离散数据求未知数据的方法。与拟合不同的是，它要求曲线通过所有的已知数据。SciPy的interpolate模块提供了许多对数据进行插值运算的函数。

一维插值
一维数据的插值运算可以通过interp1d()完成。其调用形式如下，它实际上不是函数而是一个类：

interp1d(x, y, kind='linear', ...)
其中x和y参数是一系列已知的数据点，kind参数是插值类型，它可以是字符串或整数，它给出插值的B样条曲线的阶数，可以有如下侯选值：

'zero'、'nearest'：阶梯插值，相当于0阶B样条曲线。
'slinear'、'linear'：线性插值，用一条直线连接所有的取样点，相当于1阶B样条曲线，'slinear'使用扩展库中的相关函数计算，而'linear'则直接使用Python编写的函数进行运算，其结果一样。
'quadratic'、'cubic'：2阶和3阶B样条曲线，更高阶的曲线可以直接使用整数值指定。
interp1d对象可以计算x的取值范围之内任意点的函数值。它可以像函数一样直接调用，它和NumPy的ufunc函数一样能对数组中的每个元素进行计算，并返回一个新的数组。

下面的程序演示了kind参数与其对应的插值曲线。程序中我们使用循环对相同的数据进行4种不同阶数的插值运算。❶首先使用数据点创建一个interp1d对象f，通过kind参数指定其阶数。❷调用f()计算出一系列的插值结果。本例中，决定插值曲线的数据点一共有11个，插值之后的曲线数据点有101个。

import numpy as np
from scipy import interpolate
import pylab as pl

x = np.linspace(0, 10, 11)
y = np.sin(x)

xnew = np.linspace(0, 10, 101)
plt.plot(x,y,'ro')
for kind in ['nearest', 'zero', 'slinear', 'quadratic']:
f = interpolate.interp1d(x,y,kind=kind) ❶
ynew = f(xnew) ❷
plt.plot(xnew, ynew, label=str(kind))

plt.legend(loc='lower right');

外推和Spline拟合
上节所介绍的interp1d类要求其参数x是一个递增的序列，并且只能在x的取值范围之内进行内插计算，不能用它进行外推运算，即计算x的取值范围之外的数据点。UnivariateSpline类的插值运算比interp1d更高级，它支持外推和拟合运算，其调用形式如下：

UnivariateSpline(x, y, w=None, bbox=[None, None], k=3, s=None)
x、y是保存数据点的X-Y坐标的数组，其中x必须是递增序列。
w是为每个数据点指定的权重值。
k为样条曲线的阶数。
s是平滑系数，它使得最终生成的样条曲线满足条件：\sum_{}^{}{(w \cdot (y-spline(x)))^2} \leq s ，即当s>0时，样条曲线并不一定通过各个数据点。为了让曲线通过所有数据点，必须将s参数设置为0。
下面的程序演示使用UnivariateSpline对数据进行插值、外推以及样条曲线拟合：

❶如%Fig(上)所示，UnivariateSpline能够进行外推运算，虽然输入数据中没有X轴大于10的点，但是它能计算出X轴在0到12的插值结果。在X轴大于10的部分，样条曲线仍然呈现出和正弦波类似的形状，越远离输入数据范围，误差会越大，因此外推的范围是有限的。由于s参数为0，因此插值曲线经过所有的数据点。

❷%Fig(下)则显示了s参数不为零时的结果，对于带噪声的输入数据，选择合适的s参数能够使得样条曲线接近无噪声时的波形，可以把它看作是使用样条曲线对数据进行拟合运算。

%Fig=使用UnivariateSpline进行插值：外推(上)，数据拟合(下)

x1 = np.linspace(0, 10, 20)
y1 = np.sin(x1)
sx1 = np.linspace(0, 12, 100)
sy1 = interpolate.UnivariateSpline(x1, y1, s=0)(sx1) ❶

x2 = np.linspace(0, 20, 200)
y2 = np.sin(x2) + np.random.standard_normal(len(x2))*0.2
sx2 = np.linspace(0, 20, 2000)
sy2 = interpolate.UnivariateSpline(x2, y2, s=8)(sx2) ❷

plt.figure(figsize=(8, 5))
plt.subplot(211)
plt.plot(x1, y1, ".", label=u"数据点")
plt.plot(sx1, sy1, label=u"spline曲线")
plt.legend()

plt.subplot(212)
plt.plot(x2, y2, ".", label=u"数据点")
plt.plot(sx2, sy2, linewidth=2, label=u"spline曲线")
plt.plot(x2, np.sin(x2), label=u"无噪声曲线")
plt.legend()

参数插值
前面所介绍的插值函数都需要X轴的数据是按照递增顺序排列的，就像一般的y=f(x)函数曲线一样。数学上还有一种参数曲线，它使用一个参数t和两个函数x=f(t), y=g(t)，定义二维平面上的一条曲线，例如圆形、心形等曲线都是参数曲线。参数曲线的插值可以通过splprep()和splev()实现，这组函数支持高维空间的曲线的插值，这里以二维曲线为例，介绍其用法。

❶首先调用splprep()，其第一个参数为一组一维数组，每个数组是各点在对应轴上的坐标。s参数为平滑系数，与UnivariateSpline的含义相同。splprep()返回两个对象，其中tck是一个元组，其中包含了插值曲线的所有信息。t是自动计算出的参数曲线的参数数组。

❷调用splev()进行插值运算，其第一个参数为一个新的参数数组，这里将t的取值范围等分200份，第二个参数为splprep()返回的第一个对象。实际上，参数数组t是正规化之后的各个线段长度的累计，因此t的范围位0到1。

其结果如图所示，图中比较了平滑系数为0和1e-4时的插值曲线。当平滑系数为0时，插值曲线通过所有的数据点。

x = [ 4.913,  4.913,  4.918,  4.938,  4.955,  4.949,  4.911,
   4.848,  4.864,  4.893,  4.935,  4.981,  5.01 ,  5.021]

y = [ 5.2785,  5.2875,  5.291 ,  5.289 ,  5.28  ,  5.26  ,  5.245 ,
   5.245 ,  5.2615,  5.278 ,  5.2775,  5.261 ,  5.245 ,  5.241]

plt.plot(x, y, "o")

for s in (0, 1e-4):
tck, t = interpolate.splprep([x, y], s=s) ❶
xi, yi = interpolate.splev(np.linspace(t[0], t[-1], 200), tck) ❷
plt.plot(xi, yi, lw=2, label=u"s=%g" % s)

plt.legend();

单调插值
前面介绍的几种插值方法，不能保证数据点的单调性，即曲线的最值可能出现在数据点之外的地方。PchipInterpolator类(别名pchip)使用单调三次插值，能够保证曲线的所有最值都出现在数据点之上。下面的程序用pchip()对数据点进行插值，并绘制其一阶导数曲线，由下图的导数曲线可知，所有最值数据点出的导数都为0。

x = [0, 1, 2, 3, 4, 5]
y = [1, 2, 1.5, 2.5, 3, 2.5]
xs = np.linspace(x[0], x[-1], 100)
curve = interpolate.pchip(x, y)
ys = curve(xs)
dys = curve.derivative(xs)
plt.plot(xs, ys, label=u"pchip")
plt.plot(xs, dys, label=u"一阶导数")
plt.plot(x, y, "o")
plt.legend(loc="best")
plt.grid()
plt.margins(0.1, 0.1)

作者: admin 时间: 2016-7-21 11:53
sci.interpolate的所有函数：
'Akima1DInterpolator',
'BPoly',
'BarycentricInterpolator',
'BivariateSpline',
'CloughTocher2DInterpolator',
'InterpolatedUnivariateSpline',
'KroghInterpolator',
'LSQBivariateSpline',
'LSQSphereBivariateSpline',
'LSQUnivariateSpline',
'LinearNDInterpolator',
'NearestNDInterpolator',
'PPoly', 'PchipInterpolator',
'PiecewisePolynomial',
'Rbf',
'RectBivariateSpline',
'RectSphereBivariateSpline',
'RegularGridInterpolator',
'SmoothBivariateSpline',
'SmoothSphereBivariateSpline',
'Tester',
'UnivariateSpline',
'__all__',
'__builtins__',
'__cached__',
'__doc__',
'__file__',
'__loader__',
'__name__',
'__package__',
'__path__',
'__spec__',
'_fitpack',
'_monotone',
'_ppoly',
'absolute_import',
'approximate_taylor_polynomial',
'barycentric_interpolate',
'bench',
'bisplev',
'bisplrep',
'dfitpack',
'division',
'fitpack',
'fitpack2',
'griddata',
'insert',
'interp1d',
'interp2d',
'interpn',
'interpnd',
'interpolate',
'krogh_interpolate',
'lagrange',
'ndgriddata',
'pchip',
'pchip_interpolate',
'piecewise_polynomial_interpolate',
'polyint',
'ppform',
'print_function',
'rbf',
'spalde',
'splantider',
'splder',
'splev',
'spleval',
'spline',
'splint',
'splmake',
'splprep',
'splrep',
'spltopp',
'sproot',
'test'

欢迎光临四联光电智能照明论坛 (http://5xhome.com/)