Python小白的数学建模课-22.插值方法

Posted youcans

tags:

篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了Python小白的数学建模课-22.插值方法相关的知识，希望对你有一定的参考价值。

Python小白的数学建模课-22.插值方法


  • 插值、拟合、回归和预测，都是数学建模中经常提到的概念，也经常被混淆。
  • 插值，是在离散数据的基础上补插连续函数，使得插值函数通过全部给定的离散数据点，多用于图像处理和缺失数据处理。
  • 使用 Scipy 工具包的 interpolate 插值模块，通过例程讲解一维插值、二维插值的实现方法。
  • 『Python小白的数学建模课 @ Youcans』带你从数模小白成为国赛达人。

1. 数据插值

1.1 插值与拟合

插值与拟合，不仅是基本的数学建模方法，也是最常用的数据处理方法。

不过，插值和拟合却经常会被混为一谈，所以我们首先看看这两个概念。

  • 插值，是在离散数据的基础上补插连续函数，使得插值函数通过全部给定的离散数据点。 插值是离散函数逼近的重要方法，利用它可通过函数在有限个点处的取值状况，估算出函数在其他点处的近似值。简单地说，插值是求过一组已知点的近似函数。
  • 拟合，是用一个连续函数（曲线）靠近给定的离散数据，使其与给定的数据相吻合。拟合也是根据一组已知点求近似函数，但不要求过已知点。

因此，插值和拟合都是根据一组已知数据点，求变化规律和特征相似的近似曲线的过程。但是插值要求近似曲线完全经过所有的给定数据点，而拟合只要求近似曲线在整体上尽可能接近数据点，并反映数据的变化规律和发展趋势。

插值可以看作是一种特殊的拟合，是要求误差函数为 0 的拟合。由于数据点通常都带有误差，误差为 0 往往意味着过度拟合，过拟合模型对于训练集以外的数据的泛化能力往往是较差的。因此在实践中，插值多用于图像处理和缺失数据处理，拟合多用于实验数据处理。

此外，还有一个常用而且容易混淆的概念： 回归。回归是研究一组随机变量与另一组随机变量之间关系的统计分析方法，包括建立数学模型并估计模型参数，并检验数学模型的可信度，也包括利用建立的模型和估计的模型参数进行预测或控制。

回归是一种数据分析方法，拟合是一种具体的数据处理方法。拟合侧重于曲线参数寻优，使曲线与数据相符；而回归侧重于研究两个或多个变量之间的关系。


1.2 Scipy 插值工具箱

数据插值是数据处理的常用方法，常见的插值算法有线性插值、B样条插值、临近插值等。

Scipy 工具包带有插值工具箱，提供了丰富的插值方法和函数，可以用于一维、二维和多维插值。

一维函数插值的类 interp1d，提供了多种插值方法，如样条函数插值、一维和多维插值、拉格朗日插值、泰勒多项式插值及自定义插值函数。函数 griddata 则提供了 N维插值的接口（N=1,2,3,…）。


2. Scipy 一维插值方法：内插值

2.1 一维插值类 interp1d

Scipy.interpolate 中的 interp1d 类是一种基于固定数据点创建函数的方法，可以使用函数插值在给定数据定义的域内的任何位置对其进行计算。注意 interp1d 是内插法，不能外推运算（外插值）。

该类定义调用，允许使用 x 轴值调用对象，此时应计算插值函数，并返回插值的 y 轴值。具体地说，interp1d 类生成已知数据点集的插值函数 y=f(x)，通过调用这个插值函数，可以在已知数据之间插值，得到指定 x 的函数值 f(x)。

class scipy.interpolate.interp1d(x, y, kind=‘linear’, axis=- 1, copy=True, bounds_error=None, fill_value=nan, assume_sorted=False)

主要参数：

  • x：一维数组，给定数据点集的 x 值。
  • y：N 维数组，给定数据点集的 y 值，数组长度必须与 x 相等。
  • kind：字符串或整数，可选项，指定使用的样条曲线的种类或插值方法。
    • 可选的字符串：‘linear’, ‘nearest’, ‘nearest-up’, ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’, ‘previous’, ‘next’；
    • ‘zero’, ‘slinear’, ‘quadratic’, ‘cubic’ 分别表示零次、一次、二次、三次样条插值；
    • ‘previous’, ‘next’ 分别表示只前点插值或后点插值；
    • ‘nearest’ 表示向下舍入， ‘nearest-up’ 表示向上舍入；
    • 默认值为 ‘linear’，即线性插值。

