第八章插值 拟合与回归.docx

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

2、xi, method, extrap)该命令用于找出由参量x决定的一元函数y=y(x)在点xi处的值yi。其中x,y为插值节点的横坐标和纵坐标，yi为在被插值点xi处的插值结果；x,y为向量， method表示采用的插值方法，MATLAB提供的插值方法有几种：nearest：最近邻点插值算法；linear：线性插值（默认）；spline：三次样条函数插值；pchip：分段三次Hermite插值；cubic：与pchip操作相同；v5cubic：在MATLAB 5.0中的三次插值。缺省时表示线性插值。 extrap表示对于超出x 范围的xi 中的分量将执行特殊的外插值法extrap。注意：所有的

3、插值方法都要求x是单调的。例1 对函数从0到10，步长为1进行采样得到插值节点，利用线性插值给出步长为0.25的插值函数图形。在命令窗口输入： x = 0:10;y = x.*sin(x); xx = 0:.25:10; yy = interp1(x,y,xx); plot(x,y,kd,xx,yy)输出结果如图8.1所示。例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 2

4、8 27 25 20 18 15 13;y1=interp1(x,y,a,spline)输出结果为：y1 = 27.8725若要得到一天24小时的温度曲线，则输入：xi=0:1/3600:24;yi=interp1(x,y,xi, spline);plot(x,y,o ,xi,yi)温度函数曲线图如图8.2所示。 图8.1 线性插值函数曲线图 图8.2 温度曲线图 8.1.2 interp2函数MATLAB中用于二维数据插值的函数是interp2，其调用格式为：zi = interp2(x, y, z, xi, yi, method)该命令用于找出由参量x,y决定的二元函数z=z(x,y)在点

5、（xi,yi）处的值zi。其中返回矩阵为zi，其元素为对应于参量xi与yi（可以是向量、或同型矩阵）的元素，若xi与yi 中有在x与y范围之外的点，则相应地返回NaN（Not a Number），method和interp1一样，常用的是cubic （双三次插值），缺省为linear（双线性插值算法）。例3 在命令窗口输入X,Y = meshgrid(-3:.25:3);Z = peaks(X,Y); %具有两个变量的采样函数，可产生一个凹凸有致的曲面，包含了三个局部极大点及三个局部极小点。XI,YI = meshgrid(-3:.125:3);ZZ = interp2(X,Y,Z,XI,YI

6、);surf(X,Y,Z);hold on;surf(XI,YI,ZZ+15) %为作比较，将插值曲面向上平移15单位。hold off输出结果如图8.3图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为

7、同型矩阵。其中y1,y2,y3为用命令meshgrid(XI,YI,ZI)生成的同型阵列。若插值点(xi,yi,zi)中有位于点(x,y,z)之外的点，则相应地返回特殊变量值NaN。method和interp1一样。例4 三维插值举例命令窗口输入 x,y,z,v = flow(20);xx,yy,zz = meshgrid(.1:.25:10, -3:.25:3, -3:.25:3);vv = interp3(x,y,z,v,xx,yy,zz);slice(xx,yy,zz,vv,6 9.5,1 2,-2 .2); shading interp;colormap cool输出图形如图8.4所示

8、。 图8.4 三维插值结果8.1.4 griddata函数 griddata也是一种常用的二维插值方法，其调用格式为：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（线性插值，为缺省算法）；natura

9、l（自然邻点插值），nearest（最近邻点插值）； cubic（ 基于三角形的三次插值）和v4（Matlab4中的griddata 算法）。例5 在命令窗口输入：rand(seed,0)x = rand(100,1)*4-2; y = rand(100,1)*4-2;z = x.*exp(-x.2-y.2);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 griddata插值结果例6

10、 有一组散乱数据点矩阵如下：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.439,3.970,

11、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;x=A(:,1)

