拟合与插值|线性最小二乘拟合|非线性最小二乘拟合|一维插值|二维插值

挖掘数据背后的规律是数学建模的重要任务，拟合与插值是常用的分析方法

• 掌握拟合与插值的基本概念和方法
• 熟悉Matlab相关程序实现
• 能够从数据中挖掘数学规律

拟合问题的基本提法

拟合问题的概念

已知一组数据(以二维为例)，即平面上n个点$(x_i,y_i),i=1,2,\dots ,n$，寻求一个函数$y=f(x)$，使得$f(x)$在某种准则下与所有数据点最为接近

$$\sum _{i=1}^n({\delta }_i{)}^2\to min$$

1. 确定$f(x)$表达式的形式
2. 估计$f(x)$中待定系数

拟合要解决的两个基本问题

确定$f(x)$表达式的形式

1. 数据可视化，直观判断$f(x)$的形式
   
   

2. 通过机理分析来确定$f(x)$的形式
   饮酒驾车中酒精浓度拟合问题
   $$y=k{e}^{-at}$$
   SARS的传播规律拟合问题
   $$i(t)=\frac{1}{1+\left(\frac{1}{i_0}-1\right)e^{-\lambda t}}$$
   估计$f(x)$中的待定参数
   根据$f(x)$表达式中，待定参数具体形式的不同，一般有两种方法
   线性最小二乘法：适用于$f(x)$中待定参数为线性形式
   非线性最小二乘法：适用于$f(x)$中待定参数为非线性形式

线性最小二乘拟合及其实现

线性最小二乘拟合

基本思路

1. 选定一组基本函数$r_1(x),r_2(x),\dots ,r_m(x);m<n$，确定$f(x)$表达式
   $$f(x)=a_1{r}_1(x)+a_2{r}_2(x)+\cdots +a_m{r}_m(x)$$
   其中，$a_1,a_2,\dots ,a_m$为待定参数
2. 按照最小二乘原则确定待定系数
   求$a_1,a_2,\dots ,a_m$使n个点$(x_i,y_i)$与曲线$f(x)$的距离平方和$J(a_1,a_2,\dots ,a_m)$最小
   $$mi{n}_{a_1,a_2,\dots ,a_m}J(a_1,a_2,\dots ,a_m)=\sum _{i=1}^n({\delta }_i{)}^2$$
   $$=\sum _{i=1}^n(y_i-a_1{r}_1{x}_i-a_1{r}_1{x}_i-\cdots -a_m{r}_m{x}_i{)}^2$$

