插值与拟合

插值的作用

     数学建模比赛中，常常需要根据已知的函数点进行数据、模型的处理和分析，而有时候现有的数据是极少的，不足以支撑分析的进行，这时就需要使用一些数学的方法，模拟产生一些新的但又比较靠谱的值来满足需求，这就是插值的作用。

使用的库：from scipy import interpolate

interpolate.interp1d

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

     其中，x 是一维数据；y 是 N 维数据且 y 在插值的方向轴上长度必须和 x 相等；kind 指定了插值类型（linear，nearest，zero，slinear，quadratic，cubic，previous，next），其中 zero，slinear，quadratic 和 cubic 指定了零、一、二、三阶样条插值；previous 和 next 仅仅返回当前点的前一个/后一个值。默认为 linear 线性插值。
nearest：最邻近插值法              zero：阶梯插值
slinear、linear：线性插值         quadratic、cubic：2、3阶B样条曲线插值

from scipy.interpolate import interp1d

x = np.linspace(0, 10, num=11, endpoint=True)
y = np.cos(-x**2/9.0)
f1 = interp1d(x, y)
f2 = interp1d(x, y, kind='cubic')

xnew = np.linspace(0, 10, num=41, endpoint=True)
import matplotlib.pyplot as plt
plt.plot(x, y, 'o', xnew, f1(xnew), '-', xnew, f2(xnew), '--')
plt.legend(['data', 'linear', 'cubic'], loc='best')
plt.show()

结果展示：

interpolate.interp2d

Spline interpolation 样条插值

    样条插值需要两个基本步骤：
    计算曲线的样条表示 2. 在所需的点对样条曲线求值
为了找到样条曲线表示，有两种不同的方法来表示曲线并获得（平滑）样条系数：直接法 和 参数法。

    直接法： 用函数 splrep 求二维平面上曲线的样条表示。前两个是必要参数，提供了曲线上点的 x、y 坐标，函数返回值为三元组 ( t , c , k )， t 为结点，c 为系数，k 为样条顺序（默认为三次），可以通过输入参数 k 进行更改。

    参数法：对于 N 维空间的曲线，函数 splprep 可以参数化地定义曲线。该函数第一个为必要参数，该输入为 N 维数组列表，表示了在 N 维空间中的曲线。数组长度是曲线上点的数目。函数返回值为 三元组 ( t , c , k ) 和变量 u。
    举例：

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

# Cubic-spline
x = np.arange(0, 2*np.pi+np.pi/4, 2*np.pi/8)
y = np.sin(x)
tck = interpolate.splrep(x, y, s=0)
xnew = np.arange(0, 2*np.pi, np.pi/50)
ynew = interpolate.splev(xnew, tck, der=0)

plt.figure()
plt.plot(x, y, 'x', xnew, ynew, xnew, np.sin(xnew), x, y, 'b')
plt.legend(['Linear', 'Cubic Spline', 'True'])
plt.axis([-0.05, 6.33, -1.05, 1.05])
plt.title('Cubic-spline interpolation')
plt.show()

    x 和 y 是给定的一些坐标点，通过 splrep 函数计算出 B-spline 曲线的参数，再将该参数传递给 splev 函数计算出各个取样点 xnew 的插值结果。

• (x,y)：实际给定的坐标点，表示为图中的 x 符号
• ( x n e w , y n e w )：xnew 是重新取的一系列点的横坐标，ynew 是根据前面得出的插值参数 tck，得到的一系列点对应的纵坐标值，它们共同构成了插值后的曲线点.
• ( x n e w , s i n ( x n e w ) ) ：实际 xnew 应该对应的 y 值，相当于 ground-truth（表示为图中绿色曲线）

1. interpolate.splrep

    作用：求取一维曲线的 B-spline 插值
给定一组数据点 ( x [ i ] ， y [ i ] )，确定在区间 x b ≤ x ≤ x e 上 k 次的光滑样条逼近。
语法：

scipy.interpolate.splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=None, t=None, full_output=0, per=0, quiet=1)

    其中，x 和 y 是点的横纵坐标，xb 和 xe 是间隔，k 为样条拟合度，建议使用三次样条（cubic splines），通常 1 ≤ k ≤ 5

# splrep 应用
import matplotlib.pyplot as plt
from scipy.interpolate import splev, splrep
x = np.linspace(0, 10, 10)
y = np.sin(x)
spl = splrep(x, y)
x2 = np.linspace(0, 10, 200)
y2 = splev(x2, spl)
plt.plot(x, y, 'o', x2, y2)
plt.show()

2. interpolate.splprep
语法：

scipy.interpolate.splprep(x, w=None, u=None, ub=None, ue=None, k=3, task=0, s=None, t=None, full_output=0, nest=None, per=0, quiet=1)

作用： 求取 N 维曲线的 B-spline 插值

phi = np.linspace(0, 2.*np.pi, 40)
r = 0.5 + np.cos(phi)         # polar coords
x, y = r * np.cos(phi), r * np.sin(phi)    # convert to cartesian

from scipy.interpolate import splprep, splev
tck, u = splprep([x, y], s=0)
new_points = splev(u, tck)

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(x, y, 'ro')
ax.plot(new_points[0], new_points[1], 'r-')
plt.show()

