matlab数据拟合实验

1、实验八matlab数据拟合实验在科学实验和生产实践中，往往需要从一组实验数据(x.，x)a = h2,/?)中，寻找变 量兀和y之i'可的函数关系y = /(x)的某种近似表达式$(兀)。上一个实验中介绍的插值方法 可以构造一个插值函数逼近己知函数。但是，一般來说，给定的实验数据(斥,牙) (d = l,2,仍的数量较大，且由于观测误差的原因，准确度不一定高，其至在个别点有很 大的谋差，形象地称之为“噪声”。如果用插值法來求y = fm的近似表达式，要使$(x) 满足插值条件，势必将“噪声”带进近似函数s(x),因而不能较好地描绘y = f(x)o为了 尽可能减少这种观测误差的影响，本

2、次实验我们介绍数据拟合和关的方法及其matlab实现。8.1实验目的1、熟悉w. matlab !«各种常见的拟合：8.2实验基本知识2、能够灵活编程來解决数据拟合的实际印数据的拟合主要分为曲线拟合(curve fitting)曲面拟合(surface fitting)0它试图从散点数 据中拟合出-个有规律的解析式世该解析音的某些資数(系数)是不运的未知量数据拟合 的h的就是要根据散点数据在“最小误差"的意义下确定出解析式中这些不定参数。举一个 简单的叶子，对于条直线y = kx + bi ；该式屮有比和b两个未知参数需要求岀。由基本的 数学知识我们可以知道，只要有两个点就

3、可以确疋出这两个参数但是若有更多的点，比如 实验数据往往仃很多点,这些点山于有j赛差并不一足在一条直线上，需要找'漆直线离这吐 点“最近尸这就悬嚟从上面的例？i们可拟两標兄(1) 点数(已知数据数目)即方程的个数要大于待求参数的个数；(2) 方程所代表的曲线、曲面等并不一定通过这些已知的点。在上一段中，“最近的定义冇很多种，不同的定义对应着数据拟合的不同准则。下面我 们介绍常见的儿种数据拟合的准则。切比雪夫(chebyshev)近似准则给定某种函数类型y = /(%)和m个数据点(x.，x)的一个集合，对整个集合极小化最人 绝对偏差卜厂/(兀)|，即确定函数类型y = f(x)的参数从

4、而极小化数量：maxinum |y. -/(x.)| i = 1,2,zn这种准则通常被称为切比雪夫(chebyshev)近似准则。这种准则在实际应用中通常很复 杂，应用这一准则所产生的最优化问题通常可以表示为线性规划问题，这可能需要高级的数 学方法，或者要用计算机的数值算法。所以在对冇限的数据点拟合一条曲线时不常使用该准 则，然而当极小化最大绝对偏差很重要的时候仍应考虑这一准则。另外，用一个简单函数代替在一个区间上定义的另一个复杂函数时，构成该准则的原则是极其重要的，在该区间上两 个函数间的最人差异必须达到最小，因此这一准则在函数逼近问题中具有很重要的应用。极大化绝对偏差之和极小化绝对偏差之

5、和准则可以归纳为：给定某种函数类型y = f(x)和加个数据点 (再)的集合，极小化绝对偏差|x-/(%,-)|的和，也就是确定函数类型y = /(x)的参数, 极小化：/=1如果令尺=卜厂/(忑)|,心1,2,，加，代表每一个绝对偏差，那么该准则可以解释成将一条由数量k加在一起构成的直线的长度极小化。由于这一准则里出现了绝对值，这个和式的各种微分不是连续的，要解决这个最优化问题时，将该和式对每个未知参数进行求导 时，问题会变得不可解，因此在常见的数据拟合问题中，该准则也不常用。但是需要注意的e特殊的场景卜,变的9前ib i相同的记号,ll'uj题是确定函数x=1用此方普解决产生的最优

6、化问题仅需使用儿个变量的演算，所以容易普及，matlab软 件中提供的数据拟合方法也基本都是基*该准则的。数据拟合根据0变量的及选収的 拟合函数的形式，川以分另；三元线性拟合、-見非线性拟合、多元拟合等不同的问题场最。 卜両我们就分别介绍这几种类型问题下最小二乘拟合的matlab实现，并在最后介绍matlab 曲线拟也工具箱的使用。8.2.1 一元线性拟合多项式拟合多项式拟合是线性拟合中最常见的形式，它的目标是找出一组多项式系数 即心1,2,+ 1,使得多项式j (x) = axx +a2x + + g“x + g“+能够在最小二乘的意义下最好地拟合原始数据。多项式拟合可以通过matlab提供