12、;y=A(:,2);z=A(:,3);X,Y,Z=griddata(x,y,z,linspace(min(x),max(x),20),linspace(min(y),max(y),20),v4);mesh(X,Y,Z); hold onplot3(x,y,z,o); hold off;输出结果如图8.6所示。图8.6 散乱数据插值图8.1.5 spline函数该函数是利用三次样条对数据进行插值，其调用格式为：yy = spline(x,y,xx)该命令用三次样条插值计算出由向量x与y 确定的一元函数y=f(x)在点xx 处的值。若参量y 是一矩阵，则以y 的每一列和x 配对，再分别计算由它们确

13、定的函数在点xx 处的值。则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).*sin(x);xx = 0:.25:20;yy = spline(x,y,xx);plot(x,y,o,xx,yy) 输出结果如图8.7。图8.7 三次样条插值8.2 拟合拟合就是用连续曲线近似地刻画或比拟平面上离散点组所表示的坐标之间的函数关系的一种数据处理方法。如果已知某函数的若干离散函数值f1,f2,fn，通

14、过调整该函数中若干待定系数f(1, 2,n)，使得该函数与已知点集的差别（最小二乘意义）最小。如果待定函数是线性，就叫线性拟合或者线性回归(主要在统计中)，否则叫做非线性拟合或者非线性回归。MATLAB中提供了线性最小二乘拟合和非线性最小二乘拟合函数。8.2.1 polyfit函数ployfit函数是多项式拟合函数，其调用格式为：p = polyfit(x,y,n)其中x,y为长度相同的向量，n为拟合多项式的次数，返回值p是拟合多项式的系数向量，幂次由高到低，拟合多项式在x处的值可以通过y=polyval(p,x)来计算。例8，2004年全国大学生数学建模竞赛C题，饮酒驾车问题，时间和血液中酒

15、精浓度的函数关系，可以利用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次多项式拟合f=polyval(p,time); %求多项式的值hold on plot(time,f)hold off输出7多项次的系数向量为：p = 0.0002

16、 -0.0113 0.2817 -3.6307 25.5345 -94.6710 153.8757 -1.2592拟合曲线如图8.8所示。图8.8 血液酒精浓度拟合曲线MATLAB提供了两个非线性最小二乘拟合函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的。8.2.2 lsqcurvefit函数该函数用来进行非线性拟合，其调用格式为：x = lsqcurvefit (fun,x0,xdata,ydata,options);其中，fun为事先建立的拟合函数F(x,xdata)，其中自变量x表示拟合函数

17、中的待定参数，xdata为已知拟合节点的x坐标，x0为待定参数x的迭代初始值，xdata,ydata为已知数据点的x和y坐标，options是一些控制参数。lsqcurvefit函数用来求含参数x（向量）的向量值函数F(x,xdata)=f(x,xdata1), f(x,xdata2), f(x,xdatan)中的参数x（向量），使得最小。例9 根据表1中的数据，利用lsqcurvefit函数拟合。表1 已知数据点X1002003004005006007008009001000y1034.544.995.355.655.906.106.266.396.506.59首先建立拟合函数的M文件myf

18、it1.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。然后在命令窗口输入： 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 = lsqcurvefit (myfit1,x0,xdata,ydata)f= myfit1(x,xdata) 输出结果为：x = 0.0069

19、-0.0029 0.0809f = Columns 1 through 9 0.0044 0.0048 0.0051 0.0054 0.0056 0.0058 0.0060 0.0061 0.0063 Column 10 0.0064 %拟合函数在数据点的函数值拟合曲线如图8.9。 plot(xdata,ydata,o,xdata,f)图8.9 拟合曲线8.2.3 lsqnonlin函数该函数用来进行非线性拟合，其调用格式为：x = lsqnonlin (fun,x0, options);其中fun为事先建立的拟合函数F(x)，其中自变量x表示拟合函数中的待定参数，x0为待定参数x的迭代初始值

20、，options是一些控制参数。由于lsqnonlin中定义的拟合函数的自变量是x，所以已知参数xdata,ydata应写在该函数中。lsqnonlin函数用来求含参量x（向量）的向量值函数f(x)=f1(x), f2(x),，fn(x)中的参量x，使得最小。其中fi(x)=f(x,xdatai,ydaytai)=F(x, xdatai)- ydatai。例10 根据表1中的数据，利用lsqnonlin函数拟合。首先建立拟合函数的M文件myfit2.m，其内容如下：function f=myfit2(x)xdata=100:100:1000;ydata=1e-03*4.54,4.99,5.35

