scipy.interpolate 单变量插值

作者: ClownFreeMan | 来源:发表于2021-04-13 14:24 被阅读0次

前言

1　数据预处理阶段，处理缺失值的方法分为三类，删除记录　数据插值　不处理，scipy提供了插值法的一些库接口，支持一维和多维（单变量和多元）插值类，Lagrange和Taylor多项式插值器等。本文解释一维插值法的几种python接口和使用案例（单变量插值）

２　一维插值的数据生成，主要有几个方式，生成多项式，用多项式预测未来值或者中间的值，使用导数／积分等数学上能让插值线性平滑(考虑多阶导数的平滑)的操作，来生成插值。几种出名的插值　拉格朗日插值、牛顿插值、分段线性插值、分段三次Hermite插值和样条插值（三次）

３　另外可以利用分段多项式，方便的给一维函数关系求n阶导数、n阶积分，定积分，求解方程

4 其它工具，　lagrange插值、approximate_taylor插值、pade插值

参考：https://docs.scipy.org/doc/scipy/reference/interpolate.html
https://blog.csdn.net/ddjhpxs/article/details/105655394

function name	function inturation
interp1d(x, y[, kind, axis, copy, … ])	Interpolate a 1-D function.
BarycentricInterpolator(xi[, yi, axis])	The interpolating polynomial for a set of points
KroghInterpolator(xi, yi[, axis])	Interpolating polynomial for a set of points.
barycentric_interpolate(xi, yi, x[, axis])	Convenience function for polynomial interpolation.
krogh_interpolate(xi, yi, x[, der, axis ])	Convenience function for polynomial interpolation.
pchip_interpolate(xi, yi, x[, der, axis ])	Convenience function for pchip interpolation.
CubicHermiteSpline(x, y, dydx[, axis, … ])	Piecewise-cubic interpolator matching values and first derivatives.
PchipInterpolator(x, y[, axis, extrapolate ])	PCHIP 1-D monotonic cubic interpolation.
Akima1DInterpolator(x, y[, axis ])	Akima interpolator
CubicSpline(x, y[, axis, bc_type, extrapolate ])	Cubic spline data interpolator.
PPoly(c, x[, extrapolate, axis ])	Piecewise polynomial in terms of coefficients and breakpoints
BPoly(c, x[, extrapolate, axis ])	Piecewise polynomial in terms of coefficients and breakpoints.

interp1d

功能描述：一维插值（分段线性插值），给定的　x １d-array 和　y nd-array，　生成一个近似的映射函数f:y = f(x), 可以使用　ynew = f(xnew) 来生成插值数据.
参数定义：scipy.interpolate.interp1d(x, y, kind='linear', axis=- 1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)
参数说明：
- kind: 指定插值函数的数据生成线性特征　‘linear’ ‘nearest’‘nearest-up’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’, ‘previous’, ‘next’

nearest   "snaps" to the nearest data point.
zero      is a zero order spline. It's value at any point is the last raw value seen.
linear    performs linear interpolation and slinear uses a first order spline. They use different code and can produce similar but subtly different results.
quadratic uses second order spline interpolation.
cubic     uses third order spline interpolation.


import numpy as np
import matplotlib.pyplot as plt
import scipy.interpolate as interpolate

np.random.seed(6)
kinds = ('nearest', 'zero', 'linear', 'slinear', 'quadratic', 'cubic')

N = 10
x = np.linspace(0, 1, N)
y = np.random.randint(10, size=(N,))

new_x = np.linspace(0, 1, 28)
fig, axs = plt.subplots(nrows=len(kinds)+1, sharex=True)
axs[0].plot(x, y, 'bo-')
axs[0].set_title('raw')
for ax, kind in zip(axs[1:], kinds):
    new_y = interpolate.interp1d(x, y, kind=kind)(new_x)
    ax.plot(new_x, new_y, 'ro-')
    ax.set_title(kind)

plt.show()

axis: 指定ｙ轴哪个维度插值，默认是y的最后一维上
copy: 插值是引用还是拷贝
bounds_error:　如果为True，则任何时候尝试对x范围之外的值进行插值都会引发ValueError（在此情况下需要进行插值）。如果为False，则分配超出范围的值fill_value。默认情况下，除非引发错误fill_value="extrapolate"
fill_value:　用来指定填充值
assume_sorted: 如果为False，则x的值可以按任何顺序排列，并且将首先对其进行排序。如果为True，则x必须是单调递增值的数组。
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/interpolate.py#L329-L714