interp1d 允许通过参数 bounds_error、fill_value 设置外推时的边界值，但这并不是进行外推插值。

返回值：

  • 类 interp1d() 返回一个函数，其调用方法使用插值来查找新点的值。

2.2 Python 例程：interp1d 的使用

使用示例：

# 1. 一维插值使用示例
import numpy as np
import matplotlib.pyplot as plt  # 导入 Matplotlib 工具包
from scipy.interpolate import interp1d  # 导入 scipy 中的一维插值工具 interp1d

# 已知数据点集 (x,y)
x = [0.0, 2.0, 4.0, 6.0, 8.0, 10.0]  # 已知数据 x
y = [3.1, 2.7, 1.5, 0.1, 1.0, 3.9]   # 已知数据 y
# 由给定数据点集 (x,y) 求插值函数 fx
fx = interp1d(x, y, kind='linear')  # 由已知数据 (x,y) 求出插值函数 fx
# 由插值函数 fx 计算插值点的函数值
xInterp = np.linspace(0,10,100) # 指定需插值的数据点集 xInterp
yInterp = fx(xInterp)  # 调用插值函数 fx，计算 xInterp 的函数值
# 绘图
plt.plot(xInterp, yInterp, label="linear interpolate")
plt.show()

2.3 Python 例程：一维插值方法比较

通过设置 interp1d 类的参数kind，可以指定使用的样条曲线的种类或插值方法。

上节中介绍了 kind 各种选项所指的插值方法，本节进一步通过例程来比较不同方法的插值结果。

Python 例程：

# mathmodel24_v1.py
# Demo24 of mathematical modLSing algorithm
# Demo of interpolate with Scipy.interpolate
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-01

# 2. 一维插值方法(内插)比较
import numpy as np
import matplotlib.pyplot as plt  # 导入 Matplotlib 工具包
from scipy.interpolate import interp1d  # 导入 scipy 中的一维插值工具 interp1d

# 生成已知数据点集 (x,y)，需插值的数据点集 xnew
np.random.seed(5)
x = np.linspace(0, 5, 10)  # 生成已知数据点集的 x
y = np.cos(x/10)*2 + 0.5*np.random.rand(10)  # 生成已知数据点集的 y
xnew = np.linspace(0, 5, 100)  # 指定需插值的数据点集 xnew

# 使用不同插值方法，由给定数据点集 (x,y) 求插值函数 fx
f1 = interp1d(x, y, kind="linear")  # 线性插值
f2 = interp1d(x, y, kind="zero")  # 零阶样条插值
f3 = interp1d(x, y, kind="slinear")  # 一次样条插值
f4 = interp1d(x, y, kind="quadratic")  # 二次样条插值
f5 = interp1d(x, y, kind="cubic")  # 三次样条插值
f6 = interp1d(x, y, kind="nearest")  # 临近点插值，向下舍入
# f7 = interp1d(x, y, kind="nearest-up")  # 临近点插值，向上舍入
f8 = interp1d(x, y, kind="previous")  # 前点插值
f9 = interp1d(x, y, kind="next")  # 后点插值

# 绘图
plt.figure(figsize=(8,6))
plt.suptitle("Data interpolate")  # 全局标题
plt.subplot(221)
plt.plot(x, y, "o",  label="data")  # 已知数据点
plt.plot(xnew, f2(xnew), label="0-order spline")  # 零阶样条插值
plt.plot(xnew, f3(xnew), label="1-order spline")  # 一阶样条插值
plt.legend(loc="lower left")
plt.subplot(222)
plt.plot(x, y, "o",  label="data")  # 已知数据点
plt.plot(xnew, f4(xnew), label="2-order spline")  # 二阶样条插值
plt.plot(xnew, f5(xnew), label="3-order spline")  # 三阶样条插值
plt.legend(loc="lower left")
plt.subplot(223)
plt.plot(x, y, "o",  label="data")  # 已知数据点
plt.plot(xnew, f1(xnew), label="linear")  # 线性插值
plt.plot(xnew, f6(xnew), label="nearest")  # 临近点插值，向下舍入
# plt.plot(xnew, f7(xnew), label="nearest-up")  # 临近点插值，向上舍入
plt.legend(loc="lower left")
plt.subplot(224)
plt.plot(x, y, "o",  label="data")  # 已知数据点
plt.plot(xnew, f8(xnew), label="previous")  # 前点插值
plt.plot(xnew, f9(xnew), label="next")  # 后点插值
plt.legend(loc="lower left")
plt.show()