21、,5.65,5.90,6.10,6.26,6.39,6.50,6.59;f=x(1)+x(2)*exp(-0.02*x(3)*xdata)- ydata;然后在命令窗口输入： x0=0.2,0.05,0.05;x = lsqnonlin(myfit2,x0)f = myfit2(x)输出结果为：x = 0.0069 -0.0029 0.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.2436 Column 10 -0.2313该结果与

22、lsqcurvefit函数拟合的结果相同。8.3 回归分析MATLAB中提供了一些线性和非线性回归分析函数。8.3.1 regress函数一般地，称由确定的模型为一元线性回归模型，记为固定的未知参数、称为回归系数，自变量x也称为回归变量。一般地，称为高斯马尔柯夫线性模型(k元线性回归模型)。，MATLAB提供了多元线性回归函数regress，采用的是最小二乘估计，其调用格式有：b = regress (y,x)返回值为线性模型y = x*b的回归系数向量。其中x为n(k+1) 矩阵，行对应于观测值，列对应于预测变量，y为n维向量，为因变量，一元线性回归可取k=1。b,bint,r,rint,s

23、tats = regress(y,x,alpha)其中bint是回归系数的区间估计，r是残差，rint是置信区间，stats是用于检验回归模型的统计量，有四个数值：相关系数r2，F值，与F对应的概率P和误差方差估计，alpha是显著性水平(缺省的时候为0.05)。相关系数r2越接近1，说明回归方程越显著；与F对应的概率P y=8.8818 8.9487 9.0541 9.1545 9.2693 9.4289 9.6160 9.8150 9.9825 10.1558 10.3193;k= 7.8381 7.9167 8.0048 8.1026 8.2556 8.5822 8.8287 9.075

24、6 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,r,rint,stats=regress(y,x)输出结果为b = -55.4988 0.5644 7.1254 0.0222bint = -83.6008 -27

25、.3968 0.4693 0.6594 3.5740 10.6769 -0.1914 0.2359r = -0.0262 -0.0032 0.0033 0.0261 0.0164 -0.0249 0.0006 0.0124 0.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.0190 0.0438 -0.0304 0.0577 -0.0460 0.0259 -0.0420 0

26、.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，由于p rcoplot(r,rint)输出残差图见8.10。图8.10 残差图 8.3.2 rstool函数该函数是多元二项式回归函数，其调用格式为rstool(x,y, model,alpha)其中x为nm为矩阵，y为n维列向量，model为以下4种模型：linear（线性，缺省）：；interaction（交叉）：；quadratic（完全二次）：；pu

27、requadratic（纯二次）：。alpha为显著性水平，默认值为0.05。例13 设某商品的需求量与消费者的平均收入、商品价格的统计数据如表2，建立多元二项式纯二次回归模型，并预测平均收入为1000、价格为6时的商品需求量。表2 需求量、平均收入和价格统计表需求量10075807050659010011060收入10006001200500300400130011001300300价格5766875439可以直接使用多元二项式回归，在命令窗口输入： x1=1000 600 1200 500 300 400 1300 1100 1300 300;x2=5 7 6 6 8 7 5 4 3 9;

28、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。单击图形框左边Export，则出现图8.12对话框，可以将回归参数beta、剩余标准差rmse和残差residuals传送到MATALB的工作区中。图8.12输出对话框在命令窗口输入： betarmseres

29、iduals输出结果为：beta = 110.5313 0.1464 -26.5709 -0.0001 1.8475rmse = 4.5362residuals = 5.2724 -0.7162 -4.5158 -1.9390 -3.3315 3.4566 3.4843 -3.4452 -0.0976 1.8320故回归模型为：，剩余标准差为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=ones(10,1) x1 x2 (x1.2) (x2.2);b,bint,r,rint,stats=regress(y,X);b,stats输出结果为： b = 110.5313 0.1464 -26.5709