试例:

import matplotlib.pyplot as plt
from scipy import interpolate
import numpy as np

x = np.arange(0, 10)
y = np.exp(-x/3.0)  #　e的幂指数
f = interpolate.interp1d(x, y)

xnew = np.arange(0, 9, 0.1)
ynew = f(xnew)   # use interpolation function returned by `interp1d`
plt.plot(x, y, 'o', xnew, ynew, '-')
plt.show()

x
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
y
array([1.        , 0.71653131, 0.51341712, 0.36787944, 0.26359714,
       0.1888756 , 0.13533528, 0.09697197, 0.06948345, 0.04978707])
xnew
array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1.1, 1.2,
       1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1, 2.2, 2.3, 2.4, 2.5,
       2.6, 2.7, 2.8, 2.9, 3. , 3.1, 3.2, 3.3, 3.4, 3.5, 3.6, 3.7, 3.8,
       3.9, 4. , 4.1, 4.2, 4.3, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9, 5. , 5.1,
       5.2, 5.3, 5.4, 5.5, 5.6, 5.7, 5.8, 5.9, 6. , 6.1, 6.2, 6.3, 6.4,
       6.5, 6.6, 6.7, 6.8, 6.9, 7. , 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7,
       7.8, 7.9, 8. , 8.1, 8.2, 8.3, 8.4, 8.5, 8.6, 8.7, 8.8, 8.9]
ynew
array([1.        , 0.97165313, 0.94330626, 0.91495939, 0.88661252,
       0.85826566, 0.82991879, 0.80157192, 0.77322505, 0.74487818,
       0.71653131, 0.69621989, 0.67590847, 0.65559705, 0.63528563,
       0.61497421, 0.5946628 , 0.57435138, 0.55403996, 0.53372854,
       0.51341712, 0.49886335, 0.48430958, 0.46975582, 0.45520205,
       0.44064828, 0.42609451, 0.41154074, 0.39698698, 0.38243321,
       0.36787944, 0.35745121, 0.34702298, 0.33659475, 0.32616652,
       0.31573829, 0.30531006, 0.29488183, 0.2844536 , 0.27402537,
       0.26359714, 0.25612498, 0.24865283, 0.24118068, 0.23370852,
       0.22623637, 0.21876422, 0.21129206, 0.20381991, 0.19634776,
       0.1888756 , 0.18352157, 0.17816754, 0.17281351, 0.16745947,
       0.16210544, 0.15675141, 0.15139738, 0.14604335, 0.14068932,
       0.13533528, 0.13149895, 0.12766262, 0.12382629, 0.11998996,
       0.11615363, 0.11231729, 0.10848096, 0.10464463, 0.1008083 ,
       0.09697197, 0.09422312, 0.09147426, 0.08872541, 0.08597656,
       0.08322771, 0.08047886, 0.07773001, 0.07498115, 0.0722323 ,
       0.06948345, 0.06751381, 0.06554417, 0.06357454, 0.0616049 ,
       0.05963526, 0.05766562, 0.05569598, 0.05372634, 0.05175671])

BarycentricInterpolator

功能描述：重心多项式坐标插值，给定的　x １d-array 和　y nd-array，　构造多项式.
参数定义：scipy.interpolate.BarycentricInterpolator(xi, yi=None, axis=0)
参数说明：
- xi: 多项式应通过的点的x坐标的一维数组
- yi: 多项式应通过的点的y坐标。如果为None，则稍后将通过set_y方法提供y值。
- axis: y中的维度
- add_xi(self, xi[, yi]) Add more x values to the set to be interpolated
- set_yi(self, yi[, axis]) Update the y values to be interpolated
- call(self, x):
参考: https://blog.csdn.net/silangquan/article/details/21990713
https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/polyint.py#L497-L650
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/polyint.py#L497-L650

试例:

import matplotlib.pyplot as plt
from scipy import interpolate
import numpy as np

x = np.array([1,2,3,4,5,6,7])
x = np.array([1,1.2,1.3,1.4,5,6,7])
y = np.array([1,3,5,7,4,2,1])