7、的 polyfit()函数实现。该函数的调用格式为：p/s=polyfit(x,y/n)其中，x,y是输入的实验数据，n是拟合计算用到的多项式的次数，返回值p是拟合多 项式的系数，s是用來估计误差和预测数据的结构体。与polyfit()函数配合使用的函数是polyval(),这个函数根据拟合出来的多项式系数p计 算给定数据x处的y值。其调用格式如下：y=polyval(p,x)其中x是给定的需要计算拟合值的向量，p是polyfit()函数的返回值，返回值y是根据p 计算出来的x处的多项式的值。例8j利用函数y = /+2兀+ 3,在区间xe-2,l产牛32个数据点，并对y值加入适当的随机噪声，

8、然厉用二次多项式拟合得到的数据，观察拟合效呆。编写下面的程序test81.m0clc;clear all; a=l;b=2;c=3;x=-2:0.1:l; %产生 x 值y=a*x.a2+b*x+c;%产生y的精确值figureplot(x,y, 丫) %对已知的精确值作图hold on;yrandom=y+(rand(size(x)-0.5)*0.3; %给 y 值加入随机噪声 ployrandom；*') %对要进彳亍拟合的实验数据作图 p,s=polyfit(x/yrandomz2);%对实验数据进行一.次多项式拟合，得到多项式系数 yfit=polyval(p,x);plot(

9、x,yfit/-kd')legendc精确值t实验数据t拟合值')p运行该程序可彳譬咿果p =0.96401.98013.0171当然師;我们的实验数据是加入随机噪声的匸所以该程序每次运行的结果会略有不同, 但能发现拟合效果还是相为好的。另外，拟合的效果如图81所示。：:丿图8-1多项式拟合效果图已知函数线性组合的拟合方法实际上，多项式函数的每一项的都可以被看做是一个函数，而多项式函数就是这些函数 的线性组合。对于更一般的情形，假设已知某函数的线性组合为：/w =(x) + c2f2 (x) + c3f3(x) + + ctifn (x)其中,fl （x） , f2 （x）,，

10、九（兀）为已知函数,c, c2, , c”为待定系数,这时假设已经测出数据（坷，开），（兀2，儿），（九，九），则可以建立其如卜的线性方程: ac = y£3)九（西厂具中，a =z(x2)/2(2) 九(兀2)，v =丿（心），2（心）九（兀）_.xw.且c = q,a?,。所以该方程的最小二乘解为c = /ly o例8-2假设测出了一组数据（兀，必），由下面的表格8-1给出，r已知函数的原型为y（x） = cj + c2exx +c3 cos（-2x04" +c4x2,试用已知数据求出待立系数c;的值。表8jx. 00.20.40.70.90.920.991.21.41

11、.481.5)2.882.25761.96831.95282.08622.1092.1979 2.5409 2.9627 3.155 3.2052编写不面叭程序test82.m0clc;clear all;x=00.20.40.70.90.92y=2.882.25761.96831.95282.08622.109a=ones(size(x) exp(-3*x) cos(-3*x).*exp(-4*x) x.a2;c=ay;cl=c' x0=0:0.01:1.51;al=ones(size(xo) exp(-3*x0) cos(-3*x0).*exp(-4*x0) xo.a2;yl=al

12、*c;plot(xo,yl,x,y,w)图82已知函数线性组介的拟合效果图f该无呈序可以得出拟合曲线和已知数据b如图82所示。8.2.2 一元非线性拟合在上面的例&2屮，我们能看到虽然要拟介的曲数形式是非常复杂的，但是在复杂的函数形式如指数函数、摹函数的指数上没冇未知参数，所以上面的问题会比较简单，可以直接 求解线性方程纟h.得到结果。但在实际的问题场景中，如果我们己知的函数形式为y(x) = aehx 这样在幕指数上带有未知数的形式，简单的线性求解就无能为力了。通常这种情况下，我们 有两种方法可以解决，一是将非线性函数进行适当的变换变成线性形式，另一种是有matlab 的优化工具箱中

13、提供的lsqcurvefit()函数进行求解。非线性拟合化为线性拟合对于某些非线性的变量关系可以通过简单的变量替换使之直线化,这样就可以按最小二 乘法原理求出替换厉变量的直线方程。在实际中常利用此直线方程绘制数据的标准曲线，同 时根据需要可将此直线方程还原为劇线方程，实现对数据的曲线拟合。常用的非线性函数线性化的方法如下：指数函数y = aebx,当b>0时，y随兀增大而增大；当bvo时，y随x增大 而减少。可以两边取对数，得lny = lnd + bx。因此当以in)，和x绘制的散点 图呈直线趋势时，可考虑采用指数函数來描述y与兀之间的非线性关系，in a和3、b分别为截距和斜率。对数

