第七章-数据处理.doc

_第七章 数据处理及应用【教学目的】：本章着重介绍数据的插值、拟合与多元回归分析，并利用Matlab相应工具结合具体案例进行分析。【教学重点难点】：教学重点：插值与拟合，回归分析，聚类与分类。教学难点：定性与定量的结合，如何在定性分析基础上给出适当的定量关系。【课时安排】：8学时【教学方法】：采用多媒体教学手段，配合实例教学法，通过对典型例题的讲解启发学生思维，并给与学生适当的课后思考讨论的时间，加深知识掌握的程度。安排一定课时的上机操作。【教学内容】：7.1数据插值与拟合在工程和科学实验中，当研究对象的机理不清楚的时候，经常需要从一组实验观测数据 (xi , yi ) (i = 1, 2, , n)中寻找自变量x与因变量y之间的某种函数关系 y = f(x)。比如测量了人的身高和体重的一些数据，要确定两者的函数关系，但身高与体重的机理我们不清楚，所以寻找尽量吻合这组测量数据的近似函数模型就很重要了。函数 f(x) 的产生办法因观测数据与要求的不同而异，通常可采用数据插值与数据拟合的方法。7.1.1数据插值1.插值问题的描述对给定的一组测量数据，要确定通过所有这些数据点的曲线或曲面的问题就是插值问题。对一维插值问题可以这样描述：设f(x)在区间a, b上连续，x0,x1,xn为 a, b上n+1个互不相同的点，且已知f(x)的一组实验观测数据 (xi , yi ) (i = 1, 2, , n)，要求一个性质优良、便于计算的近似函数j(x)，使得，i =0,1,n (7.1)成立，这就是一维插值问题。其中称a, b为插值区间，点x0,x1,xn为插值节点，函数j(x)为插值函数，f(x)为被插值函数，式(1)为插值条件。求插值函数j(x)的方法称为插值法。关于高维插值可类似定义，本节只介绍一维和二维插值。 2.基本插值方法简介 插值函数的取法很多，可以是代数多项式，也可以是三角多项式或有理函数；可以是a, b上任意光滑函数，也可以是分段光滑函数。对一维插值，最常用最基本的插值方法有：分段多项式插值与三次样条插值；二维插值根据数据分布规律可分为网格节点插值和散乱数据插值，相应的方法有双三次样条插值方法和改进的Shepard方法。具体的方法原理请参阅计算方法的专业书籍，这里不再详细介绍。下面我们着重介绍Matlab中如何实现数据插值。3.插值方法的Matlab实现一维数据插值MATLAB中用函数interp1()来处理一维数据插值，它提供了四种插值方法供选择：线性插值、三次样条插值、三次插值和最临近插值。命令 interp1格式 yi = interp1(x, y, xi, method) %对被插值节点xi, 用method方法进行插值.说明 (1）.输入参数说明：x,y为插值节点，均为向量；xi为任取的被插值点，可以是一个数值，也可以是一个向量；yi为被插值点xi处的插值结果；（2）.其中method是选用的插值方法，具体有： nearest表示最临近插值linear表示线性插值，默认 cubic表示三次插值 spline表示三次样条函数插值注意上述method中所有的插值方法都要求x是单调的，并且xi不能超过x的取值范围, 其中最后一种插值的曲线比较平滑； （3）三次样条插值函数的调用格式有两种等价格式： yi = interp1(x, y, xi, spline) yi = spline(x, y, xi)例1：下表给出了12名成年女子的身高与腿长的测量数据：身高143145146149153155156157158159162164腿长888588929396989796987072试研究身高与腿长的关系，并给出身高为148、150、160时腿长的预测值。解：在matlab中输入代码：x=143145146149153155156157158159162164; %插值节点y=888588929396989796987072; x1=143:0.2:164; %被插值节点，用于确定插值函数.plot(x,y,o);hold on %原始测量数据散点图.y1=interp1(x,y,x1,spline); %求被插值节点处的函数值.yp=interp1(x,y,148 150 160,spline) %求身高为148、150、160时腿长.plot(x1,y1,x,y,r:) %画出插值函数图形及测量数据的折线图.xlabel(身高),ylabel(腿长) %加坐标轴标签输出结果为：yp = 91.4561 92.1942 99.4787 %对应身高的腿长输出图形见图7.1：图7.1注意：1).matlab只会给出被插值节点处的函数值，而不会给出具体的函数解析表达式，这有点类似于我们求微分方程的数值解。需要求点对应的插值（未知的），可以将被插值节点放在xi中；2).图7.1中有三条曲线，其中圆圈点是原始测量数据点（横坐标为插值节点），实线是插值函数图形，虚线是插值节点间的连接折线段。二维数据插值针对二维插值中的插值基点为网格节点和散乱节点，MATLAB中分别提供了函数interp2()和griddata()来进行二维插值。先介绍规则区域上给定数据有规律分布的二维插值。命令 interp2格式 zi = interp2(x, y, z, xi, yi, method) %针对网格节点的二维插值.说明 (1）.输入参数说明：x，y，z为插值节点，其中x和y是自变量，x是m维向量，指明数据网格的横坐标，y是n维向量，指明数据网格的纵坐标，z是阶矩阵，表示相应于网格点的函数值；zi 为被插值点 (xi, yi) 处的插值结果；（2）.其中method是选用的插值方法，具体有： nearest表示最临近插值linear表示双线性插值，默认 cubic表示双三次插值 spline表示双三次样条函数插值注意上述method中所有的插值方法都要求x和y是单调的网格，x和y可以是等距的也可以是不等距的。xi和yi应是方向不同的向量，即一个是行向量，另一个是列向量。几种方法中最后一种插值的曲面比较平滑。例2：已知在某山区测得一些地点的高程如下表。其平面区域为，试用不同的插值方法作出该山区的地貌图。 x y12001600200024002800320036004000120011301250128012307409005007001600132014501420140013007009008502000139015001500140090017076095024001500120017013501450120011507728001500120017015501600155013807703200150015501600155016001600160015503600148015001550157143013001200980解：输入程序代码：x=1200:400:4000;y=1200:400:3600;z=1130 12501280 1230740 900500 700132014501420 14001300 700900 850139015001500 1400900 170760 95015001200170 13501450 12001150 7715001200170 15501600 15501380 770150015501600 15501600 16001600 1550148015001550 1571430 13001200 980;%原始数据的山区地貌图figure(1)meshz(x,y,z)xlabel(X),ylabel(Y),zlabel(Z)title(原始数据地貌图)%为平滑曲面，加密网格x1=1200:50:4000;y1=1200:50:3600;%最临近插值figure(2)zn=interp2(x,y,z,x1,y1,nearest);surfc(x1,y1,zn)xlabel(X),ylabel(Y),zlabel(Z)title(最临近插值地貌图)%双线性插值figure(3)zl=interp2(x,y,z,x1,y1,linear);surfc(x1,y1,zl)xlabel(X),ylabel(Y),zlabel(Z)title(双线性插值地貌图)%双三次插值figure(4)zc=interp2(x,y,z,x1,y1,cubic);surfc(x1,y1,zc)xlabel(X),ylabel(Y),zlabel(Z)title(双三次插值地貌图)%双三次样条函数插值figure(5)zs=interp2(x,y,z,x1,y1,spline);surfc(x1,y1,zs)xlabel(X),ylabel(Y),zlabel(Z)title(双三次样条函数插值地貌图)输出可视化图形分别见图7.2图7.6： 图7.2 图7.3 图7.4 图7.5图7.6从图形可以看出，原始数据地貌图是很粗糙的，因为测量点比较少。几种插值方法中最临近插值和双线性插值效果较差，而最后一种插值的曲面比较平滑，效果较好。如果给定的数据是在规则区域上的散乱数据或随机分布的数据，即数据不是在网格上取的，则可用函数griddata()来解决二维插值问题。命令 griddata格式 zi =griddata(x, y, z, xi, yi, method) %针对散乱数据的二维插值.说明 （1）.输入参数说明：x，y，z都是n维向量，分别指明所给插值节点的横坐标、纵坐标和z坐标；zi 为被插值点 (xi, yi) 处的插值结果；xi和yi应是方向不同的向量，即一个是行向量，另一个是列向量；（2）.其中method是选用的插值方法，具体有： nearest表示最临近插值linear表示双线性插值，默认 cubic表示双三次插值 v4表示matlab提供的插值方法其中v4方法比较好。针对二维散乱插值问题，在matlab中还提供了两个插值函数：e01sef()和e01sff()。通常两者要配合使用，其调用格式为： fnodes, a, rnw, b, c = e01sef(x, y, z) sz(i, j), ifail =e01sff(x, y, z, rnw, fnodes, sx(i),sy(j)其中： x, y, z 为插值节点, 均为n维向量； sx(i), sy(j) 为被插值节点； sz(i, j)为被插值点(sx(i), sy(j)处的插值结果；其他输出参数涉及插值算法。两个函数中e01sef输出fnodes和rnw为确定插值的参数，它们是e01sff 需要的输入参数，因此两函数需配合使用。例3：在某海域测得一些点（x, y）处的水深 z（单位：英尺），见下表，水深数据是在低潮时测得的。船的吃水深度为 5 英尺，问在矩形区域(75, 200)(-50, 150)内的哪些地方船要避免进入（AMCM 86A题）。x129.0140.073.588.0185.5195.075.5y7.5141.523.0147.022.5137.585.5z4868688x157.577.577.081.0162.0162.0117.5y-6.5-81.03.056.5-66.584.0-33.5z99889491. 基本假设 除了一些散乱的测量数据外，题目没有给出其他信息。为了简化问题，首先给出以下合理假设：(1).所给测量数据是精确可用的；(2).该海域海底是平滑的, 不存在珊瑚礁、水底峡谷、山脊等突变地形。2.问题分析在假设基础上，可以考虑用某种光滑的曲面去拟合逼近已知的数据点或以已知的数据点为基础，利用二维插值方法补充一些点的水深，然后作出海底曲面图和等高线图，并求出水深小于5的海域范围。3.问题求解(1).先作出测量点的分布散点图：输入代码：x=129.0 140.0 73.5 88.0 185.5195.0 75.5 157.5 77.577.0 81.0 162.0 162.0 117.5;y=7.5 141.5 23.0 147.0 22.5 137.5 85.5 -6.5 -81.0 3.0 56.5 -66.5 84.0 -33.5;z=-4 8 6 8 6 8 8 9 9 88 9 4 9; %相当于以海平面作为xy平面.plot(x,y,o)title(测量数据xy平面分布图)图7.7从图7.7可以直观看出在矩形区域(75, 200)(-50, 150)内测量点是散乱分布的，所以用matlab中的griddata函数作海底二维插值。(2).海底地貌图绘制：xi=75:1:200;yi=-50:1:150;zi=griddata(x,y,z,xi,yi,v4); %用v4方法作散乱数据的二维插值.figure(2)mesh(xi,yi,zi)xlabel(X);ylabel(Y);zlabel(Z);title(海底地貌图)rotate3d %产生可旋转的3D图形.图7.8图7.8给出了在矩形海域(75, 200)(-50, 150)内，海底的地貌图。其中明显有些区域海水深度较浅，不适合船只通过，为了更准确确定吃水深度为5英尺的船只不能进入的海域，下面我们绘出该区域的危险区域的海底地貌图。图7.9(3).危险区域，即水深小于5英尺的海底地貌图：u,v=find(zi-5); %在zi矩阵中找出值大于-5的元素并将坐标存入向量u,v,也就是找出水深小于5英尺的插值点.zzi=zeros(size(zi)-5; %产生一个与zi同型的, 元素均为-5的矩阵.for i=1:length(u) zzi(u(i),v(i)=zi(u(i),v(i); %将zi中值大于-5的元素替换zzi中对应元素.endfigure(3)mesh(xi,yi,zzi) %绘出以z坐标等于-5为基底面的危险区域地貌图.rotate3d xlabel(X);ylabel(Y);zlabel(Z);title(海底危险区域地貌图)图7.9给出了水深小于5英尺的海底区域的地貌图，它可以理解成用z=-5的平面去截图7.8所得上半部分的图形。(4).危险区域平面图：危险区域平面图，即绘制z=5的等高线，代码和等值线图如下：.figure(4)contour(xi,yi,zi,-5,-5,r) %作深度为5的海底等值线图xlabel(X);ylabel(Y)title(危险区域平面图)图7.10图7.10给出了平面矩形区域(75, 200)(-50, 150)内的危险区域，其中线条内部区域是船只要避免进入的区域。 7.1.2 数据拟合1.拟合问题的描述对给定的一组测量数据(xi , yi ) (i = 1, 2, , n)，如果要寻找变量x与y的函数关系的近似表达式，前面提到的插值方法在一定程度上可以较好解决问题，但它有明显的缺陷。一是测量数据常带有测试误差，而插值多项式又经过所有这些点，保留了误差；二是如果实验测量数据较多，则必然出现高次插值多项式，这样近似效果并不理想。而数据拟合却能较好的避免这些问题。所谓数据拟合就是从给定的一组数据出发，寻找函数的一个近似表达式y=j(x)，要求该函数在某种准则下能尽量反应数据的整体变化趋势，而不一定经过所有数据点(xi , yi )。数据拟合问题也叫曲线拟合问题，其中y=j(x)称为拟合曲线。求解曲线拟合问题要求我们先要给出拟合的函数类型，然后利用测量数据按照一定的方法求出参数，其中最常用的解法是最小二乘法。2.曲线拟合的最小二乘原理对给定的一组测量数据(xi , yi ) (i = 1, 2, , n)，拟合的函数为，其中为待定系数。为使函数在整体上尽可能与给定数据点接近，我们常采用n个已知点(,)与曲线的距离(偏差)的平方和最小，即来保证每个偏差的绝对值都很小，这一原则称为最小二乘原则，根据最小二乘原则确定拟合函数的方法称为最小二乘法，满足上述要求的参数取值称为该问题的最小二乘解。确定最小二乘解的问题涉及多元函数的极值问题，在后面回归分析中我们还会有所接触，这里我们不再详细介绍，有兴趣的同学可以参考数学分析或计算方法类教科书。我们的重点是掌握应用MATLAB工具进行最小二乘估计，即进行曲线拟合。事实上在MATLAB中已有现成的求最小二乘问题的函数polyfit，称为多项式拟合函数，并且这个函数允许多项式的次数可以是任意次的。同学们即使没有这方面的数学基础，只要软件工具熟悉了，一样可以做的很好。曲线拟合分为线性最小二乘拟合和非线性最小二乘拟合。如果拟合函数的待定系数a0, a1, a2, , am 全部以线性形式出现，我们称之为线性最小二乘拟合；如果拟合函数的待定参数a1, a2, , am不能全部以线性形式出现，如指数拟合函数等，这就非线性最小二乘拟合问题。注意：最小二乘原理中的偏差也可以是点到直线的垂直线段的长度（点到直线的距离），也可以是点沿（平行）X轴方向到直线的距离（横向距离）或点沿（平行）Y轴方向到直线的距离（纵向距离）。其中最常用是纵向距离。3.拟合函数的确定最小二乘法中，确定拟合函数类型是很关键的。常用的有两种方式：1).通过机理分析建立数学模型来确定j (x) ，比如前面提到的人口增长的logistic模型就是机理分析法推导出来的，但参数的确定需要用到统计数据进行曲线拟合；2).如果无现成的规则或事物机理不清楚，可以通过散点图，结合曲线的形状变化趋势进行分析，建立经验模型。比如图7.11中数据点基本分布在一条带型区域上，所以可以考虑用直线模型作为经验模型；对图7.12，可以看作二次曲线，也可以用反指数函数作为经验模型。图7.11 图7.12当然这具有一定的主观性，因此要求我们多观察分析，找出最适合的拟合函数。一种好的处理方法是对同一问题，分别选择不同的函数进行最小二乘拟合，比较各自误差的大小，从中选出误差较小的作为拟合函数。4.用Matlab进行曲线拟合值线性最小二乘拟合MATLAB中线性最小二乘拟合其实就是作多项式拟合，他可以看作是函数的线性组合。命令 polyfit，polyval格式 a=polyfit(x, y, m) %对给定数据作m次多项式拟合.y=polyval(a, x) %调用拟合出来的多项式计算在x处的值，即求预测值.参数说明：x,y同长度的数组，需要拟合的实验数据；a输出拟合多项式系数a = a1, a2, , am+1 (数组)，从高次到低次；m拟合多项式次数。例4：对本节例题1中数据作多项式拟合。解：输入代码：x=143145146149153155156157158159162164;y=888588929396989796987072;plot(x,y,o); %绘制散点图.a=polyfit(x,y,1); %作1次多项式拟合.z=polyval(a,x) %求预测值.plot(x,y,o,x,z,r-) %拟合效果对比图.xlabel(身高),ylabel(腿长)拟合效果对比图如下：图7.13拟合系数为：a = 0.7465 -20.4874，即a1=0.7465，a2=-20.4874。例5：已知一室模型快速静脉注射下的血药浓度数据 (t=0注射300mg)如下表：t(h)0.250.511.523468c(mg/ml)19.2118.1515.3614.712.899.327.455.243.01试求血药浓度随时间的变化规律c(t)。解：1).作出散点图，观察规律，代码和图形如下：t=0.25 0.5 1 1.5 2 3 4 6 8;c=19.21 18.15 15.36 14.7 12.89 9.32 7.45 5.24 3.01；figure(1)plot(t,c,+); %绘制散点图xlabel(时间),ylabel(药物浓度)图7.14图7.14显示血药浓度随时间的变化规律是时间越长，浓度越低，而且开始时下降速度快，逐渐减弱，两者关系大致呈现类抛物线规律，因此考虑用多项式作拟合。2).作线性拟合：a=polyfit(t,c,2) %2次多项式拟合.z=polyval(a,t) figure(2)plot(t,c,+,t,z,r-) %2次多项式拟合效果对比图.图7.15拟合多项式系数为：a =0.2564 -4.788 19.8437。故所求的二次拟合多项式为。从图7.15可以看出模型拟合效果良好，除个别点存在较大偏差外，基本反映了实验数据的变化趋势。为了寻求更好的拟合函数，下面考虑用3次多项式进行拟合，对比效果见图7.16，程序略。图7.16非线性最小二乘拟合MATLAB中提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。它们所采用的算法是一样的。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)，但两者定义f(x)的方式是不同的。这里主要介绍lsqcurvefit的应用。命令 lsqcurvefit格式 x = lsqcurvefit(fun,x0,xdata,ydata) %用xdata, ydata拟合fun中的参数并返回x.x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub) x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)x,resnorm = lsqcurvefit() x,resnorm,residual = lsqcurvefit() x,resnorm,residual,exitflag = lsqcurvefit() x,resnorm,residual,exitflag,output = lsqcurvefit() x,resnorm,residual,exitflag,output,lambda = lsqcurvefit() x,resnorm,residual,exitflag,output,lambda,jacobian =lsqcurvefit()参数说明(1).,是满足关系ydata=F(x, xdata)的已知数据点；为参数向量；x0为迭代初值；lb、ub 为解向量的下界和上界，若没有指定界，则lb= ，ub= ；options是选项；(2). fun是一个事先建立的定义拟合函数的M-文件, 自变量为x和；(3). resnorm=sum (fun(x,xdata)-ydata).2)，即在x 处残差的平方和； residual=fun(x,xdata)-ydata，即在x 处的残差； exitflag为终止迭代的条件； output为输出的优化信息； lambda为解x 处的Lagrange乘子； jacobian为解x 处拟合函数fun 的jacobian 矩阵。例6：已知上个世纪近两百年的美国人口统计数据如下表，试建立描述人口增长变化的数学模型，并用所给数据拟合出相应参数。年份1790180018718201830184018501860人口(百万)3.95.37.29.612.917.123.231.4年份1870188018901900197192019301940人口(百万)38.650.262.976.092.076.5123.2131.7年份195019601970198019902000人口(百万)150.7179.3204226.5251.4281.4解：1).数学模型描述:在第3章我们已经介绍了人口增长的经典数学模型，logistic模型：其中参数x0为人口初始数据，xm为最大人口容量，r为固定增长率，x(t)为第t年人口数量。只要估计出模型的参数，就可以进行预测。其中x0可以用第一年的统计数据代替，另两个参数需要用曲线拟合的方法进行估计。2).对参数作非线性拟合(1).先定义外部函数，文件名为usa_renkou.mfunction f=usa_renkou(x,year)f=x(1)./(1+(x(1)/3.9-1).*exp(-x(2).*year); %x(1)=xm;x(2)=r;x0=3.9(2).输入主程序year=0:1:21; %将年份用0-21替换，即0表示1790年，21对应2000年.usp=3.900 5.300 7.200 9.600 12.90 17.7 23.20 31.40 38.60 50.20 62.90 76 92 76.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 281.4;x0=70 1.8;x=lsqcurvefit(usa_renkou,x0,year,usp) z=usa_renkou(x,year); %调用函数进行人口预测.plot(year,usp,o,year,z,r)xlabel(年份);ylabel(人口数量);title(1790-2000年美国人口数量拟合效果图)(3).运行结果为：x = 342.4339 0.2735z = Columns 1 through 16 3.9000 5.786 6.6844 8.7337 1.3899 14.8180 19.2181 24.8257 31.9081 40.7517 51.6387 64.8095 80.476 98.4357 118.6714 140.6687 Columns 17 through 22 163.7595 187.1252 209.9082 231.3333 250.8067 267.9655即所求参数为xm=342.4339，r=0.2735，z为预测值向量。对比效果见图7.17。图7.17(4).预测与误差：这里我们仅用拟合出来的函数计算2005年的美国人口数据，只需执行usa_renkou(x,26)后便可得2005年预测人口为319.7976百万。而实际上2005年美国人口数为296百万，预测误差为7.4%。从图7.17也可以看出模型在19世纪中叶到20世纪中叶拟合效果不太理想，但这段时间前后拟合效果还是很好的。例7：NPK施肥问题（1992年全国大学生数学建模竞赛A题） 某地区作物生长所需的营养素主要是氮(N)、钾(K)、磷(P)。某作物研究所在某地区对土豆与生菜做了一定数量的实验，其中土豆的实验数据如下列表所示，其中ha表示公顷，t表示吨，kg表示公斤。当一个营养素的施肥量变化时，总将另两个营养素的施肥量保持在第七个水平上，如对土豆产量关于N的施肥量做实验时，P与K的施肥量分别取为196kgha与372kgha。试分析施肥量与产量之间关系，并对所得结果从应用价值与如何改进等方面做出估价。土 豆NPK施肥量(kg/ha)产量(t/ha)施肥量(kg/ha)产量(t/ha)施肥量(kg/ha)产量(t/ha)015.18033.46018.983421.362432.474727.356725.724936.069334.867132.297337.9614039.5213534.039841.0418638.4420239.4514740.0927937.7325943.1519641.2637238.4333643.4624542.1746543.8740440.8329440.3655842.7747130.7534242.7365146.22生 菜NPK施肥量(kg/ha)产量(t/ha)施肥量(kg/ha)产量(t/ha)施肥量(kg/ha)产量(t/ha)01102063901575281270499484716765614569812469316898416271471433140162411217751961771861756168225929421942791920224216339122643721797280193448921344651584336161258722075582011392141168524536511940解: 为考察氮、磷、钾三种肥料对作物的施肥效果，我们以氮、磷、钾的施肥量为自变量，土豆产量为因变量描点作图，先观察数据的分布特点和变化规律。1).输入数据，绘制施肥量与土豆产量的散点图：nshi=0 34 67 71 135 202 259 336 404 471; %N的施肥量.nchan=15.18 21.36 25.72 32.29 34.03 39.45 43.15 43.46 40.83 30.75;pshi=0 24 49 73 98 147 196 245 294 342; %P的施肥量.pchan=33.46 32.47 36.06 37.96 41.04 40.09 41.26 42.17 40.36 42.73;kshi=0 47 93 140 186 279 372 465 558 651; %K的施肥量.kchan=18.98 27.35 34.86 39.52 38.44 37.73 38.43 43.87 42.77 46.22;figure(1)plot(nshi,nchan,*) %N施肥量与产量关系图.figure(2)plot(pshi,pchan,ro) %P施肥量与产量关系图.figure(3)plot(kshi,kchan,b+) %K施肥量与产量关系图.2).运行后，输出图7.18图7.20： 图7.18 N施肥量与产量关系图 图7.19 P