x = np.arange(0, 10)
y = np.exp(-x/3.0)  #　e的幂指数

bci = interpolate.BarycentricInterpolator(x, y)   # create Barycentric interpolation object

x1 = np.arange(0.5, 9, 1)
y1 =  bci(x1)                                       # __call__(x) 
x2 = np.arange(0.9, 9, 1)
y2 =  bci(x2)

plt.plot(x, y, 'o', x1, y1, 'o', x2, y2, 'o') 
plt.show()

barycentric_interpolate

功能描述：函数重心插值法，给定的　x １d-array 和　y nd-array，　和　BarycentricInterpolator实现原理一致，　由于重心插值计算相对较慢，若多次计算插值，应使用BarycentricInterpolator。
参数定义：scipy.interpolate.barycentric_interpolate(xi, yi, x, axis=0)
参数说明：
- xi: 多项式应通过的点的x坐标的一维数组
- yi: 多项式应通过的点的y坐标。
- x: 要估计的点
- axis: y中的维度
- y: @return: 估计值
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/polyint.py#L653-L714

试例:

import matplotlib.pyplot as plt
from scipy.interpolate import barycentric_interpolate
import numpy as np

x_observed = np.linspace(0.0, 10.0, 11)
y_observed = np.sin(x_observed)
x = np.linspace(min(x_observed), max(x_observed), num=100)
y = barycentric_interpolate(x_observed, y_observed, x)
plt.plot(x_observed, y_observed, "o", label="observation")
plt.plot(x, y, label="barycentric interpolation")
plt.legend()
plt.show()

KroghInterpolator

功能描述：估值导数的插值法(克罗格插值)，给定的　x １d-array 和　y nd-array，　构造多项式.可估计某个点处的导数
参数定义：scipy.interpolate.KroghInterpolator(xi, yi, axis=0)
参数说明：
- xi: 多项式应通过的点的x坐标的一维数组
- yi: 多项式应通过的点的y坐标。
- axis: y中的维度
- call(self, x) 估值
- derivative(self, x[, der]) 估计一个点处的导数
- derivatives(self, x[, der])　估计多个点处的导数
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/polyint.py#L220-L354

试例:

import matplotlib.pyplot as plt
from scipy import interpolate
import numpy as np

x = np.arange(0, 10)
y = np.exp(-x/3.0)  #　e的幂指数
# y = x / 3.0
ki = interpolate.KroghInterpolator(x, y)

x1 = np.arange(0.5, 9, 1)
y1 =  ki(x1)                                       # __call__(x) 
x2 = np.arange(0.9, 9, 1)
y2 =  ki(x2)

derivatives = ki.derivative(x1)
derivatives

plt.plot(x, y, 'o', x1, y1, 'o', x2, y2, 'o') 
plt.show()

krogh_interpolate

功能描述：函数　估值导数的插值法，给定的　x １d-array 和　y nd-array，　构造多项式.可估计某个点处的导数，　实现原理和KroghInterpolator一致
参数定义：scipy.interpolate.krogh_interpolate(xi, yi, x, der=0, axis=0)
参数说明：
- xi: 多项式应通过的点的x坐标的一维数组
- yi: 多项式应通过的点的y坐标。
- x: 要计算的一个或多个点
- der: 估值导数的阶数　？？
- axis: y中的维度
- d: @return: 估计的插值
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/polyint.py#L357-L420

试例:

import matplotlib.pyplot as plt
from scipy.interpolate import krogh_interpolate
import numpy as np

x_observed = np.linspace(0.0, 10.0, 11)
y_observed = np.sin(x_observed)
x = np.linspace(min(x_observed), max(x_observed), num=100)
y = krogh_interpolate(x_observed, y_observed, x, [1,2,3,4,5])
plt.plot(x_observed, y_observed, "o", label="observation")
plt.plot(x, y[0], label="krogh derivative 1")
plt.plot(x, y[1], label="krogh derivative 2")
plt.plot(x, y[2], label="krogh derivative 3")
plt.plot(x, y[3], label="krogh derivative 4")
plt.plot(x, y[4], label="krogh derivative 5")
plt.legend()
plt.show()

PchipInterpolator