程序运行结果：

结果分析：

  1. 线性插值是常用的插值方法，简单地说可以理解为将相邻的数据点用线段连接，算法简单、运算速度快。一阶样条曲线插值，等效于线性插值。

  2. 最近邻点插值、前点插值、后点插值和 0阶样条插值的结果都是阶梯形状曲线，只是选点方法不同。

  3. 样条插值是重要的插值方法，用光滑曲线连接数据点。每一个样条都是用一个多项式表达的，多项式的次数就是样条曲线的阶数，决定了样条曲线的形状和性质：

    • 0 阶样条曲线，在每一区间上样条函数为常数，样条曲线呈阶梯形状；
    • 1 阶样条曲线，在每一区间上样条函数为线性函数，样条曲线呈折线段形状；
    • 2 阶样条曲线，在每一区间上样条函数为二次函数，整体一阶连续可导；依次类推。
    • 由二阶开始，插值函数不再具有局域性，改变某一节点，函数整体都会改变。
    • 2 阶和 3阶样条插值最为常用，更高阶的样条插值过于复杂，通常结果的差异并不大。图中的 2阶和 3阶样条插值结果就已经近似是重合的。

3. Scipy 一维插值方法：外插值

3.1 一维插值类 UnivariateSpline

Scipy.interpolate 中的 UnivariateSpline 类是一种基于固定数据点创建函数的方法，使用一维样条曲线拟合到给定的数据点集。

该类定义调用，允许使用 x 轴值调用对象，此时应计算样条曲线，并返回插值的 y 轴值。具体地说，UnivariateSpline 类生成已知数据点集的样条插值函数 y=spl(x)，通过调用样条插值函数，可以计算指定 x 的函数值 f(x)。

class scipy.interpolate.UnivariateSpline(x, y, w=None, bbox=[None, None], k=3, s=None, ext=0, check_finite=False)

主要参数：

  • x：一维数组，数值必须递增。
  • y：一维数组，数组长度必须与 x 相等。
  • w：一维数组，正数，可选项。每个数据点的权重，默认所有点的权重相等。
  • k：整数，可选项。样条函数的阶数， 1 ≤ k ≤ 5 1 \\leq k \\leq 5 1≤k≤5，默认值为 3。
  • s：实数，可选项，平滑参数：
    • s=0，数据插值，样条曲线必须通过所有数据点；
    • s>0，数据拟合，满足 ∑ ( w ( y − s p l ( x ) ) ) 2 ≤ s \\sum(w(y-spl(x)))^2 \\leq s ∑(w(y−spl(x)))2≤s；
    • 默认不设置 s，则 s=len(w)。
  • ext：整数或字符串，可选项。用于控制外推插值的方案：
    • ext= 0 或 “extrapolate”，返回外推值，默认值；
    • ext= 1 或 “zeros”，返回 0；
    • ext= 2 或 “raise”，抛出异常值 ValueError；
    • ext= 3 或 “const”，返回边界值。

返回值：

  • 类 UnivariateSpline() 返回一个函数，其调用方法使用插值来查找新点的值。

UnivariateSpline 可以外插值，允许通过设置参数 ext= 0 或 “extrapolate” 外推插值。但如果外推范围过大，


3.2 Python 例程：一维插值（UnivariateSpline）

Python 例程：

# 3. 一维插值方法(外插)
import numpy as np
import matplotlib.pyplot as plt  # 导入 Matplotlib 工具包
from scipy.interpolate import UnivariateSpline  # 导入 scipy 中的一维插值工具 UnivariateSpline

# 生成已知数据点集 (x,y)，需插值的数据点集 xnew
x = np.linspace(0, 10, 11)  # 生成已知数据点集的 x
y = np.cos((x)**2/30)*2+2  # 生成已知数据点集的 y
xnew = np.linspace(-0.5, 10.5, 110)  # 指定需插值的数据点集 xnew

# 使用 UnivariateSpline 插值工具，由给定数据点集 (x,y) 求插值函数 fSpl
fSpl1 = UnivariateSpline(x, y, s=0)  # 三次样条插值，s=0：插值函数经过所有数据点
y1 = fSpl1(xnew)  # 由插值函数 fSpl1 计算插值点的函数值 y1