14、函数y = a + bnx f当b>随x增大而增大，先当b<0 时，y随x增人而减少，先快后慢c因此，当以y和巾工绘制的散点图呈直线 趋势时，可考虑采用对数函数描述y-mxz间的非线性关系，式中的b和a分 3为誉和截距。数y = 6zx当>0时，y随x增大而增大；当bvo时，男x增大而减少。可以两边取对数，得到iny = ln + /?lnxo因此，当以lny和in兀绘制的散点图呈rr线趋势吋，可考贱釆用幕函 in a和方分为截距和斜率。例83已知数掘如表82所示，试观察数据给h洽适 下的最佳拟合参数，观察拟合效果。农8-2xi1.22.84.3可的非线性关系,1乘准则5.4

15、6.82.111.528.141.972.391.4开" 首先编写程序test83 l.m,将己知数据在不同朋标轴下的散点图做出来。clc;clear all;x=1.2 2.8 4.3 5.4 6.8 7.9; y=2.1 11.5 28.1 41.9 72.3 91.4; subplot(2,2,l) plot(xzy；d')titl*正常坐标系j;subplot(2,2,2) semilogx(x/y/'d,) title(*lnx-y') subplot(2,2,3) semilogyfx'd')title('x-lny'

16、;)subplot(2,2,4) loglog(x/y；d,) title('lnx-lny')得到的结果如图83所示。图8-3例83数据在不同坐标轴下的效果图通过图像我们很清楚的能看出來，以in兀为横坐标，以in y为纵坐标绘制的曲线更接近于直线，所以我们考虑选用幕函数y =作为拟合函数，下面我们用非线性拟合线性化j的思想來编写程序test83_2.m來蠢待定参数clc;clear all;x=1.2 2.8 4.3 5.4 6.8 7.9;y=2.1 11.5 28.1 41.9 72.3 91.4;logx=log(x);logy=log(y);p=polyfit(log

17、xjogy,l);a=exp(p(2);b=p(l)；x0=l:0.1:10;yo=a*xo.ab;plot(x,y,'*：xo,yo)legend('原始数据点t拟合曲线')；运行结果可以首先求的in y和in x的关系为in y = 0.374331+ 2.01486 in x进而,通过转化可以得到y = 1.45402%20'486,从而得到拟合曲线，并画图，可得到结果如图8-4所示，可以发现拟合效果还是令人满意的。160140120100&0甸4020678910图8-4例8-3的拟合效果图优化工具箱的lsqcurvefit()函数蛊要注意的是，

18、在上面的例子中，我们虽然利用非线性拟合线性化的方法，将一个非线 性的问题转化成了一个简单的线性问题，并简单的利用多项式拟合的polyfit()命令得到了拟 介参数。但是这种做法得到的拟介参数是对变换后的线性多项式的最小二乘拟介，并不是对 原始函数的最小二乘拟介，在很多问题场景下，这种做法得到的结果可能在变换冋原始函数数之后含兰生很大的误差。所以在解决实际聖非线性拟合问題的时候，我们可以通过对实际数 理进行一些变换；去观察数据可能满足的函数规律，但呦实际进行求解参数的时候，瓜量i 还是利用优化工具箱中提供的lsqcurvefit()函数或者甚至有时候需要自己运用-些智能算法 來求解这类最优化问题

19、。这里我们介甕1下matlab自带的lsqcurvefit()函数的用法， 的调用格式为：p=lsqcurvefit(f un ,po,xdata,ydata)其中，fun为非线性函数模型，p0为最优化的初始值，xdata、ydata是原始输入输出数据向量；返回值p为拟合输出的未知参数。 _公,具体见表8-2,试利用拟合方法确定参数a,b,c。例84在一次传染病屮，己知t时刻的染炳久数/(f)满足如卜模型)想一次傷共部门每曙天记录 農8-3天数10510传染者0.20.40.515302.43.13540450.9编写下而的程序test84.m,调用isqcurvefitf)函数來解决该问题。

