第六章 Matlab插值、拟合与回归

6.1.2 二维插值
Matlab提供了两个二维插值的函数interp2（对二维网格数据插值）和 griddata（对二维随机数据点的插值）。其中interp2的调用格式有：
格式1 ZI = interp2(X,Y,Z,XI,YI) 返回矩阵ZI，其元素包含对应于参量XI与YI（可以是向量、或同型矩阵）的元素。 参量X与Y必须是单调的，且相同的划分格式(与meshgrid相似）。若Xi与Yi中有在X
第六章 Matlab插值、拟合与回归
主要内容
6.1. 插值运算 6.2. 拟合运算 6.3. 回归分析
6.1. 插值运算
6.1.1一维插值
Matlab提供了interp1()函数进行一维多项式插值，该函数使用多项式
技术，用多项式函数通过所提供的数据点计算目标插值点上的插值函数值。
其调用方法有：
格式1 yi = interp1(x,Y,xi) %返回插值向量yi，每一元素对应于参量xi，同时由向量x与Y的内插值决 定。参量x指定数据Y的点。若Y为一矩阵，则按Y的每列计算。 格式2 yi = interp1(x,Y,xi,method) %用指定的算法计算插值： ’nearest’：最近邻点插值，直接完成计算； ’linear’：线性插值（缺省方式），直接完成计算； ’spline’：三次样条函数插值。 ’cubic’： 分段三次Hermite插值。 其它，如’v5cubic’ 。 对于超出x范围的xi的分量，使用方法’nearest’、’linear’、’cubic’的插 值算法，相应地将返回NaN。对其他的方法，interp1将对超出的分量执行外插值 算法。i = interp1(x,Y,xi,method,'extrap')
x = [0 2 4 5 8 12 12.8 17.2 19.9 20]; y = exp(x).*sin(x); xx = 0:.25:20; yy = spline(x,y,xx); plot(x,y,'o',xx,yy)
6.2. 拟合运算
在实际工程应用和科学实践中，经常需要寻求两个（或多个）变量间的关系，
interp2函数能够较好的进行二维插值运算，但是它只能处理以网格形
式给出的数据，如果数据是随机给出的，此时可以使用griddata函数。其 调用格式为：
格式:z=griddata(x0,y0,z0,x,y,’method’) ’ v4 ’:MATLAB4.0提供的插值算法，公认效果较好； ’linear’：双线性插值算法（缺省算法）； ’nearest’：最临近插值；
clear x=-3+6*rand(199,1);y=-2+4*rand(199,1); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); plot(x,y,'*');figure;plot3(x,y,z,'*');figure; [x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2); z1=griddata(x,y,z,x1,y1,'cubic'); surf(x1,y1,z1);figure; z2=griddata(x,y,z,x1,y1,'v4');surf(x1,y1,z2); %误差比较
clear;[x,y]=meshgrid(-3:0.6:3,-2:0.4:2); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); surf(x,y,z);figure [x1,y1]=meshgrid(-3:0.2:3,-2:0.2:2); z1=interp2(x,y,z,x1,y1); z2=interp2(x,y,z,x1,y1,'cubic'); z3=interp2(x,y,z,x1,y1,'spline'); surf(x1,y1,z1);figure %比较误差 z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1); surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.08]) figure;surf(x1,y1,abs(z0-z3)); axis([-3,3,-2,2,0,0.08])
与Y范围之外的点，则相应地返回NaN。
格式2 ZI = interp2(Z,XI,YI) 种情形进行计算。 格式3 ZI = interp2(X,Y,Z,XI,YI,method) 用指定的算法method计算二维插值： ’linear’：双线性插值算法（缺省）； ’nearest’：最临近插值； ’spline’：三次样条插值； ’cubic’：双三次插值。 表示X=1:n、Y=1:m，其中[m,n]=size(Z)。再按第一
小，在Matlab中使用polyfit()函数实现，调用格式为：p=polyfit(x,y,n)，n为选 定的多项式的次数，p是多项式按降幂排列得出的行向量，可以使用
poly2sym转换为多项式，使用polyval()计算多项式的值。
clear;clc; x0=-1+2*[0:10]/10; y0=1./(1+25*x0.^2); x=-1:0.01:1; ya=1./(1+25*x.^2); p3=polyfit(x0,y0,3);y1=polyval(p3,x); p5=polyfit(x0,y0,5);y2=polyval(p5,x); p8=polyfit(x0,y0,8);y3=polyval(p8,x);
示例1: years=1950:10:1990; service=10:10:30;
wage = [150.697 199.592 187.625; 179.323 195.072 250.287 203.212 179.092 322.767; 226.505 153.706 426.730 249.633 120.281 598.243]; w = interp2(service,years,wage,15,1975)
legend('位移曲线','速度曲线','位移点','速度点');
’cubic’：双三次插值。
还有三维插值运算函数interp3,n维网格插值interpn，其调用格式同 interp1和interp2，对应的三维网格生成函数为[x,y,z]=meshgrid(x1,y1,z1) 和非网格生成函数griddata3(),griddatan()，他们同griddata调用格式相同。
示例1:
x=0:0.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);
plot(x,y,'ro',x,y);
x1=0:0.02:1; %要插值点 y1=(x1.^2-3*x1+5).*exp(-5*x1).*sin(x1); y2=interp1(x,y,x1);y3=interp1(x,y,x1,'spline'); y4=interp1(x,y,x1,'nearest');y5=interp1(x,y,x1,'cubic'); plot(x1,[y2' y3' y4' y5'],':',x,y,'ro',x1,y1); legend('linear','spline','nearest','cubic','样本点','原函数'); %计算各插值的最大误差 [max(abs(y1-y2)),max(abs(y1-y3)),max(abs(y1-y4)),max(abs(y1y5))];
而实际去只能通过观测得到一些离散的数据点。针对这些分散的数据点，运
用某种拟合方法生成一条连续的曲线，这个过程称为曲线拟合。 曲线拟合可分为：
（1）参数拟合 — 最小二乘法
（2）非参数拟合 — 插值法
6.2.1 多项式拟合
一般多项式拟合的目标是寻找一组多项式系数ai,使得多项式
f(x)=a1xn+a2xn-1 +…+anx+an+1 能够较好的拟合原始数据，使整体拟合误差较
v=[0 .4794 .8415 .9975 .9093 .5985 .1411]; s=[1 1.5 2 2.5 3 3.5 4]; p1=polyfit(t,s,8); p2=polyfit(t,v,8); tt=0:0.1:3; s1=polyval(p1,tt); v1=polyval(p2,tt); plot(tt,s1,'r-',tt,v1,'b',t,s,'p',t,v,'d'); xlabel('t');ylabel('x(t),y(t)');
z0=(x1.^2-2*x1).*exp(-x1.^2-y1.^2-x1.*y1);
surf(x1,y1,abs(z0-z1));axis([-3,3,-2,2,0,0.15]) figure; surf(x1,y1,abs(z0-z2));axis([-3,3,-2,2,0,0.15])
[x,y,z]=meshgrid(-1:0.2:1); [x0,y0,z0]=meshgrid(-1:0.05:1); V=exp(x.^2.*z+y.^2.*x+z.^2.*y).*cos(x.^2.*y.*z+z.^2.*y.*x); V0=exp(x0.^2.*z0+y0.^2.*x0+z0.^2.*y0).*cos(x0.^2.*y0.*z0+z0.^2.*y0.*x0); V1=interp3(x,y,z,V,x0,y0,z0,'spline'); err=V1-V0; max(err(:)) slice(x0,y0,z0,V1,[-0.5,0.3, 0.9],[0.6,-0.1],[-1,-0.5,0.5,1]) title('Slives for Four Dim Figures');
p10=polyfit(x0,y0,10);y4=polyval(p10,x);
plot(x,ya,x,y1,x,y2,'-.',x,y3,'--',x,y4,':')
示例2:已知测得某质点的位移和随度时间变化如下，求该质点的速度与位
移随时间的变化曲线以及位移随速度的变化曲线。
clear;clc
t=[0 0.5 1 1.5 2 2.5 3];
6.1.3 三次样条插值
已知平面上n个点xi,yi(i=1...n)，其中x1<x2<…<xn，S(x)为三次样条 函数需要满足三个条件： (1)S(xi)=yi(i=1...n)，即函数经过样本点； (2)S(x)在每个子区间[xi,xi+1]上为三次多项式：
(3)S(x)在整个区间[x1,xn]上有连续的一阶及二阶导数。 Matlab提供了一个生成三次参数样条类的函数S=csapi(x,y) 其中x=[x1,x2,….,xn], y=[y1,y2,…,yn]为样本点。S返回样条函数对象的
x0=-3:.6:3; y0=-2:.4:2; [x,y]=ndgrid(x0,y0); z=(x.^2-2*x).*exp(-x.^2-y.^2-x.*y); sp=csapi({x0,y0},z); fnplt(sp);
示例7:对离散分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值计算
sp=csapi(x,y); fnplt(sp)
c=[sp.breaks(1:4)' sp.breaks(2:5)' sp.coefs(1:4,:),... sp.breaks(5:8)' sp.breaks(6:9)' sp.coefs(5:8,:)]
处理多个自变量的网格数据三次样条插值类的函数： S=csapi({x1,x2,…,xn},z) 其中xi为自变量标识，z是网格数据的样本点，S是三次样条曲线的对象。 函数spline可以进行三次样条数据插值 格式 yy = spline(x,y,xx)
插值结果，包括子区间点、各区间点三次多项式系数等。
可用 fnplt()绘制出插值结果，其调用格式：fnplt(Hale Waihona Puke Baidu) 对给定的向量xp，可用fnval()函数计算，其调用格式：yp=fnval(S,xp)
其中得出的yp是xp上各点的插值结果。
x=0:.12:1; y=(x.^2-3*x+5).*exp(-5*x).*sin(x);