fSpl2 = UnivariateSpline(x, y)  # 三次样条插值，默认 s= len(w)
y2 = fSpl2(xnew)  # 由插值函数 fSpl2 计算插值点的函数值 y2

fSpl2.set_smoothing_factor(0.1)  # 设置光滑因子 sf
y3 = fSpl2(xnew)   # 由插值函数 fSpl2(sf=0.1) 计算插值点的函数值 y3

# 绘图
fig, ax = plt.subplots(figsize=(8,6))
plt.plot(x, y, 'ro', ms=5, label="data")
plt.plot(xnew, y1, 'm', label="3rd spline interpolate")
plt.plot(xnew, y2, 'g', label="3rd spline fitting")
plt.plot(xnew, y3, 'b--', label="smoothing factor")
ax.set_title("Data interpolate with extrapolation")
plt.legend(loc="best")
plt.show()

程序运行结果：

结果分析：

  1. 类 UnivariateSpline 既可以进行数据插值，也可以进行数据拟合。
    • 参数 s=0 时，要求样条函数通过所有数据点，即为数据插值；
    • s>0 或不设置参数 s 时，不要求样条函数通过所有数据点，就是用样条函数拟合给定的数据。
  2. 类 UnivariateSpline 进行数据插值或数据拟合，具有一定的外推能力，可以由插值函数进行外推计算。但插值算法的外推范围十分有限，容易发散或失真，使用时要非常谨慎。
  3. 图中绿色曲线是三次样条数据拟合的结果，蓝色曲线对拟合曲线作了进一步的平滑处理。

4. Scipy 二维插值方法

4.1 二维插值类 interp2d

Scipy.interpolate 中的 interp2d 类是一种基于固定数据点集创建函数的二维插值方法，经常用于二维图像处理。

类 interp2d 的定义和使用方法与一维插值 interp1d 类似。该类定义调用，允许使用 (x,y) 值调用对象，此时计算插值函数，并返回插值的 z轴值。具体地说，interp2d 类生成已知数据点集的插值函数 z=f(x,y)，通过调用这个插值函数，可以在已知数据之间插值，得到指定 (x,y) 的函数值 f(x,y)。

class scipy.interpolate.interp2d(x,y,z,kind=‘linear’,copy=True,bounds_error=False,fill_value=None))

主要参数：

  • x,y：一维数组，给定数据点集的 x,y 值。
  • z：一维数组，给定数据点集对应的函数值 z。
  • kind：字符串或整数，可选项，指定使用的样条曲线的种类或插值方法：‘linear’ 表示线性插值，‘cubic’ 表示三次插值，‘quintic’ 表示五次插值。默认值为 ‘linear’，即线性插值。

返回值：

  • 类 interp2d() 返回一个函数，其调用方法使用插值来查找新点的值。

注意：给定数据点集 x,y,z 有两种方式：

  1. 以一维数组 x 和 y 定义给定网络点集的列坐标和行坐标序列，一维数组 x 与 y 的长度 len(x)、len(y) 不一定相等，二维数组 z 的形状为 len(x)*len(y)， z i j = f ( x i , y i ) z_{ij}=f(x_i,y_i) zij​=f(xi​,yi​) 是数据网格点 (xi,yj) 的函数值。

  2. 以一维数组 x, y 定义给定网络点集的列坐标和行坐标序列，网络点集中每个点的坐标分别由二维数组 xx 和 yy 表示， 二维数组 z 是与 (xx,yy) 对应的网格中每一个数据点的函数值，$z_{ij} = f(xx_{ij}, yy_{ij}) $。


4.2 Python 例程：二维插值（interp2d）

使用示例：

# mathmodel24_v1.py
# Demo24 of mathematical modLSing algorithm
# Demo of interpolate with Scipy.interpolate
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-01

# 4. 二维插值方法
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp2d  # 导入 scipy 中的二维插值工具 interp2d

# 生成已知数据网格点集 (xx,yy,z)
x = np.linspace(-1, 1.5, 25)  # x 是一维数组
y = np.linspace(-1, 1, 20)  # y 是一维数组
xx, yy = np.meshgrid(x, y)  # 生成网格点的坐标 xx,yy (二维数组)
z = np.sin((xx+yy+xx**2+yy**2))  # 计算数据网格点的值 z=f(xx,yy)
print("shape of original dataset:\\n\\txx:{},yy:{},z:{}".format(xx.shape,yy.shape,z.shape))