20、3.84.14.2504555604.44.5clc;clear all;t=0:5:60;y=0.2 0.4 0.5 0.9 1.5 2.4 3.1 3.8 4.1 4.2 4.5 4.4 4.5;f=(bzx) l./(b(l) + b(2)*exp( b(3)*x ); %定义拟合函数模型 x=lsqcurvefit(f,0.5z10,0,t,y); %调用 isqcurvefitf)函数进行最小二乘拟合求解 yfit=f(x,t);plotyjko'xyfit) %绘制图形 legend!1原始数据点丁拟合曲线，)得到的结果如图8-5所示。图8-5 lsqcurvefit()0

21、数非线性拟合效果图8.23多元拟合前面我们比较详细的介绍了 matlab中实现一元拟合的方法，而我们在冇些实际问题中 自变量的数量可能会冇两个或者多个，这时候问题就变成了多元拟合。多元拟合同样存在线 性拟合和非线性拟会的嚳容査多醴性劳合，matlab最优化工具箱中提msqlinoa 求解这类问题。对于多元非线性拟合，matlab拟合工具箱中没有提供函数來解决，有兴趣 的学生可以根据最力匸乘的原理石己來编写程序实现,/外l种替代的方法是使用matlab 翁计工具箱屮提供的求解多元非线性回归的函数nlinfit()来解决这类问题。卜面我们就分别 结合实例来介年吟两个函数的用法。例8-5己知喝观测数

22、拯，如表8-4虑滋进行攻酩f。 x0.10.20.30.40.5_0.60.70.80.9y1.81.71.61.51.41.21b|850.67z55.25.8635.7,5.96.zs表8-4编行下面的稈序test85.m,调用isqlinq函数来解决该问题。clc;clear all;x=0.1:0.1:0.9j,;y=1.8 1.7 1.6 1.5 1.4 1.2 11.085 0.67*;z=5 5.2 5.8 6.3 5.7 5.9 6 7 7'm=x,y/ones(size(x); %构造调用isqlin函数所需的m矩阵 x=lsqlin(m,z);%调用lsqlin()