功能描述：分段三次埃米尔特插值法，给定的　x １d-array 和　y nd-array，
参数定义：scipy.interpolate.PchipInterpolator(x, y, axis=0, extrapolate=None)
参考: https://zhuanlan.zhihu.com/p/169664296
参数说明：
- xi: 一维数组，单调递增实数值。x不能包含重复值
- yi: 一维实数数组。y沿插值轴的长度必须等于的长度x。如果是ND数组，用axis 参数选择正确的维度
- axis: y中的维度
- extrapolate: {bool, ‘periodic’, None},　默认是True ,false决定是否根据两个边缘点来延伸,　'periodic' extrapolation 周期延伸
- call(self, x[, nu, extrapolate]) Evaluate the piecewise polynomial or its derivative.
- 估值这个分段多项式／或他的导数
- derivative(self[, nu]) Construct a new piecewise polynomial representing the derivative.
- ```
 构造导数的分段多项式 返回一个　PPoly对象
```
- antiderivative(self[, nu]) Construct a new piecewise polynomial representing the antiderivative.
- 构造不定积分的分段多项式　返回一个　PPoly对象
- roots(self[, discontinuity, extrapolate]) Find real roots of the the piecewise polynomial.
- 构造原始的分段多项式
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/_cubic.py#L156-L299

试例

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import PchipInterpolator as PCHIP
#PCHIP stands for Piecewise Cubic Hermite Interpolating Polynomial

x = [1,2,3,4,5,6]
y = [-1,-1,0,1,1,1]
pchip = PCHIP(x,y)
x_new = [2.5,3.5,4.5]
y_new = pchip.__call__(x_new) # y_new=pchip(x_new)
print(y_new)

pp_d = pchip.derivative()           # 导数的ppoly 对象
pp_anti_d = pchip.antiderivative()  # 不定积分的ppoly　对象
print(pchip.roots())                # 多项式的原式子　ndarray　表示

#同时执行两个plot，才能将两个图绘制到一个图上
plt.plot(x_new, y_new, 'ro')#插值点为红色圆点
plt.plot(x, y, 'b-')#原始点为蓝色线条
plt.show()

pchip_interpolate

功能描述：pchip 插值法的便捷调用，给定的　x １d-array 和　y nd-array
参数定义：scipy.interpolate.pchip_interpolate(xi, yi, x, der=0, axis=0)
参数说明：
- xi: 多项式应通过的点的x坐标的一维数组
- yi: 多项式应通过的点的y坐标。
- x: 要计算的一个或多个点
- der: 估值导数的阶数　？？　默认０
- axis: yi中的维度
- y: @return 估计的插值
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/_cubic.py#L302-L360

试例:

import matplotlib.pyplot as plt
from scipy.interpolate import pchip_interpolate
import numpy as np

x_observed = np.linspace(0.0, 10.0, 11)
y_observed = np.sin(x_observed)
x = np.linspace(min(x_observed), max(x_observed), num=100)
y = pchip_interpolate(x_observed, y_observed, x,[0,1,2,3])
plt.plot(x_observed, y_observed, "o", label="observation")
plt.plot(x, y[0], label="pchip der 0")
plt.plot(x, y[1], label="pchip der 1")
plt.plot(x, y[2], label="pchip der 2")
plt.plot(x, y[3], label="pchip der 3")
plt.legend()
plt.show()

Akima1DInterpolator

功能描述：akima一维插值,The interpolation method by Akima uses a continuously differentiable sub-spline built from piecewise cubic polynomials. The resultant curve passes through the given data points and will appear smooth and natural.
特点是插值显得尽可能使拟合曲线平滑自然
参数定义：scipy.interpolate.Akima1DInterpolator(x, y, axis=0)
参数说明：
- xi: 一维数组，单调递增实数值。x不能包含重复值
- yi: 一维实数数组。y沿插值轴的长度必须等于的长度x。如果是ND数组，用axis 参数选择正确的维度
- axis: y中的维度
- call(self, x[, nu, extrapolate]) Evaluate the piecewise polynomial or its derivative.
- 估值这个分段多项式／或他的导数
- derivative(self[, nu]) Construct a new piecewise polynomial representing the derivative.
- ```
 构造导数的分段多项式 返回一个　PPoly对象
```
- antiderivative(self[, nu]) Construct a new piecewise polynomial representing the antiderivative.
- 构造不定积分的分段多项式　返回一个　PPoly对象
- roots(self[, discontinuity, extrapolate]) Find real roots of the the piecewise polynomial.
- 构造原始的分段多项式
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/_cubic.py#L363-L461

