scipy.interpolate 单变量插值

2021-04-13  本文已影响0人  ClownFreeMan

前言

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

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()

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

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

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

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

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

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

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

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

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

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

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

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.

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()

拉格朗日插值

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()

泰勒多项式插值

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

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)
上一篇 下一篇

猜你喜欢

热点阅读