优化问题的求解
超定方程组，方程个数大于未知量个数的方程组
$$\{\begin{array}{c}{r}_{11}{a}_1+r_{12}{a}_2+\cdots +r_{1m}{a}_m=y_1\\ r_{21}{a}_1+r_{22}{a}_2+\cdots +r_{2m}{a}_m=y_2\\ \dots \dots \\ r_{n1}{a}_1+r_{n2}{a}_2+\cdots +r_{nm}{a}_m=y_n\end{array}$$
$$R=\left(\begin{array}{ccccccc}{r}_{11}& & r_{12}& & \dots & & r_{1m}\\ r_{21}& & r_{22}& & \dots & & r_{2m}\\ \dots & & \dots & & \dots & & \dots \\ r_{n1}& & r_{n2}& & \dots & & r_{nm}\end{array}\right)$$
$$a=\left(\begin{array}{c}{a}_1\\ a_2\\ \dots \\ a_m\end{array}\right)y=\left(\begin{array}{c}{y}_1\\ y_2\\ \dots \\ y_n\end{array}\right)$$
超定方程组可改写为：$Ra=y$
一般地，超定方程组是不存在解的矛盾方程组，不存在准确解
需要在最小二乘的意义下求解
$$\sum _{i=1}^n({\delta }_i{)}^2=\sum _{i=1}^n(y_i-a_1{r}_1{x}_i-a_1{r}_1{x}_i-\cdots -a_m{r}_m{x}_i{)}^2$$
最小化的解$a_1,a_2,\dots ,a_m$，称为超定方程组的线性最小二乘解

线性最小二乘的求解

当$R^{T}R$可逆时，超定方程组$Ra=y$存在最小二乘解，且为如下正规方程组的解
$$R^{T}Ra=R^{T}y$$
$$a=(R^{T}R{)}^{-1}{R}^{T}y$$
线性最小二乘的Matlab实现
对超定方程组$R_{n\times m}{a}_{m{\times }_1}=y_{n{\times }_1}$

a = R \ y

即可得到相应的最小二乘解

例

求二次多项式
$$f(x)=a_1{x}^2+a_{2}x+a_3$$
中的系数$a_1,a_2,a_3$
$$\left(\begin{array}{ccccc}{x}_1^2& & x_1& & 1\\ x_2^2& & x_2& & 1\\ \dots & & \dots & & \dots \\ x_{11}^2& & x_{11}& & 1\end{array}\right)=\left(\begin{array}{c}{a}_1\\ a_2\\ a_3\end{array}\right)=\left(\begin{array}{c}{y}_1\\ y_2\\ \dots \\ y_{11}\end{array}\right)$$

x = 0 : 0.1 : 1;
y = [-0.447, 1.978, 3.28, 6.16, 7.08, 7.34, 7.66, 9.58, 9.48, 9.30, 11.2];
R = [(x.^2)',x',ones(11,1)];
A=R\y'

计算结果

A = [-9.8108, 20.1293, -0.0317]

$$f(x)=-9.8108x^2+20.1293x-0.0317$$

非线性最小二乘拟合及其实现

基本思路

1. 确定问题的函数表达式的形式
   $$y=f(a_1,a_2,\dots ,a_m;x)$$
   其中$a_1,a_2,\dots ,a_m$为待定参数
2. 按照最小二乘原则确定待定参数
   在最小二乘意义下，求使得
   $$mi{n}_{a_1,a_2,\dots ,a_m}J(a_1,a_2,\dots ,a_m)=\sum _{i=1}^n({\delta }_i{)}^2$$
   $$=\sum _{i=1}^n(y_i-f(a_1,a_2,\dots ,a_m;x_i){)}^2$$
   最小化的解$\widehat{a_1},\widehat{a_2},\dots ,\widehat{a_m}$，称为待定系数$a_1,a_2,\dots ,a_m$的非线性最小二乘解

非线性最小二乘拟合的Matlab实现

两个非线性最小二乘拟合函数

• lsqcurvefit
• lsqnonlin

非线性最小二乘求解命令lsqcurvefit
已知数据点
$$xdata=(xdat{a}_1,xdat{a}_2,\dots ,xdat{a}_n{)}^T$$
$$ydata=(ydat{a}_1,ydat{a}_2,\dots ,ydat{a}_n{)}^T$$
lsqcurvefit通过求解如下最小二乘问题，得到参数向量x的估计
$$min||ydata-F(x,xdata)|{|}_2^2=mi{n}_x\sum _i(ydat{a}_i-F(x,xdat{a}_i{)}^2$$
其中，$F(x,xdata)$表示预测值向量
$$F(x,xdata)=\left(\begin{array}{c}F(x,xdat{a}_1)\\ F(x,xdat{a}_2)\\ \dots \\ F(x,xdat{a}_n)\end{array}\right)$$
基本格式

x = lsqcurvefit('fun', x0, xdata, ydata);

• fun是预先定义的函数，实现待拟合非线性函数$F(x,xdata)$
• x0给定参数初始值(lsqcurvefit求解非线性优化问题的初始迭代点)
• xdata和ydata是样本数据的输入和输出向量(或矩阵)
• x为未知参数的非线性最小二乘估计

例

$$F(x,xdata)=\left(\begin{array}{c}F(x,xdat{a}_1)\\ F(x,xdat{a}_2)\\ \dots \\ F(x,xdat{a}_n)\end{array}\right)$$

• x，三个待定参数组成的向量
• xdata，所有的自变量的取值组成的向量
• 输出，是待拟合的函数在每一个样本点处对应的拟合值
  $$=\left(\begin{array}{c}a+b{e}^{-0.02k{t}_1}\\ a+b{e}^{-0.02k{t}_2}\\ \dots \\ a+b{e}^{-0.02k{t}_{11}}\end{array}\right)$$

function f = curvefun(x, tdata)
f = x(1) + x(2)*exp(-0.02*x(3)*tdata) %其中x(1)=a;x(2)=b;x(3)=k

%主程序命令
tdata = 100:100:1000
ydata = 1e-03*[4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59];
x0 = [0.1; 0.1; 0.1];
x = lsqcurvefit('curvefun', x0, xdata, ydata)
y_pred = curvefun(x, tdata)
plot(tdata, ydata, 'bo', tdata, y_pred, 'r.:') %可视化预测结果

运算结果

x = 0.0070  -0.0030  0.1012
y_pred = 0.0045  0.0050  0.0054  0.0057  0.0059  0.0061  0.0063  0.0064  0.0065  0.0066

结论
a = 0.0070, b=-0.0030, k=1012
$$f(t)=0.0070-0.0030e^{-0.02\times 0.1012\times t}$$

一维插值及其实现

一维插值的定义

已知n+1个节点$(x_i,y_i)；i=0,1,\dots ,n$，其中$x_i$互不相同，设
$$a=x_0<x_1<\cdots <x_n=b$$
构建一个函数$y=f(x)$通过全部节点，即
$$f(x_i)=y_i;i=0,1,\dots ,n$$
用函数$f(x)$计算任一被插值点$x^{\ast }$处的插值，即
$$y^{\ast }=f(x^{\ast })$$

拉格朗日插值

已知n+1个节点$(x_i,y_i)；i=0,1,\dots ,n$，构造一个n次多项式$P_n(x)$，满足
$$P_n(x_i)=y_i;i=1,2,\dots ,n$$
被称为拉格朗日插值多项式，构造方法如下
$$P_n(x)=\sum _{i=1}^n{L}_i(x)y_i$$
其中，$L_i(x)$为n次多项式
$$L_i(x)=\frac{(x-x_0)(x-x_1)\dots (x-x_{i-1})(x-x_{i+1})\dots (x-x_n)}{(x_i-x_0)(x_i-x_1)\dots (x_i-x_{i-1})(x_i-x_{i+1})\dots (x_i-x_n)}$$
称为拉格朗日插值基函数
$$L_i(x_j)=\{\begin{array}{c}0,i\ne j\\ 1,i=j\end{array}$$
$$P_n(x_j)=y_{j}j=0,1,\dots ,n$$
两点一次线性插值多项式
$$P_1(x)=\frac{(x-x_1)}{(x_0-x_1)}{y}_0+\frac{(x-x_0)}{(x_1-x_0)}{y}_1$$
三点二次抛物插值多项式
$$P_2(x)=\frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)}{y}_0+\frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)}{y}_1+\frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)}{y}_2$$
有可能产生震荡现象