试例:

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import Akima1DInterpolator as Akima1D

#x = [1,2,3,4,5,6]
#y = [-1,-1,0,1,1,1]
x = np.arange(0, 10)
y = np.exp(-x/3.0)  #　e的幂指数

akima = Akima1D(x,y)
x_new = [2.5,3.5,4.5]
y_new = akima.__call__(x_new)
print(y_new)

pp_d = akima.derivative()           # 导数的ppoly 对象
pp_anti_d = akima.antiderivative()  # 不定积分的ppoly　对象
print(akima.roots())                # 多项式的原式子　ndarray　表示

#同时执行两个plot，才能将两个图绘制到一个图上
plt.plot(x_new, y_new, 'ro')#插值点为红色圆点
plt.plot(x, y, 'b-')#原始点为蓝色线条
plt.show()

CubicHermiteSpline

功能描述：Piecewise-cubic interpolator matching values and first derivatives.分段三次埃尔米特插值
特点是可以指定每个观测点的一阶导数值，可以计算拟合多项式，积分多项式，导数多项式，可以求定积分
参数定义：scipy.interpolate.CubicHermiteSpline(x, y, dydx, axis=0, extrapolate=None)
参考：https://github.com/Galdeano/CubicHermiteSpline
　　　https://en.wikipedia.org/wiki/Cubic_Hermite_spline
参数说明：
- x: 多项式应通过的点的x坐标的一维数组
- y: 多项式应通过的点的y坐标。长度和ｘ一致
- dydx: 每个点的导数，　长度和ｘ一致
- axis: y中的维度
- extrapolate: {bool, ‘periodic’, None},　默认是True ,false决定是否根据两个边缘点来延伸,　'periodic' extrapolation 周期延伸
- call(self, x) 计算分段多项式的插值估值
- derivative(self[, nu]) Construct a new piecewise polynomial representing the derivative.
- ```
 构造导数的分段多项式 返回一个　PPoly对象
```
- antiderivative(self[, nu]) Construct a new piecewise polynomial representing the antiderivative.
- 构造不定积分的分段多项式　返回一个　PPoly对象
- integrate(self, a, b[, extrapolate]) Compute a definite integral over a piecewise polynomial.
- ```
 计算一个分段多项式的定积分　array-like
```
- roots(self[, discontinuity, extrapolate]) Find real roots of the the piecewise polynomial.
- 构造原始的分段多项式 ndarray
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/_cubic.py#L70-L153

试例:

import matplotlib.pyplot as plt
import numpy as np
import scipy.interpolate 

x_observed = np.linspace(0.0, 10.0, 11)
y_observed = np.sin(x_observed)
dydx = np.linspace(5,5,11)    ＃　指定观测点的斜率、导数
# dydx = np.linspace(0,0,11)
chs = interpolate.CubicHermiteSpline(x_observed, y_observed, dydx, True, 'periodic') #周期外推

x = np.linspace(-3, 15, num=100)
y = chs.__call__(x)

plt.plot(x_observed, y_observed, "o", label="observation")
plt.plot(x, y, label="CubicHermiteSpline 0")

plt.legend()
plt.show()

pp_d = pchip.derivative()           # 导数的ppoly 对象
pp_anti_d = pchip.antiderivative()  # 不定积分的ppoly　对象
print(chs.integrate(3, 5))          # 计算定积分 
print(chs.roots())                  # 多项式的原式子　ndarray　表示

CubicSpline