23、旳数求解二元线性拟合问题xgrid ygrid=meshgrid(x,y);zl=x(l)*xgrid+x(2)*ygrid+x(3);mesh(xgrid,ygrid,zl) %绘制网格曲面图hold on;plot3(x/y/z；kd')%绘制原始数据点alpha(0.5) %设置透明程度hold off;得到的结果如图8-6所示。9图8-6二元线性拟合效果图例8-6反应动力学中的hougen-waston模型是非线性模型的一个经典例子，其模型如下：_01兀2-花/05y 1 + 02兀1+03兀2+0內其中，y为反应速率，其三个决定因素分别为壬（氢气），尢2 （n戊烷），兀3（异

24、戊烷），表8-5是一组实验数据式建立回归表8-5知参数加=1,2氢气t4702854704704701001001 47qj100100100285、285n 戊烷3008030080801908019030030080300190异戊烷10101201201010、65655412012010120反应速率8.55噴794.820.0217514.39,2.544.35138.50.0511.32噸13编写如下的程序test86.m。data=470 285 470 470 470 100 100 470 100 100 100 285 285;300 80 300 80 80 190 80

25、 190 300 300 80 300 190;10 10 120 120 10 10 65 65 54 120 120 10 120*;%反应物数据 rate=8.55 3.79 4.82 0.02 2.75 14.39 2.54 4.35 13 8.5 0.05 11.32 3.13： %反应速率数据 f=(b,x)( b( 1 )*x(:,2)-x(:,3)/b(5)./( 1 +b(2)*x(:,l)+b(3)*x(:,2)+b(4)*x(:,3); %定义拟合函数 b0=l;0.1;0.2;0.1;2; %参数求解的初始值 betafit=nlinfit(data,rate,f,bo

26、) %调用 nlinfit 函数求解得到如下的结果：»betafit =1.25260.06280.04000.11241.19148.2.4 matlab拟合工具箱简介为了使用八更方便地进行曲线拟合，matlab还提供了一个功能强人的曲线拟合工具箱 curve fitting toolbox,它是以图形用户界而的互动操作展示给用户的，而j1最后还可以 将该图形界面操作牛成的m函数文件保存以备后用。下面我们就以例8-4为例说明该拟合 工具箱的用法。第一步：输入观测点数据在matlab命令窗口中输入如下语彳山» t=0:5:60;» y=0.2 0.4 0.5 0.

27、9 1.5 2.4 3.1 3.8 4.1 4.2 4.5 4.4 4.5;第二步：启动曲线拟合工具箱对于1111线拟合工具箱的启动matlab提供了好多方法，如直接在命令窗口输入cftool 或者单击 matlab 界lllj左下角的 start按钮，再执行操作 toolboxes->curve fitting->cftool 即可。进入拟合工具箱初始界面，如图8-7所示。图8-7 具箱create data set i:时会在初始界面中绘i 皿 i g i i 匕 j 1111i第三步：数据导入与拟合，u单,|；4 data _l按钮,弹出“data”对话框,分别由xdata和

28、ydata的下拉菜单导入create data set数据t,y,修改数据集名data set name为“初始数据点”，然后单注 钮，创建一个新蛟据菜"肓帖f我钮"close屈护iata”对话框;这h 制出数据的散点图,如图88所示。图&8初始数据点的散点图单击i f咖g_j按钮,弹出“fitting”对话框，如图89所示，该对话框分为两部分,上面为 “fit editor”，下面为 “table of fits”。| mkhi" f 傑qo w _ ji 二【 m> 】图8-9 u fitting”对话框单击i加丄£泄按钮并将“fit

29、name"该为“拟合曲线”，然后通过下拉菜单“type offit”选择拟合曲线的类型，该工具箱提供的拟合曲线的类型有以下几种。custom equations:用户自定义的函数类型。exponential: e 指数函数。fourier:傅里叶函数。gaussian:高斯函数。interpolant:插值函数。 polynonmial:多项式函数： power:幕函数。rational;有理函数。 smoothing spline:平滑样？sun of sin functions:正弦函数类。weibull:威布尔函数。在每个函数类型对应的列表框中都会显示函数的解析式形式或相关说明

30、，以更好地选 择拟合函数类型。另外，还可以通过章|审尺讥按钮设置拟合算法、待估参数 o 如果选择"custom equations单击“ new customequation”对话框,（如图8j0所示。.图 8-10 “ new custom equation ” 对话楼该对话框包含“linear equations”线性等式和“general equations"构造等式两种自 定义函数的方法，下面分别介绍这两种方法。线性等式实现：首先单击一按钮,获得适当的输入项，再在相应的输入框中输入待处系数和拟合基函数。构造等式实现：直接在相关的输入框中输入拟合函数。函数定义好之厉单

31、击“0k”按钮退出该对话框，单-击“apply”按钮即可在“results”中 显示待定参数的数值及其他和关信息。刃外，在初始图形界面上也直接绘制出拟合曲线图。 在拟合线上单击任意处即可得到该点的朋标，右击可以改变线的颜色、类型等，执行"view->residuals->line plot”操作可以绘制出其残差图。在木例中采用构造等式实现方式，在 v 对应的输入框中输入"l./(a*exp(.b*x)+c)单 击“0k”按钮退出该对话框，再单击“apply”按钮，则在“results”框中得到如下信息:general model:f(x) = l./(a*exp(-b*x)+c)coefficie nts (with 95% con fide nee boun ds):a = b = c =7.38 (5.18,9.579)0.144 (0.1312,0.1568)0.2214 (0.2162,0.2266)goodness of fit:sse: 0.06175r-square: 0.9983adjusted r-square: 0.9979rmse: 0.0