MATLAB回归、插值、逼近、拟合【转载】
原文链接：
https://blog.csdn.net/daaikuaichuan/article/details/73870209
https://www.cnblogs.com/vanker/p/4268292.html
一、回归、插值、逼近、拟合的区别
1、总体概况
插值和拟合都是函数逼近或者数值逼近的重要组成部分。他们的共同点都是通过已知一些离散点集M上的约束，求取一个定义在连续集合S(M包含于S)的未知连续函数，从而达到获取整体规律的目的，即通过"窥几斑"来达到"知全豹"。
（1）简单的讲
- 所谓拟合是指已知某函数的若干离散函数值{f1,f2,…,fn}，通 过调整该函数中若干待定系数f(λ1, λ2,…,λ3), 使得该函数与已知点集的差别(最小二乘意义)最小。如果待定函数是线性，就叫线性拟合或者线性回归(主要在统计中)，否则叫作非线性拟合或者非线性回归。其表达式也可以是分段函数，这种情况下叫作样条拟合；
- 而插值是指已知某函数的在若干离散点上的函数值或者导数信息，通过求解该函数中待定形式的插值函数以及待定系数，使得该函数在给定离散点上满足约束。插值函数又叫作基函数，如果该基函数定义在整个定义域上，叫作全域基，否则叫作分域基。如果约束条件中只有函数值的约束，叫作Lagrange插值，否则叫作Hermite插值。
（2）从几何意义上讲
- 拟合是给定了空间中的一些点，找到一个已知形式未知参数的连续曲面来最大限度地逼近这些点
- 而插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。
2、分点探讨
- 回归一般指线性回归，是求最小二乘解的过程。在求回归前，已经假设所有型值点同时满足某一曲线方程，计算只要求出该方程的系数
- 多项式插值：用一个多项式来近似代替数据列表函数，并要求多项式通过列表函数中给定的数据点。（插值曲线要经过型值点。）
- 多项式逼近：为复杂函数寻找近似替代多项式函数，其误差在某种度量意义下最小。（逼近只要求曲线接近型值点，符合型值点趋势。）
- 多项式拟合：在插值问题中考虑给定数据点的误差，只要求在用多项式近似代替列表函数时，其误差在某种度量意义下最小。
3、注意
- 表列函数：给定n+1个不同的数据点（x0,y0）,(x1,y1)...,(xn,yn)，称由这组数据表示的函数为表列函数。
- 逼近函数：求一函数，使得按某一标准，这一函数y=f（x）能最好地反映这一组数据即逼近这一表列函数，这一函数y=f（x）称为逼近函数
- 插值函数：根据不同的标准，可以给出各种各样的函数，如使要求的函数y=f(x) 在以上的n+1个数据点出的函数值与相应数据点的纵坐标相等，即yi=f(x1) (i=0, 1, 2, ..., n) 这种函数逼近问题称为插值问题，称函数y=f(x) 为数据点的插值函数，xi称为插值点。
二、经典插值方法
1、拉格朗日插值
1 function y=lagrange(x0,y0,x) 2 %拉格朗日插值函数 3 %n 个节点数据以数组 x0, y0 输入(注意 Matlat 的数组下标从1开始), 4 %m 个插值点以数组 x 输入,输出数组 y 为 m 个插值 5 n=length(x0);m=length(x); 6 for i=1:m 7 z=x(i); 8 s=0.0; 9 for k=1:n 10 p=1.0; 11 for j=1:n 12 if j~=k 13 p=p*(z-x0(j))/(x0(k)-x0(j)); 14 end 15 end 16 s=p*y0(k)+s; 17 end 18 y(i)=s; 19 end
应用实例：
1 x0=1:1:20; 2 y0=x0.^2-20*x0-5; 3 x=1:0.1:20; 4 z=lagrange(x0,y0,x); 5 plot(x,z,':',x0,y0,'ko');
2、分段线性插值
MATLAB现成的插值函数为interp1，其调用格式为： yi= interp1(x,y,xi,'method')
其中x，y为插值点，yi为在被插值点xi处的插值结果；x,y为向量， 'method'表示采用的插值方法，包括：
'method'：是最近项插值;
'linear'：线性插值;（默认）
'spline'：逐段3次样条插值; (下面的三次样条插值会用到)
'cubic'：保凹凸性3次插值
'pchip'：分段三次Hermite 插值。
例如：
在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别为
12，9，9，1，0，18 ，24，28，27，25，20，18，15，13，
推测中午12点（即13点）时的温度．
1 x=0:2:24; 2 y=[12,9,9,10,18,24,28,27,25,20,18,15,13]; 3 x1=0:0.5:24; 4 y1=interp1(x,y,x1,'linear'); 5 plot(x,y,'bo',x1,y1,'r:');
3、埃米尔特插值
如果要求插值函数不仅在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至高阶导数值，这就是埃尔米特插值问题。
已知f（x）的n+1个节点的函数值f（xi）以及导数值f`（xi），可得一个至多n+1次的多项式H（x），即hermite插值多项式。新建以下这个函数：
1 function y = hermite( x0,y0,y1,x ) 2 %埃尔米特插值多项式 3 %x0为点横坐标 4 %y0为函数值 5 %y1为导数值 6 %m个插值点用数组x输入 7 n=length(x0);m=length(x); 8 for k=1:m 9 yy=0.0; 10 for i=1:n 11 h=1.0; 12 a=0.0; 13 for j=1:n 14 if j~=i 15 h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2; 16 a=1/(x0(i)-x0(j))+a; 17 end 18 end 19 yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i)); 20 end 21 y(k)=yy; 22 end
4、样条插值
所谓样条（ Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并使连接点处有连续的曲率。数学上将具有一定光滑性的分段多项式称为样条函数。在实际中最常用的是二次样条函数和三次样条函数：
（1）二次样条函数插值
首先，我们注意到s2 (x)中含有 n + 2 个特定常数，故应需要 n + 2 个插值条件，因此，二次样条插值问题可分为两类：
- 已知插值节点xi 和相应的函数值 yi (i = 0,1,…,n) 以及端点 x0 （或 xn ）处的导数值y'0（或y'n）
- 已知插值节点xi 和相应的导数值 y'i (i = 0,1,…,n) 以及端点 x0 （或 xn ）处的函数值y0 （或yn ）
（2）三次样条函数插值
由于 s3 (x)中含有n + 3 个待定系数，故应需要 n + 3 个插值条件，已知插值节点xi 和相应的函数值 f(xi ) = yi (i = 0,1,…,n) ，这里提供了 n + 1 个条件，还需要 2 个边界条件。因此，三次样条插值问题可分为三类：
- s'3 (a) = y'0 ,s'3 (b) = y'n 。由这种边界条件建立的样条插值函数称为 f(x) 的完备三次样条插值函数。特别地,y0' = yn`= 0时，样条曲线在端点处呈水平状态。如果 f' (x) 不知道，我们可以要求 s'3 (x) 与 f' (x) 在端点处近似相等。这时以x0 , x1 , x2 , x3 为节点作一个三次 Newton 插值多项式 Na (x) ，以 xn, xn−1, xn−2, xn−3 作一个三次 Newton 插值多项式 Nb (x) ，要求s' (a) = N'a (a), s' (b) = N'b (b)由这种边界条件建立的三次样条称为 f(x) 的 Lagrange 三次样条插值函数。
- s"3 (a) = y"0 ,s"3 (b) = y"3 。特别地 y"n = y"n = 0 时，称为自然边界条件。
- s'3 ( a + 0) = s'3 ( b − 0), s"3 (a + 0) = s"3 (b − 0) ， (这里要求 s3 (a + 0) =s3 (b − 0) )此条件称为周期条件。
（3）Matlab实现（三次样条插值）
Matlab中的函数：
1、y=interp1(x0,y0,x,`spline`);%(spline改成linear，则变成线性插值)
2、y=spline(x0,y0,xi);%这个是根据己知的x，y数据，用样条函数插值出xi处的值。即由x,y的值计算出xi对应的函数值。
3、pp=spline(x0,y0);%是由根据己知的x，y数据，求出它的样条函数表达式，不过该表达式不是用矩阵直接表示，要求点x`的值，要用函数y`=ppval(pp,x`);
4、pp=csape(x,y,'变界类型','边界值conds');生成各种边界条件的三次样条插值. 其中,(x,y)为数据向量,边界类型可为:
'complete'：给定边界一阶导数,即默认的边界条件,Lagrange边界条件
'not-a-knot'：非扭结条件,不用给边界值.
'periodic'：周期性边界条件,不用给边界值.
'second'：给定边界二阶导数.
'variational'：自然样条(边界二阶导数为[0,0]
边界值conds可用1x2矩阵表示,矩阵元素取值为1,2,此时，使用命令pp=csape(x0,y0_ext,conds)其中 y0_ext=[left, y0, right]，这里 left 表示左边界的取值， right 表示右边界的取值。conds(i)=j 的含义是给定端点 i 的 j 阶导数， 即 conds 的第一个元素表示左边界的条件，第二个元素表示右边界的条件， conds=[2,1]表示左边界是二阶导数，右边界是一阶导数，对应的值由 left 和 right 给出。
例子：
表 1
x 0 3 5 7 9 11 12 13 14 15
y 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6
要求用 Lagrange、分段线性和三次样条三种插值方法计算。
编程实现：
1 clear,clc 2 x0=[0,3,5,7,9,11,12,13,14,15]; 3 y0=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6]; 4 t=0:0.05:15; 5 %拉格朗日插值函数 6 y1=lagrange(x0,y0,t);%调用编写的lagrange函数 7 dy1=(lagrange(x0,y0,0.0001)-lagrange(x0,y0,0))/0.0001%x=0处斜率 8 min1=min(lagrange(x0,y0,13:0.001:15))%13到15最小值 9 subplot(2,2,1); 10 plot(x0,y0,'ro',t,y1);%画出曲线 11 title('拉格朗日插值函数'); 12 %分段线性插值 13 y2=interp1(x0,y0,t,'spline');%注意区分spline与linear 14 Y2=interp1(x0,y0,t);%默认linear 15 dy2=(interp1(x0,y0,0.0001,'spline')-interp1(x0,y0,0,'spline'))/0.0001%x=0处斜率 16 min2=min(interp1(x0,y0,13:0.001:15,'spline'))%13到15最小值 17 subplot(2,2,2); 18 plot(t,y2,'b',t,Y2,'r',x0,y0,'ro');%画出曲线 19 title('分段线性插值'); 20 legend('边条','线性');%显示图形图例 21 %三次线条插值A 22 y3=spline(x0,y0,t); 23 dy3=(spline(x0,y0,0.0001)-spline(x0,y0,0))/0.0001%x=0处斜率 24 min3=min(spline(x0,y0,13:0.001:15))%13到15最小值 25 subplot(2,2,3); 26 plot(x0,y0,'ro',t,y3);%画出曲线 27 title('三次线条插值A'); 28 %三次线条插值B 29 pp1=csape(x0,y0);%默认的边界条件,即给定边界一阶导数 30 pp2=csape(x0,y0,'second');%给定边界二阶导数 31 y4=ppval(pp1,t); 32 Y4=ppval(pp2,t); 33 dy4=(ppval(pp1,0.0001)-ppval(pp1,0))/0.0001%x=0处斜率 34 min4=min(ppval(pp1,13:0.001:15))%13到15最小值 35 subplot(2,2,4); 36 plot(t,y4,'b',t,Y4,'r',x0,y0,'ro');%画出曲线 37 title('三次线条插值B'); 38 legend('一阶','二阶');
综上，可以看出，拉格朗日插值函数根本不能应用，分段线性函数的光滑性较差，推荐三次样条插值。
同时，可以看出，interp1(x0,y0,’spline’)等价于spline(x0,y0)。
最后，将上述所有情况封装起来，变成下列函数：
1 function y = showAllInterp( x0,y0,s,t) 2 %显示x0,y0之间所有不同类型的插值情况 3 %字符串s选择要输出的插值类型： 4 %all：全部类型 lagrange：拉格朗日插值函数 5 %linear：分段线性插值 spline：三次线条插值A 6 %csape：三次线条插值B 7 if(nargin<4) 8 t=linspace(x0(1),x0(length(x0)),500);%默认 9 end 10 switch s 11 case 'lagrange' %拉格朗日插值函数 12 y1=lagrange(x0,y0,t);%调用编写的lagrange函数 13 plot(x0,y0,'ro',t,y1);%画出曲线 14 title('拉格朗日插值函数'); 15 if(nargout==1) 16 y=y1; 17 end 18 case 'linear' %分段线性插值 19 y2=interp1(x0,y0,t);%默认linear 20 plot(x0,y0,'ro',t,y2,'b');%画出曲线 21 title('分段线性插值'); 22 if(nargout==1) 23 y=y2; 24 end 25 case 'spline' %三次线条插值A 26 y3=spline(x0,y0,t); %等价于interp1(x0,y0,t,'spline'); 27 plot(x0,y0,'ro',t,y3);%画出曲线 28 title('三次线条插值A'); 29 if(nargout==1) 30 y=y3; 31 end 32 case 'csape' %三次线条插值B 33 pp1=csape(x0,y0);%默认的边界条件,即给定边界一阶导数 34 pp2=csape(x0,y0,'second');%给定边界二阶导数 35 y4=ppval(pp1,t); 36 Y4=ppval(pp2,t); 37 plot(t,y4,'b',t,Y4,'r',x0,y0,'ro');%画出曲线 38 title('三次线条插值B'); 39 legend('一阶','二阶'); 40 if(nargout==1) 41 y=y4; 42 end 43 case 'all' %显示全部 44 y1=lagrange(x0,y0,t);%调用编写的lagrange函数 45 subplot(2,2,1); 46 plot(x0,y0,'ro',t,y1);%画出曲线 47 title('拉格朗日插值函数'); 48 49 y2=interp1(x0,y0,t);%默认linear 50 subplot(2,2,2); 51 plot(x0,y0,'ro',t,y2);%画出曲线 52 title('分段线性插值'); 53 54 y3=spline(x0,y0,t); %等价于interp1(x0,y0,t,'spline'); 55 subplot(2,2,3); 56 plot(x0,y0,'ro',t,y3);%画出曲线 57 title('三次线条插值A'); 58 59 pp1=csape(x0,y0);%默认的边界条件,即给定边界一阶导数 60 pp2=csape(x0,y0,'second');%给定边界二阶导数 61 y4=ppval(pp1,t); 62 Y4=ppval(pp2,t); 63 subplot(2,2,4); 64 plot(t,y4,'b',t,Y4,'r',x0,y0,'ro');%画出曲线 65 title('三次线条插值B'); 66 legend('一阶','二阶'); 67 end
5、二维插值之插值节点为网格节点
已知m x n个节点：（xi，yj，zij）(i=1…m,j=1…n)，且xi,yi递增。求（x，y）处的插值z。 Matlab可以直接调用interp2(x0,y0,z0,x,y,`method`)其中 x0,y0 分别为 m 维和 n 维向量，表示节点， z0 为 n × m 维矩阵，表示节点值， x,y为一维数组，表示插值点， x 与 y 应是方向不同的向量，即一个是行向量，另一个是列向量， z 为矩阵，它的行数为 y 的维数，列数为 x 的维数，表示得到的插值， 'method'的用法同上面的一维插值。
如果是三次样条插值，可以使用命令pp=csape({x0,y0},z0,conds,valconds), z=fnval(pp,{x,y})其中 x0,y0 分别为 m 维和 n 维向量， z0 为 m × n 维矩阵， z 为矩阵，它的行数为 x 的维数，列数为 y 的维数，表示得到的插值，具体使用方法同一维插值。
eg:
1 x=100:100:500; 2 y=100:100:400; 3 z=[636 697 624 478 450 712 630 478 420 674 598 412 400 626 552 334 310]; 4 p=100:1:500; 5 q=100:1:400; 6 q=q';%须为列向量 7 z0=interp2(x,y,z,p,q);%分段线性插值 8 z1=interp2(x,y,z,p,q,'spline');%三次线条插值 9 subplot(2,1,1); 10 mesh(p,q,z0); 11 title('分段线性插值'); 12 subplot(2,1,2); 13 mesh(p,q,z1); 14 title('三次线条插值'); 15 %可以观察出，三次线条插值的图像更平滑
用csape函数插值：
1 x=100:100:500; 2 y=100:100:400; 3 z=[636 697 624 478 450 4 712 630 478 420 5 674 598 412 400 6 626 552 334 310]; 7 p=100:1:500; 8 q=100:1:400; 9 q=q'; 10 %三次线条插值 11 pp=csape({x,y},z');%注意跟interp2的区别,有个转置 12 z0=fnval(pp,{p,q}); 13 mesh(p,q,z0');%注意跟interp2的区别,有个转置 14 title('三次线条插值');
6、二维插值之插值节点为散乱节点
已知 n 个节点： ( xi , yi , zi )(i = 1,2,…, n) ，求点 (x, y) 处的插值 z 。对上述问题， Matlab 中提供了插值函数 griddata，其格式为： ZI = GRIDDATA(X,Y,Z,XI,YI)其中 X、 Y、 Z 均为 n 维向量，指明所给数据点的横坐标、纵坐标和竖坐标。向量 XI、YI 是给定的网格点的横坐标和纵坐标，返回值 ZI 为网格（ XI， YI）处的函数值。 XI与 YI 应是方向不同的向量，即一个是行向量，另一个是列向量。
eg:
1 %散乱节点的二维插值 2 x=[129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5]; 3 y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5]; 4 z=-[4 8 6 8 6 8 8 9 9 8 8 9 4 9]; 5 x0=[75:1:200]; 6 y0=[-85:1:145]'; 7 z0=griddata(x,y,z,x0,y0,'cubic');%保凹凸性3次插值 8 %[xx,yy]=meshgrid(x0,y0);无需采样，故不需要该函数 9 mesh(x0,y0,z0);
在上述问题中,补上寻找最大值的程序:
%max(z0)返回一个行向量，向量的第i个元素是矩阵A的第i列上的最大值 %find(A) 寻找矩阵A非零元素下标,返回矩阵A中非零元素所在位置 %[i,j,v]=find(A)返回矩阵A中非零元素所在的行i,列j,和元素的值v(按所在位置先后顺序输出) [p,q]=find(z0==max(max(z0))); zmax=z0(p,q)
三、最小二乘法实现曲线拟合
（1）用最小二乘法求一个形如 y = a + bx^ 2 的经验公式：
1 %等价于[1,x^2][a;b]=y,转换成解超定方程问题，超定方程的解是根据最小二乘法得来的 2 3 x=[19 25 31 38 44]'; 4 y=[19.0 32.3 49.0 73.3 97.8]'; 5 r=[ones(5,1),x.^2] 6 ab=r\y 7 x0=19:0.1:44; 8 y0=ab(1)+ab(2)*x0.^2; 9 plot(x,y,'o',x0,y0,'r')
（2）多项式拟合
a=polyfit(x,y,n)用多项式求过已知点的表达式，其中x为源数据点对应的横坐标，可为行向量、矩阵，y为源数据点对应的纵坐标，可为行向量、矩阵，n为你要拟合的阶数，一阶直线拟合，二阶抛物线拟合，并非阶次越高越好，看拟合情况而定，a为m+1的行向量。 polyfit函数的数学基础是最小二乘法曲线拟合原理，所得到的函数值在基点处的值与原来点的坐标偏差最小，常用于数据拟合，polyfit 做出来的值从左到右表示从高次到低次的多项式系数。如果要求拟合函数在x`点的函数值，可以调用polyval(a,x`)函数
eg：
1 x0=[1990 1991 1992 1993 1994 1995 1996]; 2 y0=[70 122 144 152 174 196 202]; 3 %画出散点图 4 plot(x0,y0,'ro'); 5 hold on 6 %用线性拟合 7 p=polyfit(x0,y0,1); 8 z0=polyval(p,x0); 9 plot(x0,z0);
四、最小二乘优化 (最小二乘:least square)
1 <span style="color:#000000;">%拟合形如y=a+bx^2的函数 2 %采样点 3 x=[19 25 31 38 44]'; 4 y=[19 32.3 49 73.3 97.8]'; 5 r=[ones(5,1),x.^2]; 6 ab=lsqlin(r,y) 7 x0=19:0.1:44; 8 y0=ab(1)+ab(2)*x0.^2; 9 plot(x,y,'o',x0,y0,'r')</span>
五、曲线拟合与函数逼近
eg：
求 f(x) =cos x, （-pi/2<=x<=pi/2） 在H = Span{1, x^2 , x^4} 中的最佳平方逼近多项式。
程序如下：
1 syms x%定义符号数值 2 base=[1,x^2,x^4]; 3 y1=base.'*base 4 y2=cos(x)*base.' 5 r1=int(y1,-pi/2,pi/2) 6 r2=int(y2,-pi/2,pi/2) 7 a=r1\r2%a为符号数值 8 xishu1=double(a)%化简符号数值 9 digits(8)%设置符号数值的精度 10 xishu2=vpa(a)%任意精度（符号类）数值