功能描述：三次样条数据插值器，给定的　x １d-array 和　y nd-array，　构造多项式.可估计某个点处的导数
参数定义：scipy.interpolate.CubicSpline(x, y, axis=0, bc_type='not-a-knot', extrapolate=None)
参数说明：
- x: 多项式应通过的点的x坐标的一维数组
- y: 多项式应通过的点的y坐标。
- axis: y中的维度
- bc_type: 边界条件类型。需要两个附加的方程式(由边界条件给出)来确定每个线段上多项式的所有系数
- “ not-a-knot”(默认)：曲线末端的第一段和第二段是相同的多项式。当没有边界条件的信息时，这是一个很好的默认设置。
  ‘periodic’：假定插值函数是周期的周期性x[-1] - x[0]。 y的第一个和最后一个值必须相同：y[0] == y[-1]。这种边界条件将导致y'[0] == y'[-1]和y''[0] == y''[-1]。
  ‘clamped’：曲线末端的一阶导数为零。假设y为1Dbc_type=((1, 0.0), (1, 0.0))是相同的条件。
  ‘natural’：曲线末端的二阶导数为零。假设y为1Dbc_type=((2, 0.0), (2, 0.0))是相同的条件。
  如果bc_type是2元组，则第一个和第二个值将分别应用于曲线的起点和终点。元组值可以是前面提到的字符串之一(‘periodic’除外)，也可以是元组(顺序为deriv_values)，允许在曲线末端指定任意导数：
  ‘order’：导数阶1或2。
  ’deriv_value‘：数组包含微分值，形状必须与y相同，但不包括axis尺寸。例如，如果y是一维，那么deriv_value必须是标量。如果y是形状为(n0，n1，n2)且轴= 2的3D，则deriv_value必须是2D并具有形状(n0，n1)。
extrapolate: {bool, ‘periodic’, None},　默认是True ,false决定是否根据两个边缘点来延伸,　'periodic' extrapolation 周期延伸
call(self, x) 计算分段多项式的插值估值
derivative(self[, nu]) Construct a new piecewise polynomial representing the derivative.

 构造导数的分段多项式 返回一个　PPoly对象

antiderivative(self[, nu]) Construct a new piecewise polynomial representing the antiderivative.
构造不定积分的分段多项式　返回一个　PPoly对象
integrate(self, a, b[, extrapolate]) Compute a definite integral over a piecewise polynomial.

 计算一个分段多项式的定积分　array-like

roots(self[, discontinuity, extrapolate]) Find real roots of the the piecewise polynomial.
构造原始的分段多项式 ndarray
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/_cubic.py#L464-L855

试例:

from scipy.interpolate import CubicSpline
import numpy as np
import matplotlib.pyplot as plt

x = np.arange(10)
y = np.sin(x)
cs = CubicSpline(x, y)
xs = np.arange(-0.5, 9.6, 0.1)

fig, ax = plt.subplots(figsize=(6.5, 4))
ax.plot(x, y, 'o', label='data')
ax.plot(xs, np.sin(xs), label='true')
ax.plot(xs, cs(xs), label="S")        #  0 Order of derivative to evaluate
ax.plot(xs, cs(xs, 1), label="S'")    #  1 Order of derivative to evaluate
ax.plot(xs, cs(xs, 2), label="S''")   #  2 Order of derivative to evaluate
ax.plot(xs, cs(xs, 3), label="S'''")
ax.set_xlim(-0.5, 9.5)
ax.legend(loc='lower left', ncol=2)
plt.show()

from scipy.interpolate import CubicSpline
import numpy as np
import matplotlib.pyplot as plt

theta = 2 * np.pi * np.linspace(0, 1, 5)
y = np.c_[np.cos(theta), np.sin(theta)]
cs = CubicSpline(theta, y, bc_type='periodic')
print("ds/dx={:.1f} ds/dy={:.1f}".format(cs(0, 1)[0], cs(0, 1)[1]))

xs = 2 * np.pi * np.linspace(0, 1, 100)
fig, ax = plt.subplots(figsize=(6.5, 4))
ax.plot(y[:, 0], y[:, 1], 'o', label='data')
ax.plot(np.cos(xs), np.sin(xs), label='true')
ax.plot(cs(xs)[:, 0], cs(xs)[:, 1], label='spline')
ax.axes.set_aspect('equal')
ax.legend(loc='center')
plt.show()

PPoly

功能描述：Piecewise polynomial in terms of coefficients and breakpoints　基于系数和断点的分段多项式.
分段多项式的类定义, 对象可以表示一个分段多项式，通常可用于表示　CubicSpline CubicHermiteSpline PchipInterpolator 等插值多项式的一阶导数或者不定积分.
一般不会直接创建对象，　会根据观测值的插值算法derivative、antiderivative来构造对象
参数定义：scipy.interpolate.PPoly(c, x, extrapolate=None, axis=0)
参数说明：
- c: 多项式系数, order k and m intervals.
- x: 多项式分段点，　一维数组，单调递增/递减实数值
- axis: y中的维度
- extrapolate: {bool, ‘periodic’, None},　默认是True ,false决定是否根据两个边缘点来延伸,　'periodic' extrapolation 周期延伸
- call(self, x[, nu, extrapolate])　Evaluate the piecewise polynomial or its derivative.
  计算多项式在　ｘ处的值
