|
一、回归、插值、逼近、拟合的区别 1、回归一般指线性回归，是求最小二乘解的过程。在求回归前，已经假设所有型值点同时满足某一曲线方程，计算只要求出该方程的系数
2、多项式插值：用一个多项式来近似代替数据列表函数，并要求多项式通过列表函数中给定的数据点。（插值曲线要经过型值点。）
3、多项式逼近：为复杂函数寻找近似替代多项式函数，其误差在某种度量意义下最小。（逼近只要求曲线接近型值点，符合型值点趋势。）
4、多项式拟合：在插值问题中考虑给定数据点的误差，只要求在用多项式近似代替列表函数时，其误差在某种度量意义下最小。
注意：
表列函数：给定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称为插值点。
插值和拟合都是函数逼近或者数值逼近的重要组成部分
他们的共同点都是通过已知一些离散点集M上的约束，求取一个定义 在连续集合S(M包含于S)的未知连续函数，从而达到获取整体规律的目的，即通过'窥几斑'来达到'知全豹'。
简单的讲，所谓拟合是指已知某函数的若干离散函数值{f1,f2,…,fn}，通 过调整该函数中若干待定系数f(λ1, λ2,…,λ3), 使得该函数与已知点集的 差别(最小二乘意义)最小。如果待定函数是线性，就叫线性拟合或者线性回归(主要在统计中)，否则叫作非线性拟合或者非线性回归。表 达式也可以是分段函数，这种情况下叫作样条拟合。而插值是指已知某函数的在若干离散点上的函数值或者导数信息，通过求解该函数中待定形式的插值函数以及待定系数，使得该函数在给 定离散点上满足约束。插值函数又叫作基函数，如果该基函数定义在整个定义域上，叫作全域基，否则叫作分域基。如果约束条件中只有函数值的约束，叫作Lagrange插值，否则叫作Hermite插值。
从几何意义上将，拟合是给定了空间中的一些点，找到一个已知形式未知参数的连续曲面来最大限度地逼近这些点；而插值是找到一个(或几个分片光滑的)连续曲面来穿过这些点。
二、经典插值方法
1、拉格朗日插值
function y=lagrange(x0,y0,x) %拉格朗日插值函数 %n 个节点数据以数组 x0, y0 输入(注意 Matlat 的数组下标从1开始), %m 个插值点以数组 x 输入,输出数组 y 为 m 个插值 n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end
应用实例：x0=1:1:20;y0=x0.^2-20*x0-5;x=1:0.1:20;z=lagrange(x0,y0,x);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点）时的温度． x=0:2:24; y=[12,9,9,10,18,24,28,27,25,20,18,15,13]; x1=0:0.5:24; y1=interp1(x,y,x1,'linear'); plot(x,y,'bo',x1,y1,'r:');
3、埃米尔特插值 如果要求插值函数不仅在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至高阶导数值，这就是埃尔米特插值问题。
已知f（x）的n+1个节点的函数值f（xi）以及导数值f`（xi），可得一个至多n+1次的多项式H（x），即hermite插值多项式。新建以下这个函数： function y = hermite( x0,y0,y1,x ) %埃尔米特插值多项式 %x0为点横坐标 %y0为函数值 %y1为导数值 %m个插值点用数组x输入 n=length(x0);m=length(x); for k=1:m yy=0.0; for i=1:n h=1.0; a=0.0; for j=1:n if j~=i h=h*((x(k)-x0(j))/(x0(i)-x0(j)))^2; a=1/(x0(i)-x0(j))+a; end end yy=yy+h*((x0(i)-x(k))*(2*a*y0(i)-y1(i))+y0(i)); end y(k)=yy; end
4、样条插值 所谓样条（ Spline）本来是工程设计中使用的一种绘图工具，它是富有弹性的细木条或细金属条。绘图员利用它把一些已知点连接成一条光滑曲线（称为样条曲线），并使连接点处有连续的曲率。数学上将具有一定光滑性的分段多项式称为样条函数。在实际中最常用的是二次样条函数和三次样条函数：
二次样条函数插值 
首先，我们注意到s2 (x)中含有 n + 2 个特定常数，故应需要 n + 2 个插值条件，因此，二次样条插值问题可分为两类：
（1）已知插值节点xi 和相应的函数值 yi (i = 0,1,…,n) 以及端点 x0 （或 xn ）处的导数值y'0（或y'n）
（2）已知插值节点xi 和相应的导数值 y'i (i = 0,1,…,n) 以及端点 x0 （或 xn ）处的函数值y0 （或yn ）
三次样条函数插值 
由于 s3 (x)中含有n + 3 个待定系数，故应需要 n + 3 个插值条件，已知插值节点xi 和相应的函数值 f(xi ) = yi (i = 0,1,…,n) ，这里提供了 n + 1 个条件，还需要 2 个边界条件。因此，三次样条插值问题可分为三类：
（1）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 三次样条插值函数。
（2）s'3 (a) = y'0 ,s'3 (b) = y'3 。特别地 y'n = y'n = 0 时，称为自然边界条件。
（3）s'3 ( a + 0) = s'3 ( b − 0), s'3 (a + 0) = s'3 (b − 0) ， (这里要求 s3 (a + 0) =s3 (b − 0) )此条件称为周期条件。
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、分段线性和三次样条三种插值方法计算。
编程实现： clear,clc x0=[0,3,5,7,9,11,12,13,14,15]; y0=[0,1.2,1.7,2.0,2.1,2.0,1.8,1.2,1.0,1.6]; t=0:0.05:15; %拉格朗日插值函数 y1=lagrange(x0,y0,t);%调用编写的lagrange函数 dy1=(lagrange(x0,y0,0.0001)-lagrange(x0,y0,0))/0.0001%x=0处斜率 min1=min(lagrange(x0,y0,13:0.001:15))%13到15最小值 subplot(2,2,1); plot(x0,y0,'ro',t,y1);%画出曲线 title('拉格朗日插值函数'); %分段线性插值 y2=interp1(x0,y0,t,'spline');%注意区分spline与linear Y2=interp1(x0,y0,t);%默认linear dy2=(interp1(x0,y0,0.0001,'spline')-interp1(x0,y0,0,'spline'))/0.0001%x=0处斜率 min2=min(interp1(x0,y0,13:0.001:15,'spline'))%13到15最小值 subplot(2,2,2); plot(t,y2,'b',t,Y2,'r',x0,y0,'ro');%画出曲线 title('分段线性插值'); legend('边条','线性');%显示图形图例 %三次线条插值A y3=spline(x0,y0,t); dy3=(spline(x0,y0,0.0001)-spline(x0,y0,0))/0.0001%x=0处斜率 min3=min(spline(x0,y0,13:0.001:15))%13到15最小值 subplot(2,2,3); plot(x0,y0,'ro',t,y3);%画出曲线 title('三次线条插值A'); %三次线条插值B pp1=csape(x0,y0);%默认的边界条件,即给定边界一阶导数 pp2=csape(x0,y0,'second');%给定边界二阶导数 y4=ppval(pp1,t); Y4=ppval(pp2,t); dy4=(ppval(pp1,0.0001)-ppval(pp1,0))/0.0001%x=0处斜率 min4=min(ppval(pp1,13:0.001:15))%13到15最小值 subplot(2,2,4); plot(t,y4,'b',t,Y4,'r',x0,y0,'ro');%画出曲线 title('三次线条插值B'); legend('一阶','二阶');
综上，可以看出，拉格朗日插值函数根本不能应用，分段线性函数的光滑性较差，推荐三次样条插值。
同时，可以看出，interp1(x0,y0,’spline’)等价于spline(x0,y0)。
最后，将上述所有情况封装起来，变成下列函数：
function y = showAllInterp( x0,y0,s,t)%显示x0,y0之间所有不同类型的插值情况%字符串s选择要输出的插值类型：%all：全部类型 lagrange：拉格朗日插值函数%linear：分段线性插值 spline：三次线条插值A%csape：三次线条插值Bif(nargin<4) t=linspace(x0(1),x0(length(x0)),500);%默认endswitch s case 'lagrange' %拉格朗日插值函数 y1=lagrange(x0,y0,t);%调用编写的lagrange函数 plot(x0,y0,'ro',t,y1);%画出曲线 title('拉格朗日插值函数'); if(nargout==1) y=y1; end case 'linear' %分段线性插值 y2=interp1(x0,y0,t);%默认linear plot(x0,y0,'ro',t,y2,'b');%画出曲线 title('分段线性插值'); if(nargout==1) y=y2; end case 'spline' %三次线条插值A y3=spline(x0,y0,t); %等价于interp1(x0,y0,t,'spline'); plot(x0,y0,'ro',t,y3);%画出曲线 title('三次线条插值A'); if(nargout==1) y=y3; end case 'csape' %三次线条插值B pp1=csape(x0,y0);%默认的边界条件,即给定边界一阶导数 pp2=csape(x0,y0,'second');%给定边界二阶导数 y4=ppval(pp1,t); Y4=ppval(pp2,t); plot(t,y4,'b',t,Y4,'r',x0,y0,'ro');%画出曲线 title('三次线条插值B'); legend('一阶','二阶'); if(nargout==1) y=y4; end case 'all' %显示全部 y1=lagrange(x0,y0,t);%调用编写的lagrange函数 subplot(2,2,1); plot(x0,y0,'ro',t,y1);%画出曲线 title('拉格朗日插值函数'); y2=interp1(x0,y0,t);%默认linear subplot(2,2,2); plot(x0,y0,'ro',t,y2);%画出曲线 title('分段线性插值'); y3=spline(x0,y0,t); %等价于interp1(x0,y0,t,'spline'); subplot(2,2,3); plot(x0,y0,'ro',t,y3);%画出曲线 title('三次线条插值A'); pp1=csape(x0,y0);%默认的边界条件,即给定边界一阶导数 pp2=csape(x0,y0,'second');%给定边界二阶导数 y4=ppval(pp1,t); Y4=ppval(pp2,t); subplot(2,2,4); plot(t,y4,'b',t,Y4,'r',x0,y0,'ro');%画出曲线 title('三次线条插值B'); legend('一阶','二阶');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:
x=100:100:500; y=100:100:400; z=[636 697 624 478 450 712 630 478 420 674 598 412 400 626 552 334 310]; p=100:1:500; q=100:1:400; q=q';%须为列向量 z0=interp2(x,y,z,p,q);%分段线性插值 z1=interp2(x,y,z,p,q,'spline');%三次线条插值 subplot(2,1,1); mesh(p,q,z0); title('分段线性插值'); subplot(2,1,2); mesh(p,q,z1); title('三次线条插值'); %可以观察出，三次线条插值的图像更平滑
(2)、用csape函数插值：
x=100:100:500; y=100:100:400; z=[636 697 624 478 450 712 630 478 420 674 598 412 400 626 552 334 310]; p=100:1:500; q=100:1:400; q=q'; %三次线条插值 pp=csape({x,y},z');%注意跟interp2的区别,有个转置 z0=fnval(pp,{p,q}); mesh(p,q,z0');%注意跟interp2的区别,有个转置 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:
%散乱节点的二维插值 x=[129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5]; y=[7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5]; z=-[4 8 6 8 6 8 8 9 9 8 8 9 4 9]; x0=[75:1:200]; y0=[-85:1:145]'; z0=griddata(x,y,z,x0,y0,'cubic');%保凹凸性3次插值 %[xx,yy]=meshgrid(x0,y0);无需采样，故不需要该函数 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,x^2][a;b]=y,转换成解超定方程问题，超定方程的解是根据最小二乘法得来的x=[19 25 31 38 44]'; y=[19.0 32.3 49.0 73.3 97.8]'; r=[ones(5,1),x.^2] ab=r\y x0=19:0.1:44; y0=ab(1)+ab(2)*x0.^2; 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：
x0=[1990 1991 1992 1993 1994 1995 1996]; y0=[70 122 144 152 174 196 202]; %画出散点图 plot(x0,y0,'ro'); hold on %用线性拟合 p=polyfit(x0,y0,1); z0=polyval(p,x0); plot(x0,z0);
四、最小二乘优化 (最小二乘:least square)
<span style='color:#000000;'>%拟合形如y=a+bx^2的函数 %采样点 x=[19 25 31 38 44]'; y=[19 32.3 49 73.3 97.8]'; r=[ones(5,1),x.^2]; ab=lsqlin(r,y) x0=19:0.1:44; y0=ab(1)+ab(2)*x0.^2; 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} 中的最佳平方逼近多项式。 程序如下： syms x%定义符号数值 base=[1,x^2,x^4]; y1=base.'*base y2=cos(x)*base.' r1=int(y1,-pi/2,pi/2) r2=int(y2,-pi/2,pi/2) a=r1\r2%a为符号数值 xishu1=double(a)%化简符号数值 digits(8)%设置符号数值的精度 xishu2=vpa(a)%任意精度（符号类）数值
转自http://www.cnblogs.com/vanker/p/4268292.html
|