例

随着节点个数的增加，多项式阶数不断增大，$P_n(x)$产生较大的震荡现象

分段线性插值

基本思想：构造一条一条的直线段进行插值计算
$$L_n(x)=\sum _{i=1}^n{l}_i(x)y_i$$
其中
$$l_i(x)=\{\begin{array}{c}\frac{x-x_{i-1}}{(x_i-x_{i-1})},x_{i-1}\le x<x_i\\ \frac{x-x_{i+1}}{(x_i-x_{i+1})},x_i\le x<x_{i+1}\\ 0,otherwise\end{array}$$

• 思路简单易实现
• 计算量小
• 不会产生震荡现象
• 光滑性不好
  

三次样条插值

基本思想
构造一段一段的三次多项式曲线进行插值计算
$$s_i(x)=a_i(x{)}^3+b_i(x{)}^2+c_{i}x+d_i,i=1,2,\dots ,n$$

三次样条插值函数
S(x)={si(x),x∈[xi−1,xi]∣i=1,2,…,n} S(x)=\left \{ s_{i}(x),x\in [x_{i-1},x_{i}]| i=1,2,\dots,n \right \} S(x)={si​(x),x∈[xi−1​,xi​]∣i=1,2,…,n}
利用如下三类条件

1. 插值条件n+1
   $$S(x_i)=y_i,i=0,1,2,\dots ,n$$
2. 二阶光滑3n-3
   $$S(x)\in C^2[x_0,x_n]$$
3. 自然边界条件2
   $$S^{\prime \prime }(x_0)=S^{\prime \prime }(x_n)=0$$
   可求得三次样条插值函数的所有待定系数

一维插值的Matlab实现

yi = interp1(x, y, xi, 'method')

• x，y，插值节点对应的坐标
• xi，被插值节点的坐标
• yi，与xi对应的插值结果
• ‘method’，插值方法
  nearest：最邻近插值
  linear：分段线性插值
  spline：三次样条插值
  cubic：立方插值
  缺省时：分段线性插值
• x是单调的
• xi不能够超过x的范围(内插)

例

hours = 1:12;
temps =[5 8 9 15 25 29 31 30 22 25 27 24];
h= 1 : 0.1 : 12；
t_nearest = interp1(hours, temps, h, 'nearest');%邻近插值
t_linear = interp1(hours, temps, h, 'linear');%分段线性插值
t_spline = interp1(hours, temps, h, 'spline');%三次样条插值

• hours，原始的插值节点，横坐标的取值向量
• temps，因变量温度的取值
• h，被插值节点对应的向量，每隔1/10小时产生一个时间点对应的等差数列

subplot(2,2,1), plot(hours, temps, 'bo')
subplot(2,2,2), plot(hours, temps, 'bo', h, t_nearest,'r:')%邻近插值
subplot(2,2,3), plot(hours, temps, 'bo', h, t_linear,  'r:')%线性插值
subplot(2,2,4), plot(hours, temps, 'bo', h, t_spline, 'r:')%三次样条插值

二维插值及其实现

二维插值定义

1. 网格节点
   已知mxn个节点
   $$(x_i,y_i,z_{ij});i=1,\dots ,m;j=1,\dots ,n$$
   其中$x_i,y_i$互不相同，设
   $$a=x_1<x_2<\cdots <x_m=b$$
   $$c=y_1<y_2<\cdots <y_n=d$$
   构造一个二元函数$z=f(x,y)$，通过已知节点，即
   $$f(x_i,y_j)=z_{ij};i=1,\dots ,m;j=1,\dots ,n$$
   并由此计算插值
   $$z^{\ast }=f(x^{\ast },y^{\ast })$$
   

2. 散乱节点
   已知n个节点
   $$(x_i,y_i,z_i);i=1,\dots ,n$$
   其中$x_i,y_i$互不相同
   构造一个二元函数$z=f(x,y)$，通过已知节点，即
   $$f(x_i,y_j)=z_i;i=1,\dots ,n$$
   并由此计算插值
   $$z^{\ast }=f(x^{\ast },y^{\ast })$$
   

二维插值方法

1. 最邻近插值
   基本思想
   将与被插值点最邻近节点的函数值作为插值结果
   

• 思路直观，易理解
• 计算量小，易于实现
• 不连续

2. 分片线性插值
   基本思想
   构造一片一片的空间三角平面进行插值计算
   

• 直观，易于实现
• 连续
• 光滑性不好

3. 双线性插值
   基本思想
   构造一片一片的双线性空间二次曲面进行插值计算
   双线性插值函数表达式为
   $$f(x,y)=(ax+b)(cy+d)$$
   四个待定系数的确定：
   利用矩形四个顶点(插值节点)处的函数值，建立四个方程进行求解
   

• 易理解实现
• 函数连续
• 光滑性不好

4. 三次卷积插值，三次样条插值
   保证光滑

二维插值的Matlab实现

1. 网格节点

z = interp2(x0, y0, z0, x, y, 'method')

• x0，y0，z0，指定插值节点对应的坐标
• x，y，指定被插值节点的信息
• z，被插值节点的插值结果
• method，插值方法
  nearest：最邻近插值
  linear：线性插值$C^0$
  spline：三次样条插值$C^2$
  cubic：三次卷积插值$C^1$
  缺省时：线性插值
• x0，y0单调
• x，y为同型矩阵，或一个行向量另一个列向量
• x，y不能够超过x0，y0的范围(内插)

2. 散乱节点

z = gridata(x0, y0, z0, x, y, 'method')

• x0，y0，z0，指定插值节点对应的坐标
• x，y，指定被插值节点的信息
• z，被插值节点的插值结果
• method，插值方法
  nearest：最邻近插值
  linear：线性插值$C^0$
  cubic：三次插值$C^2$
  v4：双调和样条插值$C^2$
  缺省时：线性插值
• x，y为同型矩阵，或一个行向量另一个列向量
• x，y不能够超过x0，y0的范围(内插)

例1

画出原始数据的温度分布图

x=1:5;
y=1:3;
temps=[82 81 80 82 84;
	79 63 61 65 81;
	84 84 82 85 86];
mesh(x,y, temps)

然后，在x，y方向上每隔0.2个单位处进行插值，
再输入以下命令：

xi = 1: 0.2 : 5;
yi = 1: 0.2 : 3;
zi = interp2(x, y, temps, xi', yi, 'cubic');
mesh(xi, yi, zi)

画出插值后的温度分布曲面图

例2

画出原始海拔高度数据分布图

x=1200: 400: 4000;
y=[1200 1600 2400 2800 3200 3600];
[X,Y] = meshgrid(x,y);
Z=[1130 1250 1280 1230 1040 900 500 700;
1320 1450 1420 1400 1300 700 900 850;
1390 1500 1500 1400 900 1100 1060 950;
1500 1200 1100 1350 1450 1200 1150 1010;
1500 1550 1600 1550 1600 1600 1600 1550;
1480 1500 1550 1510 1430 1300 1200 980];
surfc(X, Y, Z);

• meshgrid，对分割向量进行网格交叉，得到网格矩阵
• surfc，三维空间曲面绘制
  

然后，对x，y方向上的取值向量加密，定义被插值数据点

xi = linspace(1200, 4000, 50);
yi = linspace(1200, 3600, 50);
[XI,YI]=meshgrid(xi,yi);

• linspace，从1200到4000，产生50个均匀分布的点，构成等差数列
  进行插值计算，并可视化

methods = {'nearest', 'linear', 'spline', 'cubic'};
for i = 1 : 4
	ZI = interp2(X, Y, Z, XI, YI, methods{i});
	subplot(2, 2, i);
	surfc(XI, YI, ZI);%这个函数直接就是把等高线画在下面
	title(['方法:',methods{i}]);
	xlabel('x');
	ylabel('y');
end