scipy--插值

1 scipy.interpolate

scipy.interpolate是插值模块，插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。与拟合不同的是，要求曲线通过所有的已知数据。计算插值有两种基本的方法：

• 对一个完整的数据集去拟合一个函数；
• 仿样内插法：对数据集的不同部分拟合出不同的函数，而函数之间的曲线平滑对接。

SciPy的interpolate模块提供了许多对数据进行插值运算的函数，范围涵盖简单的一维插值到复杂多维插值求解。

当样本数据变化归因于一个独立的变量时，就使用一维插值；反之样本数据归因于多个独立变量时，使用多维插值。

2 一维插值

2.1 内插值 interp1d()

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

其中，x和y参数是一系列已知的数据点，kind参数是插值类型，可以是字符串或整数，它给出插值的B样条曲线的阶数，候选值及作用下表所示：

候选值	作用
‘zero’ 、‘nearest’	阶梯插值，相当于0阶B样条曲线
‘slinear’ 、‘linear’	线性插值，用一条直线连接所有的取样点，相当于一阶B样条曲线
‘quadratic’ 、‘cubic’	二阶和三阶B样条曲线，更高阶的曲线可以直接使用整数值指定

interp1d比Matlab的interp有些优势，因为返回的是函数，不需要在事先设定需要求解的点，而是在需要使用时调用函数。

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

x=np.linspace(0,10,5)
y=np.sin(x)
xnew=np.linspace(0,10,101)
plt.plot(x,y,'ro', label="interpolating point")
list1=['linear','nearest']

for kind in list1:
    f=interpolate.interp1d(x,y,kind=kind)
    #f是一个函数，用这个函数就可以找插值点的函数值了：
    ynew=f(xnew)
    plt.plot(xnew,ynew,label=kind)
plt.plot(xnew,np.sin(xnew),label="real f(x)",linestyle="--")
plt.legend(loc='lower right')
plt.show()

out:

但是！interp1d不能外推运算（外插值）

2.2 外插值

我们可以使用UnivariateSpline进行外插值
调用方式如下：
UnivariateSpline(x,y,w=None,bbox=[None,None],k=3,s=None)

• x,y是X-Y坐标数组
• w是每个数据点的权重值
• k为样条曲线的阶数
• s为平滑参数。
• • s=0，样条曲线强制通过所有插值点
• • s>0,满足$\sum (w(y-spline(x)){)}^2\le s$

x1=np.linspace(0,10,20)
y1=np.sin(x1)

sx1=np.linspace(0,15,100)
func1=interpolate.UnivariateSpline(x1,y1,s=0)
sy1=func1(sx1)
plt.plot(x1,y1,'o', label="interpolation points")
plt.plot(sx1,sy1, label="interpolation function")
plt.plot(sx1,np.sin(sx1), label="real function")
plt.legend()
plt.show()

out:
我们可以发现，在插值区间 [0,10] 之间，interpolation function 和 real function极为接近，但是一旦超出该区间，即外插值，插值的效果就会非常差。

3 二维插值

2.1 interp2d()

interp2d(x,y,z,kind='linear')
这里有几个注意事项：

1. interp2d()中，输入的x,y,z先用ravel()被转成了一维数组
2. func()的输入必须是一维的，输出是二维的（有点奇怪，感觉完成度不高）
3. 插值的源数据必须是等距网格。不然的话，运行不保存但结果不对。

def func(x,y):
    # return (x+y)*np.exp(-5*(x**2+y**2))
    return np.sin(pow(x,2)+pow(y,2))
x,y=np.mgrid[-2:2:8j,-2:2:8j] # 等距网络体现在这里
z=func(x,y)
choose = np.random.rand(8,8)
f = interpolate.interp2d(x,y,z,kind='cubic') # note x and y are 2-dim
x_real,y_real=np.mgrid[-2:2:1000j,-2:2:1000j]
z_real = func(x_real,y_real)
xnew=np.linspace(-2,2,100)
ynew=np.linspace(-2,2,100)
znew=f(xnew,ynew)# xnew, ynew是一维的，输出znew是二维的
xnew,ynew=np.mgrid[-2:2:100j,-2:2:100j]#统一变成二维，便于下一步画图

ax=plt.subplot(111,projection='3d')
ax.plot_surface(xnew,ynew,znew, color="orange")
ax.plot_surface(x_real,y_real,z_real,color="g")
ax.scatter(x,y,z,c='r',marker='^')
plt.show()

out：

Rbf()

def func(x,y):
    # return (x+y)*np.exp(-5*(x**2+y**2))
    return np.sin(1/x+1/y)
x,y = np.random.randn(50), np.random.randn(50)
z=func(x,y)
choose = np.random.rand(8,8)
f = interpolate.Rbf(x,y,z)

x_real,y_real = np.mgrid[-2:2:1000j,-2:2:1000j]
z_real = func(x_real,y_real)


xnew, ynew = np.mgrid[-2:2:1000j, -2:2:1000j]
znew = f(xnew,ynew) # 此时 xnew, ynew是二维的，输出znew是二维的

ax = plt.subplot(111,projection='3d')
ax.plot_surface(xnew,ynew,znew, color="orange")
ax.plot_surface(x_real,y_real,z_real,color="g")
ax.scatter(x,y,z,c='r',marker='^')
plt.show()

参考：
https://www.jianshu.com/p/b306095309db
https://www.guofei.site/2017/06/06/scipyinterpolate.html