用MATLAB进行数据拟合

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,*,Click to edit Master title style,Click to edit Master text styles,Second level,Third level,Fourth level,Fifth level,*,用MATLAB进行数据拟合,拟 合,1.,拟合问题引例,拟 合 问 题 引 例 1,温度,t(,0,C) 20.5 32.7 51.0 73.0 95.7,电阻,R(,) 765 826 873 942 1032,已知热敏电阻数据：,求60,0,C,时的电阻R。,设,R=at+b,a,b,为待定系数,拟 合 问 题 引 例 2,t (h) 0.25 0.5 1 1.5 2 3 4 6 8,c (,g/ml),19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01,已知一室模型快速静脉注射下的血药浓度数据(t=0,注射300mg),求血药浓度随时间的变化规律,c(t).,作半对数坐标系,(semilogy),下的图形,MATLAB(aa1),曲 线 拟 合 问 题 的 提 法,一组二维数据，即平面上 n个点xi,yi) i=1,n, 寻求一个函数曲线y=f(x), 使 f(x) 在某种准那么下与所有数据点最为接近，即曲线拟合得最好。,+,+,+,+,+,+,+,+,+,x,y,y=f(x),(x,i,y,i,),i,i 为点xi,yi) 与曲线 y=f(x) 的距离,拟合与插值的关系,函数插值与曲线拟合都是要根据一组数据构造一个函数作为近似，由于近似的要求不同，二者的数学方法上是完全不同的。,实例：,下面数据是某次实验所得，希望得到X和 f之间的关系？,MATLAB(cn),问题：,给定一批数据点，需确定满足特定要求的曲线或曲面,解决方案：,假设不要求曲线面通过所有数据点，而是要求它反映对象整体的变化趋势，这就是数据拟合，又称曲线拟合或曲面拟合。,假设要求所求曲线面通过所给所有数据点，就是插值问题；,最临近插值、线性插值、样条插值与曲线拟合结果：,曲线拟合问题最常用的解法线性最小二乘法的根本思路,第一步:先选定一组函数 r1(x), r2(x), rm(x), mn, 令,f(x)=a1r1(x)+a2r2(x)+ +amrm(x) 1,其中 a1,a2, am 为待定系数。,第二步: 确定a1,a2, am 的准那么最小二乘准那么：,使n个点xi,yi) 与曲线 y=f(x) 的距离i 的平方和最小 。,记,问题归结为，求,a,1,a,2, a,m,使,J(a,1,a,2, a,m,),最小。,线性最小二乘法的求解：预备知识,超定方程组,：方程个数大于未知量个数的方程组,即 Ra=y,其中,超定方程一般是不存在解的矛盾方程组。,如果有向量a使得 达到最小，,则称a为上述,超定方程的最小二乘解,。,线性最小二乘法的求解,定理：当RTR可逆时，超定方程组3存在最小二乘解，且即为方程组,RTRa=RTy,的解：a=(RTR)-1RTy,所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。,其中,Ra=y （3）,线性最小二乘拟合 f(x)=a,1,r,1,(x)+ +a,m,r,m,(x),中函数,r,1,(x),r,m,(x),的选取,1. 通过机理分析建立数学模型来确定 f(x)；,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,f=a,1,+a,2,x,f=a,1,+a,2,x+a,3,x,2,f=a,1,+a,2,x+a,3,x,2,f=a,1,+a,2,/x,f=ae,bx,f=ae,-bx,2. 将数据,(x,i,y,i,) i=1, n,作图，通过直观判断确定,f(x)：,用MATLAB解拟合问题,1、线性最小二乘拟合,2、非线性最小二乘拟合,用,MATLAB,作线性最小二乘拟合,1. 作多项式,f(x)=a,1,x,m,+ +a,m,x+a,m+1,拟合,可利用已有程序,:,a=polyfit(x,y,m),2. 对超定方程组,可得最小二乘意义下的解。,，用,3.多项式在x处的值y可用以下命令计算：,y=polyvala，x,输出拟合多项式系数,a=,a,1, a,m,a,m+1,(,数组)）,输入同长度,的数组X，Y,拟合多项,式次数,即要求 出二次多项式:,中 的,使得:,例1 对下面一组数据作二次多项式拟合,1输入以下命令：,x=0:0.1:1;,y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;,R=(x.2) x ones(11,1)；,A=Ry,MATLAB(zxec1),解法1,用解超定方程的方法,2计算结果,1输入以下命令：,x=0:0.1:1;,y=-0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2;,A=polyfit(x,y,2),z=polyval(A,x);,plot(x,y,k+,x,z,r) %作出数据点和拟合曲线的图形,解法,2,用多项式拟合的命令,MATLAB(zxec2),例2 观测数据点如表所示,x,y,0,-0.447,0.1,1.978,0.2,3.28,0.3,6.16,0.4,7.08,0.5,7.34,0.6,7.66,0.7,9.56,0.8,9.48,0.9,9.3,1,11.2,分别用3次和6次多项式曲线拟合这些数据点.,x=0:0.1:1,y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2,plot(x,y,k.,markersize,25),axis(0 1.3 -2 16),p3=polyfit(x,y,3),p6=polyfit(x,y,6),编写Matlab程序如下:,s=polyval(p3,t),s1=polyval(p6,t),hold on,plot(t,s,r-,linewidth,2),plot(t,s,b-,linewidth,2),grid,x=0:0.1:1,y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2,plot(x,y,k.,markersize,25),axis(0 1.3 -2 16),p3=polyfit(x,y,3),p6=polyfit(x,y,6),例3 用切削机床进展金属品加工时, 为了适当地调整机床, 需要测定刀具的磨损速度. 在一定的时间测量刀具的厚度, 得数据如表所示:,切削时间 t/h,0,30.0,1,29.1,2,28.4,3,28.1,4,28.0,5,27.7,6,27.5,7,27.2,8,27.0,刀具厚度 y/cm,切削时间 t/h,9,26.8,10,26.5,11,26.3,12,26.1,13,25.7,14,25.3,15,24.8,16,24.0,刀具厚度 y/cm,解:,描出散点图, 在命令窗口输入:,t=0:1:16,y=30.0 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0,plot(t,y,*),解:,描出散点图, 在命令窗口输入:,t=0:1:16,y=30.0 29.1 28.4 28.1 28.0 27.7 27.5 27.2 27.0 26.8 26.5 26.3 26.1 25.7 25.3 24.8 24.0,plot(t,y,*),a =,hold on,plot(t, y1), hold off,a=polyfit(t,y,1),例3 用切削机床进展金属品加工时, 为了适当地调整机床, 需要测定刀具的磨损速度. 在一定的时间测量刀具的厚度, 得数据如表所示:,切削时间 t/h,0,1,2,3,4,5,6,7,8,刀具厚度 y/cm,切削时间 t/h,9,10,11,12,13,14,15,16,刀具厚度 y/cm,拟合曲线为:,yt,例4,2. 93,3. 03,应力-应变关系可以用一条指数曲线来描述, 即假设,式中, 表示应力, 单位是 N/m,2,; 表示应变.,应力-应变关系可以用一条指数曲线来描述, 即假设,式中, 表示应力, 单位是 N/m,2,; 表示应变.,解,选取,指数函数,作拟合时, 在拟合前需作变量代换,化为,k,1,k,2,的线性函数.,于是,令,即,在命令窗口输入:,x=500*1.0e-6 1000*1.0e-6 1500*1.0e-6 2000*1.0e-6 2375*1.0e-6,y=3.103*1.0e+3 2.465*1.0e+3 1.953*1.0e+3 1.517*1.0e+3 1.219*1.0e+3,z=log(y),a=polyfit(x,z,1),k1=exp(8.3009),w=1.55 2.47 2.93 3.03 2.89,plot(x,w,*),y1=exp(8.3009)*x.*exp( -494.5209*x),plot(x,w,*,x,y1,r-),应力-应变关系可以用一条指数曲线来描述, 即假设,式中, 表示应力, 单位是 N/m,2,; 表示应变.,拟合曲线为:,令,那么,求得,于是,在实际应用中常见的拟合曲线有:,直线,多项式,一般,n,=2, 3, 不宜过高.,双曲线(一支),指数曲线,二、人口预测线性模型,对于开场提出的实验问题, 代如数据，计算得,从而得到人口数与年份的函数关系为,把x=1999代如，估算出1999年的人口数为,y,1999年实际人口数量为.亿。,线性预测模型,英国人口学家Malthus根据百余年的人口统计资料，于1798年提出了著名的人口自然增长的,指数增长模型,。,三、人口预测的Malthus模型,根本假设 : 人口(相对)增长率 r 是常数,设,x,(,t,) 为,t,时刻的人口数量,t,=0时人口数为,x,0,指数增长模型,实际中，常用,1. 由前100年的数据求出美国的人口增长,Malthus,模型,。,2. 预测后100年每隔10年的人口状况。,3. 根据预测的人口状况和实际的人口数量,讨论人口模型的改进情况。,美国1790年1980年每隔10年的人口记录,226.5,204.0,179.3,150.7,131.7,123.2,106.5,92.0,76.0,62.9,人口(百万),1980,1970,1960,1950,1940,1930,1920,1910,1900,1890,年份,50.2,38.6,31.4,23.2,17.1,12.9,9.6,7.2,5.3,3.9,人口(百万),1880,1870,1860,1850,1840,1830,1820,1810,1800,1790,年份,例,解：,取得最小值.,其中,表示人口数量,。,表示年份,解方程组,:,即得参数,的值,.,使得,问题转化为求参数, prog41.m %,% This program is to predict the number of population %,format long,t1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;,t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;,x1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;,x2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;,lnx1=log(x1); lnx2=log(x2);,a12=sum(t1);a11=10;a21=a12;a22=sum(t1.2); d1=sum(lnx1);d2=sum(lnx1.*t1);A=a11,a12;a21,a22;D=d1;d2;ab=inv(A)*D;,disp(a=);disp(ab(1);disp(b=);disp(ab(2);,for,i=1:10 xx1(i)=exp(ab(1)+ab(2)*t1(i);,end,for,i=1:10 xx2(i)=exp(ab(1)+ab(2)*t2(i);,end,plot(t1,x1,r*-,t1,xx1,b+-, t2,x2,g*-,t2,xx2,m+-);,仿真结果说明： 人口增加的指数模型在短期内根本上能比较准确地反映人口自然增长的规律，但长期预测误差很大，需要修正预测模型。,拟合曲线,原始数据曲线,四、人口预测的Logistic模型,人口增长到一定数量后，增长率下降的原因：,资源、环境等因素对人口增长的阻滞作用,且阻滞作用随人口数量增加而变大,假设,r,固有增长率,(,x,很小时,),k人口容量资源、环境能容纳的最大数量,r,是,x,的减函数,例的Logistic模型留给同学们练习,用多项式拟合人口模型,% This program is to predict the model of population by 4-degree polynomial%,%prog42.m%,format long,t1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;,t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;,t=t1;t2;,P1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;,P2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;,P=P1;P2;,n=4;,% The degree of the fitting polynomial%,a,s=polyfit(t1,P1,n);,y=polyval(a,t);,%,a,is the coefficients vector from n-degree to 0-degree%,plot(t,P,r*-,t,y,b+-);,23,a =1.0e+006 *,仿真结果说明, 人口增加的模型用多项式拟合能比较准确地反映人口自然增长的规律，对长期预测具有指导意义。,例2: 海底光缆线长度预测模型,B,2,4,6,8,10,12,14,16,18,20,9,8,6,4,2,0,A,D,C,探测到一组等分点位置的深度数据如下表所示.,25,步探测,从而估计所需光缆的长度,为工程预算提供依据.根本情况如下图.,10.93,10.80,9.81,8.86,7.95,7.95,9.15,10.22,11.29,12.61,13.32,20,19,18,17,16,15,14,13,12,11,10,13.28,12.26,11.18,10.13,9.05,8.02,7.96,7.96,8.96,9.01,深度,(m),9,8,7,6,5,4,3,2,1,0,分点,21个等分点处的深度,(1) 预测通过这条河沟所需光缆长度的近似值.,(2) 作出铺设沟底光缆的曲线图.,解： 用12次多项式函数拟合光缆走势的曲线图如下,仿真结果说明,拟合曲线能较准确地反映光缆的走势图.,The length of the label is L= 26.3809 (m),假设所铺设的光缆足够柔软,在铺设过程中光缆触地走势光滑,紧贴地面,并且忽略水流对光缆的冲击.,% prog45.m This program is to fit the data by polynomial %,format long,t=linspace(0,20,21);,x=linspace(0,20,100);,P=9.01,8.96,7.96,7.97,8.02,9.05,10.13,11.18,12.26,13.28,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93;,a,s=polyfit(t,P,12);,yy=polyval(a,x);,disp(yy=);disp(yy);,plot(x,yy,r*-,t,P,b+-);,L=0;,for,i=2:100,L=L+sqrt(x(i)-x(i-1)2+(yy(i)-yy(i-1)2);,end,disp(The length of the label is L=);disp(L);,1. lsqcurvefit,数据点： xdata=xdata1，xdata2，xdatan，,ydata=ydata1，ydata2，ydatan,用,MATLAB,作非线性最小二乘拟合,Matlab的提供了两个求非线性最小二乘拟合的函数：,lsqcurvefit,和,lsqnonlin,。两个命令都要先建立M-文件，在其中定义函数f(x)，但两者定义f(x)的方式是不同的,可参考例题,.,lsqcurvefit,用以求含参量x（向量）的向量值函数,F(x,xdata)=（F（x，xdata,1,），,，F（x，xdata,n,）,T,中的参变量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,);,(6),x, options,funval, Jacob = lsqcurvefit (fun,x0,xdata,ydata,);,fun是一个事先建立的定义函数,F(x,xdata) 的M-文件,自变量为,x和xdata,说明：,x = lsqcurvefit (fun,x0,xdata,ydata,options);,迭代初值,已知数据点,选项见无,约束优化,lsqnonlin,用以求含参量x（向量）的向量值函数,f(x),=(f,1,(x),f,2,(x),f,n,(x),T,中的参量x，使得,最小。,其中,f,i,（x）=,f,（x，xdata,i,，ydata,i,）,=F(x,xdata,i,)-ydata,i,2. lsqnonlin,数据点： xdata=xdata1，xdata2，xdatan,ydata=ydata1，ydata2，ydatan,输入格式为：,1 x=lsqnonlinfun，x0；,2 x= lsqnonlin fun，x0，options；,3 x= lsqnonlin fun，x0，options，grad；,4 x，options= lsqnonlin fun，x0，；,5 x，options，funval= lsqnonlin fun，x0，；,说明：x= lsqnonlin fun，x0，options；,fun是一个事先建立的定义函数,f(x),的M-文件，,自变量为x,迭代初值,选项见无,约束优化,例6 观测数据点如表所示,x,y,0,3.1,0.1,3.27,0.2,3.81,0.3,4.5,0.4,5.18,0.5,6,0.6,7.05,0.7,8.56,0.8,9.69,0.9,11.25,1,13.17,求三个参数 a, b, c的值, 使得曲线 f(x)=aex+bx2+cx3 与数据点在最小二乘意义上充分接近.,首先编写存储拟合函数的函数文件.,function,f=nihehanshu(x,xdata),f=x(1)*exp(xdata)+x(2)*xdata.2+x(3)*xdata.3,例6 观测数据点如表所示,x,y,0,6,1,求三个参数 a, b, c的值, 使得曲线 f(x)=aex+bx2+cx3 与数据点在最小二乘意义上充分接近.,编写下面的程序调用拟合函数.,xdata=0:0.1:1;,ydata=3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17;,x0=0,0,0;,x,resnorm=lsqcurvefit(nihehanshu,x0,xdata,ydata),编写下面的程序调用拟合函数.,xdata=0:0.1:1;,ydata=3.1,3.27,3.81,4.5,5.18,6,7.05,8.56,9.69,11.25,13.17;,x0=0,0,0;,x,resnorm=lsqcurvefit(nihehanshu,x0,xdata,ydata),程序运行后显示,x =,resnorm =,例6 观测数据点如表所示,x,y,0,6,1,求三个参数 a, b, c的值, 使得曲线 f(x)=aex+bx2+cx3 与数据点在最小二乘意义上充分接近.,说明: 最小二乘意义上的最正确拟合函数为,f,(,x,)= 3e,x,x,2,+ 0.94,x,3,.,此时的残差是: 0.0912.,f,(,x,)= 3e,x,x,2,+ 0.94,x,3,.,拟合函数为:,最小二乘拟合：lsqcurvefit(),例7 假设有一组实测数据,xi 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0,yi 2.3201 2.647 2.9707 3.2885 3.6008 3.909 4.2147 4.5191 4.8232 5.1275,假设该数据可能满足的原型函数为,试求出满足下面数据的最小二乘解的a,b,c,d的值。,例8 用下面一组数据拟合,中的参数a，b，k,该问题即解最优化问题：,MATLAB(fzxec1),function,f=curvefun1(,x,tdata,),f=x(1)+x(2)*exp(-0.02*x(3)*tdata),%,其中 x(1)=a; x(2)=b；x(3)=k;,2输入命令,tdata=100:100:1000,cdata=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 (curvefun1,x0,tdata,cdata),f= curvefun1(x,tdata),F(x，tdata)= ，x=(a，b，k),解法1,.,用命令,lsqcurvefit,3运算结果为：,f =0.0043 0.0051 0.0056 0.0059 0.0061,4结论：,MATLAB(fzxec2),解法,2,用命令,lsqnonlin,f(x)=F(x,tdata,ctada)=,x=（a，b，k）,1,function f=curvefun2(x),tdata=100:100:1000;,cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;,f=x(1)+x(2)*exp(-0.02*x(3)*tdata)- cdata,2输入命令:,x0=0.2,0.05,0.05;,x=lsqnonlin(curvefun2,x0),f= curvefun2(x),函数curvefun2的自变量是x，cdata和tdata是参数，故应将cdata tdata的值写在,3运算结果为,f =1.0e-003 *(0.2322 -0.1243 -0.2495 -0.2413,可以看出,两个命令的计算结果是一样的.,4结论：即拟合得a=,说明：拟合与统计回归,区别与联系,统计回归,线性回归 regress命令,非线性回归,非线性回 归,1确定回归系数的命令：,beta，r，J=nlinfitx，y，model, beta0,2非线性回归命令：nlintoolx，y，model, beta0，alpha,1、回归：,残差,Jacobian矩阵,回归系数的初值,是事先用m-文件定义的非线性函数,估计出的回归系数,输入数据x、y分别为,矩阵和n维列向量，对一元非线性回归，x为n维列向量。,2、预测和预测误差估计：,Y，DELTA=nlpredci（,model, x，beta，r，J,）,求nlinfit 或nlintool所得的回归函数在x处的预测值Y及预测值的显著性为1-alpha的置信区间Y DELTA.,例题的求解,：,2、输入数据：,x=2:16;,y=6.42 8.20 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.80 10.60,10.90 10.76;,beta0=8 2;,3、,求回归系数：,beta,r ,J=nlinfit(x,y,volum,beta0)；,beta,得结果：beta =,即得回归模型为：,To MATLAB(liti41),4、预测及作图：,YY,delta=nlpredci(volum,x,beta,r ,J)；,plot(x,y,k+,x,YY,r),作业与练习,练习1 用给定的多项式，如y=x3-6x2+5x-3，产生一组数据(xi,yi，i=1,2,n),再在yi上添加随机干扰(可用rand产生(0,1)均匀分布随机数,或用rands产生N(0,1)分布随机数)，然后用xi和添加了随机干扰的yi作的3次多项式拟合，与原系数比较。,分别作1、2、4、6次多项式拟合，比较结果，体会欠拟合、过拟合现象。,练习2,用电压V=10伏的电池给电容器充电，电容器上t时刻的电压为 ，其中V,0,是电容器的初始电压， 是充电常数。试由下面一组t，V数据确定V,0,， 。,分别应用非线性最小二乘拟合以及非线性回归命令求解，并作比较，体会统计回归与拟合方法的区别。,练习3 在某海域测得一些点(x,y)处的水深z由下表给出，船的吃水深度为5英尺，估计在矩形区域75，200*-50，150里的哪些地方船要防止进入。,用插值方法作海底曲面图.作出水深小于5的海域范围,即z=5的等高线.,Thank You !,不尽之处，恳请指正！,