第四讲数据分析方法

第四讲数据分析方法第一节、数据拟合问题：给定一批数据点（输入变量与输出变量的数据），需确定满足特定要求的曲线或曲面。如果输入变量和输出变量都只有一个，则属于一元函数的拟合和插值；而若输入变量有多个，则为多元函数的拟合和插值（有点回归分析的意思）解决方案：（1）若要求所求曲线（面）通过所给所有数据点，就是插值问题；（2）若不要求曲线（面）通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合。注意：插值和拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。而面对一个实际问题，究竟应该用插值还是拟合，有时容易确定，有时则并不明显。例1：下面数据是某次实验所得，希望得到X和f之间的关系？xf12479121315171.53.96.611.715.618.819.620.621.1曲线拟合问题最常用的解法——最小二乘法的基本思路第一步：确定拟合的函数类型yf(x;a1,a2,,am)，其中a1,a2,,am为待定系数。（函数类型的确定可以根据内在的规律确定，如果无现成的规则，则可以通过散点图，联系曲线的形状进行分析）第二步：确定a1,a2,,am的最小二乘准则：要求n个已知点(xi,yi)与曲线yf(x)n的距离di的平方和(yifx2最小。())ii1用MATLAB作拟合1．多项式拟合。作多项式ya0xma1xm1am拟合，可利用a=polyfit(x,y,m)—其中x,y为给出的数据，m为多项式的次数。多项式在x处的值y可用以下命令计算：y=polyval（a,x）2．用MATLAB作非线性最小二乘拟合Matlab的提供了两个求非线性最小二乘拟合的函数：lsqcurvefit和lsqnonlin。两个命令都要先建立M-文件fun.m，在其中定义函数f(x)。（1）x=lsqcurvefit(‘fun’,x0,xdata,ydata);（2）x=lsqcurvefit(‘fun’,x0,xdata,ydata,options);（3）x=lsqcurvefit(‘fun’,x0,xdata,ydata,options,’grad’);（4）[x,options]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);（5）[x,options,funval]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);1（6）[x,options,funval,Jacob]=lsqcurvefit(‘fun’,x0,xdata,ydata,…);（7）x=lsqnonlin（‘fun’，x0）；（8）x=lsqnonlin（‘fun’，x0，options）；（9）x=lsqnonlin（‘fun’，x0，options，‘grad’）；（10）[x，options]=lsqnonlin（‘fun’，x0，…）；（11）[x，options，funval]=lsqnonlin（‘fun’，x0，…）；在使用lsqcurvefit与lsqnonlin命令时，共同的问题是要先知道函数的类型，而拟合其实是决定函数中的待定系数。第二节插值插值的基本问题：给出n个数对，(Pi,f(Pi)),i1,2,,n，求点P处对应的函数值f(P)。一、一维插值已知n1个节点(xi,yi),i0,1,,n，求任意点*处的函数值*。常用的插值方法xy有拉格朗日多项式插值、牛顿插值、分段线性插值、Hermite插值和三次样条插值。分段线性插值：将各数据点用折线连接起来多项式插值：求一个多项式通过所有数据点，可以假设出多项式的系数，最后通过求解方程得到每个系数（拉格朗日插值，用n次多项式描述n1个点）样条插值：分段多项式的光滑连接（三次样条插值）牛顿插值：利用节点之间的各阶差商和差分构造多项式Hermite插值：对插值函数，不仅要求它在节点处与函数同值，而且要求它与函数有相同的一阶、二阶甚至更高阶的导数值（1）MATLAB命令：y=interp1(x0,y0,x,'method')method指定插值的方法，默认为线性插值。其值可为：'nearest'最近项插值'linear''spline''cubic'线性插值立方样条插值立方插值。所有的插值方法要求x0是单调的。当x0为等距时可以用快速插值法，使用快速插值法的格式为'*nearest'、'*linear'、'*spline'、'*cubic'。（2）三次样条插值在Matlab中的实现在Matlab中数据点称之为断点。如果三次样条插值没有边界条件，最常用的方法，就是采用非扭结（not-a-knot）条件。这个条件强迫第1个和第2个三次多项式的三阶导数相等。对最后一个和倒数第2个三次多项式也做同样地处理。Matlab中三次样条插值也有现成的函数：y=interp1(x0,y0,x,'spline')；y=spline(x0,y0,x)；pp=csape(x0,y0,conds),pp=csape(x0,y0,conds,valconds)，y=ppval(pp,x)。其中x0,y0是已知数据点，x是插值点，y是插值点的函数值。对于三次样条插值，我们提倡使用函数csape，csape的返回值是pp形式，要求插值点2的函数值，必须调用函数ppval。pp=csape(x0,y0)：使用默认的边界条件，即Lagrange边界条件。pp=csape(x0,y0,conds,valconds)中的conds指定插值的边界条件，其值可为：'complete'边界为一阶导数，一阶导数的值在valconds参数中给出，若忽略valconds参数，则按缺省情况处理。'not-a-knot'非扭结条件'periodic''second'周期条件边界为二阶导数，二阶导数的值在valconds参数中给出，若忽略valconds参数，二阶导数的缺省值为[0,0]。'variational'设置边界的二阶导数值为[0,0]。对于一些特殊的边界条件，可以通过conds的一个12矩阵来表示，conds元素的取值为0，1，2。conds(i)=j的含义是给定端点i的j阶导数，即conds的第一个元素表示左边界的条件，第二个元素表示右边界的条件，conds=[2,1]表示左边界是二阶导数，右边界是一阶导数，对应的值由valconds给出。例2机床加工待加工零件的外形根据工艺要求由一组数据(x,y)给出（在平面情况下），用程控铣床加工时每一刀只能沿x方向和y方向走非常小的一步，这就需要从已知数据得到加工所要求的步长很小的(x,y)坐标。表中给出的x,y数据位于机翼断面的下轮廓线上，假设需要得到x坐标每改变0.1时的y坐标。试完成加工所需数据，画出曲线，并求出x0处的曲线斜率和13x15范围内y的最小值。xy03579112.0121.8131.2141.0151.601.21.72.02.1要求用分段线性和三次样条计算。解编写以下程序：x0=[035791112131415];y0=[01.21.72.02.12.01.81.21.01.6];x=0:0.1:15;y2=interp1(x0,y0,x);y3=interp1(x0,y0,x,'spline');pp1=csape(x0,y0);y4=ppval(pp1,x);pp2=csape(x0,y0,'second');y5=ppval(pp2,x);[x',y1',y2',y3',y4',y5']subplot(2,2,1)plot(x0,y0,'+',x,y2)title('Piecewiselinear')subplot(2,2,2)plot(x0,y0,'+',x,y3)3title('Spline1')subplot(2,2,3)plot(x0,y0,'+',x,y4)title('Spline2')dx=diff(x);dy=diff(y3);dy_dx=dy./dx;dy_dx0=dy_dx(1)ytemp=y3(131:151);ymin=min(ytemp);index=find(y3==ymin);xmin=x(index);[xmin,ymin]二、二维插值：前面讲述的都是一维插值，即节点为一维变量，插值函数是一元函数（曲线）。若节点是二维的，插值函数就是二元函数，即曲面。如在某区域测量了若干点（节点）的高程（节点值），为了画出较精确的等高线图，就要先插入更多的点（插值点），计算这些点的高程（插值）。（1）（2）网格节点：已知mn个节点(xi,yj,zij),i1,2,,m;j1,2,,n，构造一个二元函数zf(x,y)通过全部已知节点，即zijf(xi,yj)，再利用zf(x,y)插值。散乱节点：已知n个节点(xi,yi,zi),i1,2,,n，构造一个二元函数zf(x,y)通过全部已知节点，即zif(xi,yi)，再利用zf(x,y)插值。二维插值方法：（1）最邻近插值：二维或高维情形的最邻近插值，与被插值点最邻近的节点的函数值即为所求。（2）分片线性插值：将四个插值点（矩形的四个顶点）处的函数值依次简记为：f(xi,yj)=f1，f(xi+1,yj)=f2，f(xi+1,yj+1)=f3，f(xi,yj+1)=f4分两片的函数表达式如下：第一片(下三角形区域)：f(x,y)f1(f2f1)(xxi)(f3f2)(yyj)第二片(上三角形区域)：f(x,y)f1(f4f1)(yyj)(f3f4)(xxi)（3）双线性插值：双线性插值是一片一片的空间二次曲面构成。双线性插值函数的形式如下：f(x,y)(axb)(cyd)。其中有四个待定系数，利用该函数在矩形的四个顶点（插值节点）的函数值，得到四个代数方程，正好确定四个系数。用MATLAB作网格节点数据的插值：4z=interp2(x0,y0,z0,x,y,’method’)—method可以选择‘nearest’：最邻近插值；‘linear’：双线性插值；‘cubic’：双三次插值；缺省时,双线性插值。如果是三次样条插值，可以使用命令pp=csape({x0,y0},z0,conds,valconds)，z=fnval(pp,{x,y})其中x0,y0分别为m维和n维向量，z0为mn维矩阵，z为矩阵，它的行数为x的维数，列数为y的维数，表示得到的插值，具体使用方法同一维插值。例3在一丘陵地带测量高程，x和y方向每隔100米测一个点，得高程如下表，试插值一曲面，确定合适的模型，并由此找出最高点和该点的高程。xy100200300400100200697712674626300624630598552400500478450478420412400334310636698680662解编写程序如下：x=100:100:500;y=100:100:400;z=[636697624478450;698712630478420;680674598412400;662626552334310];pp=csape({x,y},z')xi=100:10:500;yi=100:10:400cz1=fnval(pp,{xi,yi})cz2=interp2(x,y,z,xi,yi','spline')[i,j]=find(cz1==max(max(cz1)))x=xi(i),y=yi(j),zmax=cz1(i,j)用MATLAB作散点数据的插值计算cz=griddata（x,y,z,cx,cy,‘method’）—method可以选择‘nearest’：最邻近插值；‘linear’：双线性插值；‘cubic’：双三次插值；'v4'：Matlab提供的插值方法；缺省时,双线性插值例4在某海域测得一些点（x,y）处的水深z由下表给出，在矩形区域（75,200）×(-50,150)内画出海底曲面的图形。x129y7.5z4140141.58103.5236881478185.522.56195137.58105157.