四联光电智能照明论坛

 找回密码
 立即注册
搜索
热搜: 活动 交友 discuz
查看: 4630|回复: 1
打印 上一主题 下一主题

Python插值函数学习笔记

[复制链接]
  • TA的每日心情
    开心
    2018-12-28 16:25
  • 817

    主题

    1556

    帖子

    1万

    积分

    管理员

    Rank: 9Rank: 9Rank: 9

    积分
    14941
    跳转到指定楼层
    楼主
    发表于 2016-7-21 11:49:55 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
    插值-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)
  • TA的每日心情
    开心
    2018-12-28 16:25
  • 817

    主题

    1556

    帖子

    1万

    积分

    管理员

    Rank: 9Rank: 9Rank: 9

    积分
    14941
    沙发
     楼主| 发表于 2016-7-21 11:53:33 | 只看该作者
    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'
    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    QQ|Archiver|手机版|小黑屋|Silian Lighting+ ( 蜀ICP备14004521号-1 )

    GMT+8, 2024-4-26 03:44 , Processed in 1.078125 second(s), 22 queries .

    Powered by Discuz! X3.2

    © 2001-2013 Comsenz Inc.

    快速回复 返回顶部 返回列表