# 由给定数据网格点集 (xx,yy,z) 求插值函数 fInterp: xx,yy,z 都是形状相同的二维数组
fInterp = interp2d(xx, yy, z, kind='cubic')  # 三阶样条插值

# 由插值函数 fInterp 计算需插值的数据点的函数值
xnew = np.linspace(-1, 1.5, 150)  # xnew 是一维数组
ynew = np.linspace(-1, 1, 100)  # ynew 是一维数组
znew = fInterp(xnew, ynew)  # 计算插值函数 fInterp 在 (xnew, ynew) 所描述网格点集的函数值
print("shape of interpolation dataset:\\n\\txnew:{},ynew:{},znew:{}".format(xnew.shape,ynew.shape,znew.shape))

# 绘图
fig = plt.figure(figsize=(10, 6))
ax1 = plt.subplot(1, 2, 1, projection='3d')
ax1.set_title("2-D original data")
# ax1.plot_wireframe(xx, yy, z, rstride=2, cstride=2, linewidth=1)
surf = ax1.plot_surface(xx, yy, z, rstride=2, cstride=2, cmap=plt.cm.coolwarm)
ax1.set_zlabel('zData')
xxnew, yynew = np.meshgrid(xnew, ynew)  # 将一维数组 xnew, ynew 转换为网格点集（二维数组）
print("\\txxnew:{},yynew:{},znew:{}".format(xxnew.shape,yynew.shape,znew.shape))
ax2 = plt.subplot(1, 2, 2, projection='3d')  # 3D 绘图要求 x,y,z 都是 n*m 二维数组
ax2.set_title("2-D interpolation data")
ax2.plot_wireframe(xxnew, yynew, znew, rstride=2, cstride=2,linewidth=1)
surf2 = ax2.plot_surface(xxnew, yynew, znew, rstride=2, cstride=2, cmap=plt.cm.coolwarm)
ax2.set_zlabel('zInterp')
plt.show()

运行结果：

shape of original dataset:
	xx:(20, 25),yy:(20, 25),z:(20, 25)
shape of interpolation dataset:
	xnew:(150,),ynew:(100,),znew:(100, 150)
	xxnew:(100, 150),yynew:(100, 150),znew:(100, 150)

注意：
由给定数据点集 (xx,yy,z) 求插值函数 fInterp 时 xx,yy,z 都是形状相同的二维数组（矩阵），由插值函数 fInterp 计算插值时所使用的 xnew,ynew 是一维数组，插值函数 fInterp 在 (xnew, ynew) 所构造的栅格点计算并返回插值函数值，返回值 z 是二维数组（矩阵）。


4.3 Python 例程：二维插值

通过设置 interp1d 类的参数kind，可以指定使用的样条曲线的种类或插值方法。

上节中介绍了 kind 各种选项所指的插值方法，本节进一步通过例程来比较不同方法的插值结果。

Python 例程：

# mathmodel24_v1.py
# Demo24 of mathematical modLSing algorithm
# Demo of interpolate with Scipy.interpolate
# Copyright 2021 YouCans, XUPT
# Crated：2021-08-01

# 5. 二维插值方法
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import interp2d  # 导入 scipy 中的二维插值工具 interp2d

# 生成已知数据网格点集 (xx,yy,z)
yy, xx = np.mgrid[-2:2:20j,-3:3:30j]  # 生成网格点 30x20 = 600
z = (1-0.5*xx+xx**5+yy**3) * np.exp(-xx**2-2*yy**2)  # 计算网格点的值 z
x, y = xx[0,:], yy[:,0]  # 由数据网格点 xx,yy 转换一维数组 x, y
print("shape of original dataset:\\n\\txx:{},yy:{},z:{}".format(xx.shape,yy.shape,z.shape))
print("\\tx:{},y:{},z:{}".format(x.shape, y.shape, z.shape)Python小白的数学建模课-12.非线性规划


 Python小白的数学建模课-20.网络流优化案例


 Python小白的数学建模课-03.线性规划


 Python小白的数学建模课-03.线性规划


 Python小白的数学建模课-03.线性规划


 Python小白的数学建模课-07 选址问题