x：列表类型 [x,y]，x 和 y 分别是表示曲线上坐标点的横纵坐标。
w：数组类型，长度与 x[0] 相同，默认为 ones(len(x[0]))。
k：spline 等级，推荐使用 Cubic splines（k=3），默认为 3。
s：float 类型，平滑条件，较大的 s 意味着曲线比较光滑，建议 s 的取值在 ( m − s q r t ( 2 ∗ m ) , m + s q r t ( 2 ∗ m ) ) 范围内。

     下图中，( x , y ) 是实际定义的坐标点，表示为红色圆点；通过 splprep 函数拟合 x，y，得到拟合参数 tck 和 u，得到新的拟合曲线 new_points。

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

x = np.linspace(0, 2*np.pi+np.pi/4, 10)
y = np.sin(x)

x_new = np.linspace(0, 2*np.pi+np.pi/4, 100)
f_linear = interpolate.interp1d(x, y)
tck = interpolate.splrep(x, y)
y_bspline = interpolate.splev(x_new, tck)

plt.xlabel(u'安培/A')
plt.ylabel(u'伏特/V')

plt.plot(x, y, "o",  label=u"原始数据")
plt.plot(x_new, f_linear(x_new), label=u"线性插值")
plt.plot(x_new, y_bspline, label=u"B-spline插值")

pl.legend()
pl.show()

实例分析：

import numpy as np
from scipy import interpolate
import pylab as pl
#创建数据点集并绘制
pl.figure(figsize=(12,9))
x = np.linspace(0, 10, 11)
y = np.sin(x)
ax=pl.plot()

pl.plot(x,y,'ro')
#建立插值数据点
xnew = np.linspace(0, 10, 101)
for kind in ['nearest', 'zero','linear','quadratic']:
    #根据kind创建插值对象interp1d
    f = interpolate.interp1d(x, y, kind = kind)
    ynew = f(xnew)#计算插值结果
    pl.plot(xnew, ynew, label = str(kind))

pl.xticks(fontsize=20)
pl.yticks(fontsize=20)

pl.legend(loc = 'lower right')
pl.show()

参考博客一
参考博客二

---

最小二乘法

最小二乘法拟合直线的原理：

使用的库 : from scipy.optimize import leastsq
函数调用形式 :

scipy.optimize.leastsq(func, x0, args=(), Dfun=None, full_output=0, col_deriv=0, ftol=1.49012e-08, xtol=1.49012e-08, gtol=0.0, maxfev=0, epsfcn=None, factor=100, diag=None)

参数说明 : （一般指定前三个参数即可）

• func：callable(可调用的)。
           应该至少接受一个(可能为N个向量)长度的参数，并返回M个浮点数。它不能返回NaN，否则拟合可能会失败。
         func 是我们自己定义的一个计算误差的函数。
• x0：ndarray
           最小化的起始估算。      x0 是计算的初始参数值。
• args：tuple, 可选参数。 函数的所有其他参数都放在此元组中。
        args 是指定func的其他参数。

参考博客

案例分析：

import numpy as np
from scipy import linalg
from scipy.optimize import leastsq

# 需要拟合的函数func
def func(x, p):
    a, b, c = p
    return a * np.power(x, 2) + b * x + c

    下面代码是定义计算误差函数，计算实验数据x、y和拟合函数之间的差，参数p为拟合需要找到的系数，同时添加噪声数据。

# 误差函数
def error(p, y, x):
    return y - func(x, p)

# 读入训练数据
data = pd.read_excel('../datas/data3.xlsx')
X = data['压力(MPa)'].values
Y = data['弹性模量(MPa)'].values
p0 = [7, 0.2, 0]  # 第一次猜测的函数拟合参数,可随机给出

    调用leastsq拟合，error为计算误差函数，p0为拟合参数初始值，args为需要拟合的实验数据。

Para = leastsq(error, p0, args=(X, Y))
print("拟合参数:", Para[0])

画图·：

import matplotlib.pyplot as plt
import pylab as pl

plt.rcParams['font.sans-serif'] = ['SimHei']          # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False            # 用来正常显示负号
pl.plot(X, Y, marker='+', label=u"真实数据")
pl.plot(X, func(X, Para[0]), label=u"拟合数据")
pl.legend()
pl.show()

---

多项式拟合

from skimage.metrics import mean_squared_error  # 求解均方误差
# 显示中文（因为横纵坐标想要显示中文的）
plt.rcParams['font.sans-serif'] = ['SimHei']
# 读入训练数据
data = pd.read_excel('../datas/附件1-凸轮边缘曲线.xlsx')

train_x = data['极角'].values
train_y = data['极径'].values

# 设置 x 轴的刻度大小
plt.xticks(np.arange(8))
# 设置 X 轴的网格线，风格为 点画线
plt.grid(axis='x', linestyle='-')

# 设置 y 轴的刻度大小
plt.yticks(np.arange(9))
# 设置 y 轴的网格线，风格为 点画线
plt.grid(axis='y', linestyle='-')
# 设置拟合的多项式的阶数
f = np.polyfit(train_x, train_y, 6)  # 返回多项式各项系数
p = np.poly1d(f)  # 得到拟合的多项式表达式
# 拟合y值
yvals = p(train_x)
# 拟合度
print(1 - mean_squared_error(train_y, yvals))
print('=' * 50)
# 绘图
plt.plot(train_x, train_y, 'o', label='原始数据')
plt.plot(train_x, yvals, 'r', label='拟合数据')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4)  # 指定legend的位置右下角
plt.show()