- derivative(self[, nu])　Construct a new piecewise polynomial representing the derivative.
  构造导函数的多项式　返回PPoly 对象
- antiderivative(self[, nu])　Construct a new piecewise polynomial representing the antiderivative.
  构造不定积分的多项式　　返回PPoly　对象
- integrate(self, a, b[, extrapolate])　Compute a definite integral over a piecewise polynomial.
- solve(self[, y, discontinuity, extrapolate])　Find real solutions of the the equation pp(x) == y.
  已知y,反解　pp(x) == y
- roots(self[, discontinuity, extrapolate])　Find real roots of the the piecewise polynomial.
  返回分段多项式的系数数组
- extend(self, c, x[, right])　Add additional breakpoints and coefficients to the polynomial.
  在分段多项式上增加分段点和系数
- from_spline(tck[, extrapolate])　Construct a piecewise polynomial from a spline
- from_bernstein_basis(bp[, extrapolate])　Construct a piecewise polynomial in the power basis from a polynomial in Bernstein basis.
- construct_fast(c, x[, extrapolate, axis])　Construct the piecewise polynomial without making checks.
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/interpolate.py#L942-L1348

试例:

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import PchipInterpolator as PCHIP

#x = [1,2,3,4,5,6]
#y = [-1,-1,0,1,1,1]
x = np.arange(0, 10)
y = np.exp(-x/3.0)  #　e的幂指数

pchip = PCHIP(x,y)
x_new = [2.5,3.5,4.5]
y_new = pchip.__call__(x_new) # y_new=pchip(x_new)
print(y_new)

pp_x = np.arange(-2, 12, 0.3)

pp_d = pchip.derivative()           # 导数的ppoly 对象
pp_antid = pchip.antiderivative()  # 不定积分的ppoly　对象
y_ppd = pp_d(pp_x)
y_ppantid = pp_antid(pp_x)

#同时执行两个plot，才能将两个图绘制到一个图上
plt.plot(x_new, y_new, 'ro')#插值点为红色圆点
plt.plot(x, y, 'b-')#原始点为蓝色线条
plt.plot(pp_x, y_ppd, 'y-')#原始点为yellow线条
plt.plot(pp_x, y_ppantid, 'r-')#原始点为red线条
plt.show()

BPoly

功能描述：(Bernstein polynomial 伯恩斯坦多项式 )Piecewise polynomial in terms of coefficients and breakpoints　基于系数和断点的分段多项式.
参数定义：scipy.interpolate.BPoly(c, x, extrapolate=None, axis=0)
参考：https://en.wikipedia.org/wiki/Bernstein_polynomial
http://www.idav.ucdavis.edu/education/CAGDNotes/Bernstein-Polynomials.pdf

The polynomial between x[i] and x[i + 1] is written in the Bernstein polynomial basis:

S = sum(c[a, i] * b(a, k; x) for a in range(k+1)),
where k is the degree of the polynomial, and:

b(a, k; x) = binom(k, a) * t**a * (1 - t)**(k - a),
with t = (x - x[i]) / (x[i+1] - x[i]) and binom is the binomial coefficient.

参数说明：
- c: 多项式系数, order k and m intervals.
- x: 多项式分段点，　一维数组，单调递增/递减实数值
- axis: y中的维度
- extrapolate: {bool, ‘periodic’, None},　默认是True ,false决定是否根据两个边缘点来延伸,　'periodic' extrapolation 周期延伸
- call(self, x[, nu, extrapolate])　Evaluate the piecewise polynomial or its derivative.
  计算多项式在　ｘ处的值
- extend(self, c, x[, right]) Add additional breakpoints and coefficients to the polynomial.
  添加断点和系数
- derivative(self[, nu])　Construct a new piecewise polynomial representing the derivative.
  构造导函数的多项式　返回 BPoly 对象
- antiderivative(self[, nu])　Construct a new piecewise polynomial representing the antiderivative.
  构造不定积分的多项式　　返回 BPoly　对象
