第八章插值拟合与回归

1、第八章 插值 拟合与回归在生产和实验中，由于函数 f (x) 的表达式不便于计算或者没有表达式而只 有在给定点的函数值(或其导数值) ，为此，我们希望建立一个简单的而且便于 计算的近似函数g(x)，来逼近函数f(x)，这就用到插值和拟合方法。8.1 插值 插值就是在离散数据的基础上补插连续函数，使得这条连续曲线 通过全部给定的 离散数据点 。插值是离散函数逼近的重要方法，利用它可通过函数在有 限个点处的取值状况，估算出函数在其他点处的近似值。MATLAB 中的插值函数主要有以下几个。8.1.1 interp1 函数MATLAB 中用于一维数据插值的函数是 interp1 ，其调用格式为：yi

2、= interp1(x, y, xi, 'method, 'extrap')该命令用于找出由参量x决定的一元函数y=y(x)在点xi处的值yi。其中x,y 为插值节点的横坐标和纵坐标， yi 为在被插值点 xi 处的插值结果； x,y 为向量， method '表示采用的插值方法，MATLAB提供的插值方法有几种：neares：t '最近邻点插值算法；'linear：线性插值(默认)； spline三次样条函数插值；'pchip：'分段三次Hermite插值；'cubic '与'pchip操作相同；'

3、;v5cubic在MATLAB 5.0 中的三次插值。 缺省时表示线性插值。'extrap '表示对于超出 x 范围的 xi 中的分量将执行特殊的外插值法 extrap。 注意：所有的插值方法都要求x是单调的。例 1 对函数 y xsin x 从 0 到 10，步长为 1 进行采样得到插值节点， 利用线 性插值给出步长为 0.25的插值函数图形。在命令窗口输入：>> x = 0:10;y = x.*s in( x);xx = 0:.25:10;yy = in terp1(x,y,xx);plot(x,y,'kd',xx,yy)输出结果如图8.1所示。

4、例2在一天24小时内，从零点开始每间隔2小时测得的环境温度数据分别 为12，9，9，1，0, 18，24，28，27，25，20，18，15，13推测中午13时的温度。我们采用三次样条插值，在命令窗口输入：>> x=0:2:24;a=13;y=12 9 9 10 18 24 28 27 25 20 18 15 13;y1=i nterp1(x,y,a,'spli ne')输出结果为：y1 =27.8725若要得到一天24小时的温度曲线，则输入：xi=0:1/3600:24;yi=in terp1(x,y,xi, 'spli ne');plot(x,y

5、,'o' ,xi,yi)温度函数曲线图如图8.2所示。-4图8.2温度曲线图301“1125-J20j|!I115M-1 /10-r5图8.1线性插值函数曲线图8.1.2 interp2 函数MATLAB中用于二维数据插值的函数是interp2，其调用格式为：zi = in terp2(x, y, z, xi, yi, method '该命令用于找出由参量x,y决定的二元函数z=z(x,y)在点(xi,yi )处的值zi。 其中返回矩阵为zi，其元素为对应于参量xi与yi (可以是向量、或同型矩阵) 的元素，若xi与yi中有在x与y范围之外的点，贝U相应地返回NaN (

6、 Not aNumber)，method和interp1 一样，常用的是 cubic'(双三次插值)，缺省为 linear (双线性插值算法)。例3在命令窗口输入X,Y = meshgrid(-3:.25:3);Z = peaks(X,Y);%具有两个变量的采样函数，可产生一个凹凸有致的曲面，包含了三个局部极大点及三个局部极小点。XI,YI = meshgrid(-3:.125:3);ZZ = in terp2(X,Y,Z,XI,YI);surf(X,Y,Z);hold on;surf(XI,YI,ZZ+15)%为作比较，将插值曲面向上平移 15单位hold off输出结果如图8.32

7、0020-20-4-430 10 .4图8.3二维插值曲面8.1.3 interp3 函数MATLAB中用于三维数据插值的函数是interp3，其调用格式为:vi = interp3(x, y, z, v, xi, yi, zi, method '该命令用于找出由参量x,y,z决定的三元函数v=v(x,y,z)在点(xi,yi,zi )处的值vi。参量xi,yi,zi是同型阵列或向量。若向量参量xi,yi,zi是不同长度，不同方向（行或列）的向量，这时输出参量 vi与y1,y2,y3为同型矩阵。其中y1,y2,y3 为用命令meshgrid（XI,YI,ZI）生成的同型阵列。若插值点（

8、xi,yi,zi）中有位于点(x,y,z )之外的点，则相应地返回特殊变量值NaN。method和in terpl 一样。例4三维插值举例命令窗口输入x,y,z,v = flow(20);xx,yy,zz = meshgrid(.1:.25:10, -3:.25:3, -3:.25:3);vv = in terp3(x,y,z,v,xx,yy,zz);slice(xx,yy,zz,vv,6 9.5,1 2,-2 .2);shad ing in terp;colormap cool输出图形如图8.4所示。-4010图8.4三维插值结果8.1.4 griddata 函数griddata也是一种常用

9、的二维插值方法，其调用格式为：zi =griddata （x, y, z, xi, yi, method '该命令用于找出由参量x,y决定的二元函数z=z（x,y）在点（xi,yi）处的值zi。 它和interp2的区别在于，interp2的插值数据必须是 矩形域，即已知数据点（x,y） 组成规则的矩阵，可使用meshgid生成。而griddata函数的已知数据点（x,y） 不要求规则排列，特别是对试验中随机没有规律采取的数据进行插值具有很好 的效果。method包括：linear (线性插值，为缺省算法)；natural '(自然邻点插 值)，nearest('最近邻

10、点插值);cubic (基于三角形的三次插值)和 v4'(Matlab4 中的 griddata 算法)。例5在命令窗口输入：ran d('seed',0)x = ran d(100,1)*4-2; y = ran d(100,1)*4-2;z = x.*exp(-x.A2-y.A2);ti = -2:.125:2;XI,YI = meshgrid(ti,ti);ZI = griddata(x,y,z,XI,YI);mesh(XI,YI,ZI), hold;plot3(x,y,z,'o'), hold off;输出结果如图8.5所示。图8.5 gridd

11、ata插值结果例6有一组散乱数据点矩阵如下：A=1.486,3.059,0.1;2.121,4.041,0.1;2.570,3.959,0.1;3.439,4.396,0.1;4.505,3.012,0.1;3.402,1.604,0.1;2.570,2.065,0.1;2.150,1.970,0.1;1.794,3.059,0.2;2.121,3.615,0.2;2.570,3.473,0.2;3.421,4.160,0.2;4.271,3.036,0.2;3.411,1.876,0.2;2.561,2.562,0.2;2.179,2.420,0.2;2.757,3.024,0.3;3.43

12、9,3.970,0.3;4.084,3.036,0.3;3.402,2.077,0.3;2.879,3.036,0.4;3.421,3.793,0.4;3.953,3.036,0.4;3.402,2.219,0.4;3.000,3.047,0.5;3.430,3.639,0.5;3.822,3.012,0.5;3.411,2.385,0.5;3.103,3.012,0.6;3.430,3.462,0.6;3.710,3.036,0.6;3.402,2.562,0.6;3.224,3.047,0.7;3.411,3.260,0.7;3.542,3.024,0.7;3.393,2.763,0.7;

13、x=A(:,1);y=A(:,2);z=A(:,3);X,Y,Z=griddata(x,y,z,li nspace(mi n(x),max(x),20)',li nspace(mi n(y),max(y),20),'v4');mesh(X,Y,Z); hold onplot3(x,y,z,'o'); hold off;输出结果如图8.6所示。1 1图8.6散乱数据插值图8.1.5 spline 函数该函数是利用三次样条对数据进行插值，其调用格式为：yy = spli ne(x,y,xx)该命令用三次样条插值计算出由向量x与y确定的一元函数y=f(x)在点

14、xx处的值。若参量y是一矩阵，则以y的每一列和x配对，再分别计算由它们确 定的函数在点xx处的值。贝U yy是一个阶数为length(xx)*size(y,2)的矩阵。例7对离散地分布在y=exp(x)sin(x)函数曲线上的数据点进行样条插值。在命令窗口输入：x = 0 2 4 5 8 12 12.8 17.2 19.9 20;y = exp(x).*si n( x);xx = 0:.25:20;yy = spli ne(x,y,xx);plot(x,y,'o',xx,yy)输出结果如图8.7。8图8.7三次样条插值8.2拟合拟合就是用连续曲线近似地刻画或比拟平面上离散点组所

15、表示的坐标之间 的函数关系的一种数据处理方法。如果已知某函数的若干离散函数值 f1,f2,通过调整该函数中若干待定系数f(入U 2,使得该函数与已知点集的差别(最小二乘意义)最小。如果待定函数是线性，就叫线性拟合或 者线性回归(主要在统计中)，否则叫做非线性拟合或者非线性回归。MATLAB中提供了线性最小二乘拟合和非线性最小二乘拟合函数。8.2.1 polyfit 函数ployfit函数是多项式拟合函数，其调用格式为：p = polyfit(x,y, n)其中x,y为长度相同的向量，n为拟合多项式的次数，返回值 p是拟合多项式的 系数向量，幕次由高到低，拟合多项式在 x处的值可以通过y=pol

16、yval(p,x)来计 算。例8，2004年全国大学生数学建模竞赛 C题，饮酒驾车问题，时间和血液 中酒精浓度的函数关系，可以利用polyfit进行拟合。可输入：>> time=0.25, 0.5 0.75, 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 11 12 13 14 15 16' vol=20 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4' plot(time,vol,'o')p=polyfit(time,vol,7) % 7 次多

17、项式拟合f=polyval(p,time);% 求多项式的值hold onplot(time,f) hold off输出7多项次的系数向量为:0.0002 -0.0113 0.2817 -3.6307 25.5345 -94.6710 153.8757-1.2592拟合曲线如图8.8所示。908070 .60 -50 P40 .30 .20 r10 L02468101214160 L图8.8血液酒精浓度拟合曲线MATLAB提供了两个非线性最小二乘拟合函数：Isqcurvefit和lsqnonlin。 两个命令都要先建立 M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方 式是不

18、同的。8.2.2 lsqcurvefit 函数该函数用来进行非线性拟合，其调用格式为：x = lsqcurvefit ('fu n',x0,xdata,ydata,optio ns);其中，fun为事先建立的拟合函数F(x,xdata)，其中自变量x表示拟合函数中的 待定参数，xdata为已知拟合节点的x坐标，x0为待定参数x的迭代初始值， xdata,ydata为已知数据点的x和y坐标，options是一些控制参数。lsqcurvefit函数用来求含参数x (向量)的向量值函数F(x,xdata)=f(x,xdata 1), f(x,xdata 2),f(x,xdata n)

19、中的参数x (向量)，使得n2f (x,xdatai) ydataji 10.02kx最小。例9根据表1中的数据，利用lsqcurvefit函数拟合y( x) a be表1已知数据点X1002003004005006007008009001000y X1034.544.995.355.655.906.106.266.396.506.59首先建立拟合函数的 M文件myfitl.m，其内容如下:function f = myfit1(x,xdata)f = x(1)+x(2)*exp(-0.02*x (3)*xdata);其中x(1),x(2),x(3)分别表示拟合曲线中的参数a,b,k。然后在命

20、令窗口输入：>> xdata = 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.2,0.05,0.05;x = Isqcurvefit ('myfit1',x0,xdata,ydata)f= myfit1(x,xdata)输出结果为：x =0.0069-0.00290.08090.00440.0060 0.0061Colum n 100.0064拟合曲线如图8.9。Columns 1 through 90.00560.00580.00480.0051

21、0.00540.0063%拟合函数在数据点的函数值>> plot(xdata,ydata,'o',xdata,f)6.55.54.510004100图8.9拟合曲线8.2.3 Isqnonlin 函数该函数用来进行非线性拟合，其调用格式为：x = Isqnonlin ('fun',x0, options);其中fun为事先建立的拟合函数F(x)，其中自变量x表示拟合函数中的待定参数， x0为待定参数x的迭代初始值，options是一些控制参数。由于lsqnonlin中定义的拟合函数的自变量是x，所以已知参数xdata,ydata应写在该函数中。lsq

22、nonlin函数用来求含参量 x (向量)的向量值函数f(x)=f 1(x), f2(x),，fn(x)n中的参量 x，使得fi2(x)最小。其中fi(x)=f(x,xdata i,ydaytai)=F(x, xdatai)-i 1ydatai。例10根据表1中的数据，利用lsqnonlin函数拟合y(x) a be 0.02kx。首先建立拟合函数的 M文件myfit2.m，其内容如下：function f=myfit2(x)xdata=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;f=x(

23、1)+x(2)*exp(-0.02*x (3) *xdata)- ydata;然后在命令窗口输入：>> x0=0.2,0.05,0.05;x = Isqnon li n('myfit2',x0) f = myfit2(x)输出结果为:0.0069-0.00290.0809f =1.0e-03 *Columns 1 through 9-0.0971-0.1739-0.2164-0.2463-0.2666-0.2713-0.2651-0.2538-0.2436Colum n 10-0.2313该结果与Isqcurvefit函数拟合的结果相同。8.3回归分析MATLAB中

24、提供了一些线性和非线性回归分析函数。8.3.1 regress 函数般地，称由y 01x确定的模型为一元线性回归模型，记为y 01XE 0,D2固定的未知参数1称为回归系数，自变量X也称为回归变量。般地，称E( ) 0,COV(,)2In为高斯一马尔柯夫线性模型（k11元线性回归模型）y1X11屜1X12X22X1kX2kynXuXn2XnkMATLAB提供了多元线性回归函数regress，采用的是最小二乘估计， 其调用格式有:b = regress （y,x）返回值为线性模型y = x*b的回归系数向量。其中x为nx（k+1）矩阵，行对 应于观测值，列对应于预测变量，y为n维向量，为因变量，

25、一元线性回归可取 k=1。b,bint,r,rint,stats = regress(y,x,alpha)其中bint是回归系数的区间估计，r是残差，rint是置信区间，stats是用于检验 回归模型的统计量，有四个数值：相关系数 rA2，F值，与F对应的概率P和误 差方差估计，alpha是显著性水平(缺省的时候为0.05)。相关系数宀2越接近1, 说明回归方程越显著；与 F对应的概率Pvalpha时候拒绝H0，回归模型成立。例11输入>> y=8.8818 8.9487 9.0541 9.1545 9.2693 9.4289 9.6160 9.8150 9.9825 10.155

26、8 10.3193'k= 7.8381 7.9167 8.0048 8.1026 8.2556 8.5822 8.8287 9.0756 9.2175 9.4148 9.6198'l= 8.3871 8.3872 8.3935 8.3971 8.4025 8.4048 8.4079 8.4141 8.4261 8.4377 8.4444't= 9.9551 9.9057 10.0972 9.9537 9.9370 9.9449 9.9636 10.1291 10.1573 10.2944 10.2093'x= ones(size(k) k l t;b,bint,

27、r,rint,stats=regress(y,x)输出结果为b =-55.49880.56447.12540.0222bint =-83.6008 -27.39680.46930.65943.574010.6769-0.19140.2359-0.0262 -0.00320.00330.02610.0164-0.02490.00060.01240.0137-0.0101-0.0081rint =-0.0620 0.0097-0.0479 0.0415-0.0306 0.0372-0.0111 0.0632-0.0218 0.0546-0.0621 0.0124-0.0385 0.0397-0.0

28、190 0.0438-0.0304 0.0577-0.0460 0.0259-0.0420 0.0258stats =1.0e+003 *0.0010 2.1022 0.0000 0.0000 该结果说明 y=-55.4988+0.5644k+7.1254l+0.0222t，stats 中的数据说明 r2=1，F=2102.2, p=0，由于pv0.05可知回归模型成立。可以利用 rcoplot 函数画出残差及其置信区间，红色的表示超出期望值的数据，圆圈代表残差的值，竖线代表置信区间的范围输入>> rcoplot(r,ri nt)输出残差图见8.10。RResidual Case

29、Order Plot0.060.040.020-0.02-0.04-0.061234567891011Case Number图8.10残差图8.3.2 rstool 函数该函数是多元二项式回归函数，其调用格式为rstool（x,y, 'model',alpha）其中x为nXm为矩阵，y为n维列向量，'model'为以下4种模型:'linear'（线性，缺省）：y ° 必 LmXm;'interaction'（交叉）：y ° 必 LmXmjkXjXk;1 j k m'quadratic'（完全二次

30、）：y 01X1 LmXmjkXjXk ;1 j,k mn'purequadratic'（纯二次）：y 01x1 LmxmjjX2。j 1alpha为显著性水平，默认值为0.05。例13设某商品的需求量与消费者的平均收入、商品价格的统计数据如表2, 建立多元二项式纯二次回归模型，并预测平均收入为 1000、价格为6时的商品 需求量。表2需求量、平均收入和价格统计表需求量10075807050659010011060收入10006001200500300400130011001300300价格5766875439可以直接使用多元二项式回归，在命令窗口输入:>> x仁1

31、000 600 1200 500 300 400 1300 1100 1300 300;x2=5 7 6 6 8 7 5 4 3 9;y=100 75 80 70 50 65 90 100 110 60' x=x1' x2'rstool(x,y,'purequadratic')其输出如图8.11。图8.11多元二项式回归在图8.11中x1上面的方框中输入1000, x2上面的方框中输入6,在图形框左侧的“ Predicted Y1 ”下方的数据变为88.47981，即预测出平均收入为1000, 价格为6时的商品需求量为88.4791。单击图形框左边Exp

32、ort，则出现图8.12对话框，可以将回归参数 beta、剩余标准差rmse和残差residuals传送到MATALB的工作区中。图8.12输出对话框在命令窗口输入:>> betarmseresiduals输出结果为：beta =110.53130.1464-26.5709-0.00011.8475rmse =4.5362residuals =5.2724-0.7162-4.5158-1.9390-3.33153.45663.4843-3.4452-0.09761.8320故回归模型为： y 110.5313 0.1464x1 26.5709 x2 0.0001x12 1.8475

33、x22 ，剩 余标准差为 4.5362, 说明此回归模型的显著性较好。还可将该模型转化为多元线性回归模型，利用 regress 函数进求解。可输入：>> x1=1000 600 1200 500 300 400 1300 1100 1300 300;x2=5 7 6 6 8 7 5 4 3 9;X=o nes(10,1) x1' x2' (x"2)' (x2.A2)'b,bint,r,rint,stats=regress(y,X);b,stats输出结果为：b =110.53130.1464-26.5709-0.00011.8475stat

34、s =0.9702 40.6656 0.0005 20.5771 输出的回归参数和 rstool 函数的结果相同。8.3.3 nlinfit 函数和 nlintool 函数nlinfit 函数用来确定非线性回归系数，调用格式为：beta, r, J = nlinfit (x, y, 'modelfun', beta0)其中输入数据x,y分别为nxp维矩阵和n维列向量，对于一元非线性回归，取 p=1,即可；modelfun'为事先定义的非线性回归函数的M文件，是回归系数beta和 x 的函数； beta0 是回归系数的初值，输出参数 beta 是估计出的回归系数， r 为残差， J 为 Jacobain 矩阵。例 14 根据表 1 中的数据，利用 nlinfit 函数进行非线性回归，回归函数为y(x)a be 0.02kx首先建立回归函数的 M 文件 myfit3.m ，内容如下： function f = myfit3(beta,xdata)f = beta(1)+beta(2)*exp(-0.02*beta(3)*xdata); 在命令窗口输入：>> xdata = 100:100:1000;ydata = 1e-03*4.54,4.99,5.35,5.65