- integrate(self, a, b[, extrapolate]) Compute a definite integral over a piecewise polynomial.
  计算定积分
- construct_fast(c, x[, extrapolate, axis]) Construct the piecewise polynomial without making checks.
  工厂方法，　与默认构造方法区别是不做检测
- from_power_basis(pp[, extrapolate]) Construct a piecewise polynomial in Bernstein basis from a power basis polynomial.
  构造 PPoly　对象
- from_derivatives(xi, yi[, orders, extrapolate]) Construct a piecewise polynomial in the Bernstein basis, compatible with the specified values and derivatives at breakpoints.
  工厂方法，用导数系数构造　BPoly　对象
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/interpolate.py#L1351-L1902

试例:

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import BPoly

x = [0, 1]
c = [[1], [3], [1]]
bp = BPoly(c, x)

bp_x = np.arange(-3, 10, 0.1)
bp_y = bp(bp_x)

plt.plot(bp_x, bp_y, 'r-')#原始点为red线条
plt.show()

拉格朗日插值

功能描述：拉格朗日插值给定两个一维数组x和w，返回通过点的拉格朗日插值多项式.
参数定义：scipy.interpolate.lagrange(x, w)
参考: https://numpy.org/devdocs/reference/generated/numpy.poly1d.html#numpy.poly1d
参数说明：
- x: 代表一组数据点的x坐标.
- w: w代表一组数据点的y坐标
- return: numpy.poly1d实例
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/interpolate.py#L25-L82

试例:

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import lagrange

x = np.array([0, 1, 2])
y = x**3
poly = lagrange(x, y)

lagrange_x = np.arange(-10, 10, 0.1)
lagrange_y = poly(lagrange_x)

plt.plot(lagrange_x, lagrange_y, 'r-')#原始点为red线条
plt.show()

泰勒多项式插值

功能描述：近似泰勒多项式插值给定两个一维数组x和w，返回通过点的拉格朗日插值多项式.
参数定义：scipy.interpolate.approximate_taylor_polynomial(f, x, degree, scale, order=None)
参考: https://numpy.org/devdocs/reference/generated/numpy.poly1d.html#numpy.poly1d
参数说明：
- f: 寻求泰勒多项式的函数。应该接受x值的向量.
- x: 求多项式的点。
- degree: 泰勒多项式的次数
- scale：　拟合中使用的多项式的阶数；f将被评估order+1时间
- return: numpy.poly1d实例
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/polyint.py#L423-L494

试例:

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import approximate_taylor_polynomial

x = np.linspace(-10.0, 10.0, num=100)
plt.plot(x, np.sin(x), label="sin curve")

for degree in np.arange(1, 15, step=2):
    sin_taylor = approximate_taylor_polynomial(np.sin, 0, degree, 1,order=degree + 2)
    plt.plot(x, sin_taylor(x), label="degree={0}".format(degree))

plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left',borderaxespad=0.0, shadow=True)
plt.tight_layout()
plt.axis([-10, 10, -10, 10])
plt.show()

Pade

功能描述：Pade近似值作为两个多项式之比返回多项式.
参数定义：scipy.interpolate.pade(an, m, n=None)
参考: https://numpy.org/devdocs/reference/generated/numpy.poly1d.html#numpy.poly1d
参数说明：
- an: 泰勒级数系数.
- m: 近似多项式的阶数
- p q: @turn 所限定的多项式的Pade逼近的是 p(x)/q(x)
源码:

https://github.com/scipy/scipy/blob/v1.6.2/scipy/interpolate/_pade.py#L6-L65

试例:

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import pade

e_exp = [1.0, 1.0, 1.0/2.0, 1.0/6.0, 1.0/24.0, 1.0/120.0]
p, q = pade(e_exp, 2)

e_exp.reverse()
e_poly = np.poly1d(e_exp)

e_poly(1)
p(1)/q(1)

scipy.interpolate 单变量插值

前言

目录

interp1d

BarycentricInterpolator

barycentric_interpolate

KroghInterpolator

krogh_interpolate

PchipInterpolator

pchip_interpolate

Akima1DInterpolator

CubicHermiteSpline

CubicSpline

PPoly

BPoly

拉格朗日插值

泰勒多项式插值

Pade

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读