用MATLAB进行数据拟合教学提纲.ppt

用MATLAB进行数据拟合拟拟 合合 问问 题题 引引 例例 1 1温度温度t(0C)20.5 32.7 51.0 73.0 95.7电阻电阻R()765 826 873 942 1032已知热敏电阻数据：已知热敏电阻数据：求求60600C时的电阻时的电阻R。设设 R=at+ba,b为待定系数为待定系数拟拟 合合 问问 题题 引引 例例 2 2 t(h)0.25 0.5 1 1.5 2 3 4 6 8c(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)在某种准则下与所有在某种准则下与所有数据点最为接近，即曲线拟合得最好。数据点最为接近，即曲线拟合得最好。+xyy=f(x)(xi,yi)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 的平方和最小的平方和最小。记记问题归结为，求问题归结为，求 a1,a2,am 使使 J(a1,a2,am)最小。最小。线性最小二乘法的求解：预备知识线性最小二乘法的求解：预备知识超定方程组超定方程组：方程个数大于未知量个数的方程组：方程个数大于未知量个数的方程组即即 Ra=y其中其中超定方程一般是不存在解的矛盾方程组。超定方程一般是不存在解的矛盾方程组。如果有向量如果有向量a使得使得 达到最小，达到最小，则称则称a为上述为上述超定方程的最小二乘解超定方程的最小二乘解。线性最小二乘法的求解线性最小二乘法的求解 定理：定理：当当R RT TR R可逆时，超定方程组（可逆时，超定方程组（3 3）存在最小二乘解，）存在最小二乘解，且即为方程组且即为方程组 R RT TRa=RRa=RT Ty y的解：的解：a=(Ra=(RT TR)R)-1-1R RT Ty y 所以，曲线拟合的最小二乘法要解决的问题，实际上就是所以，曲线拟合的最小二乘法要解决的问题，实际上就是求以下超定方程组的最小二乘解的问题。求以下超定方程组的最小二乘解的问题。其中其中Ra=y （3）线性最小二乘拟合线性最小二乘拟合 f(x)=a1r1(x)+amrm(x)中函数中函数rr1 1(x),(x),r rm m(x)(x)的选取的选取 1.1.通过机理分析建立数学模型来确定通过机理分析建立数学模型来确定 f(x)f(x)；+f=a1+a2xf=a1+a2x+a3x2f=a1+a2x+a3x2f=a1+a2/xf=aebxf=ae-bx 2.2.将数据将数据 (xi,yi)i=1,n 作图，通过直观判断确定作图，通过直观判断确定 f(x)：用用MATLAB解拟合问题解拟合问题1 1、线性最小二乘拟合、线性最小二乘拟合2 2、非线性最小二乘拟合、非线性最小二乘拟合用用MATLAB作线性最小二乘拟合作线性最小二乘拟合1.1.作多项式作多项式f(x)=a1xm+amx+am+1拟合拟合,可利用已有程序可利用已有程序:a=polyfit(x,y,m)2.2.对超定方程组对超定方程组可得最小二乘意义下的解。可得最小二乘意义下的解。，用，用3.3.多项式在多项式在x x处的值处的值y y可用以下命令计算：可用以下命令计算：y=polyval y=polyval（a a，x x）输出拟合多项式系数输出拟合多项式系数a=a1,am,am+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=RyMATLAB(zxec1)解法解法1 1用解超定方程的方法用解超定方程的方法2）计算结果）计算结果：=-9.8108 20.1293 -0.03171）输入以下命令：）输入以下命令：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）计算结果：）计算结果：=-9.8108 20.1293 -0.0317解法解法2用多项式拟合的命令用多项式拟合的命令MATLAB(zxec2)例例2 已知观测数据点如表所示已知观测数据点如表所示xy0-0.4470.11.9780.23.280.36.160.47.080.57.340.67.660.79.560.89.480.99.3111.2分别用分别用3次和次和6次多项式曲线拟合这些数据点次多项式曲线拟合这些数据点.x=0:0.1:1y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2plot(x,y,k.,markersize,25)axis(0 1.3-2 16)p3=polyfit(x,y,3)p6=polyfit(x,y,6)编写编写Matlab程序如下程序如下:t=0:0.1:1.2s=polyval(p3,t)s1=polyval(p6,t)hold onplot(t,s,r-,linewidth,2)plot(t,s,b-,linewidth,2)gridx=0:0.1:1y=-0.447,1.978,3.28,6.16,7.08,7.34,7.66,9.56,9.48,9.3,11.2plot(x,y,k.,markersize,25)axis(0 1.3-2 16)p3=polyfit(x,y,3)p6=polyfit(x,y,6)例例3 用切削机床进行金属品加工时用切削机床进行金属品加工时,为了适当地调整为了适当地调整机床机床,需要测定刀具的磨损速度需要测定刀具的磨损速度.在一定的时间测量在一定的时间测量刀具的厚度刀具的厚度,得数据如表所示得数据如表所示:切削时间切削时间 t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm解解:描出散点图描出散点图,在命令窗口输入在命令窗口输入:t=0:1:16y=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.0plot(t,y,*)解解:描出散点图描出散点图,在命令窗口输入在命令窗口输入:t=0:1:16y=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.0plot(t,y,*)a=-0.3012 29.3804hold onplot(t,y1),hold offa=polyfit(t,y,1)y1=-0.3012*t+29.3804例例3 用切削机床进行金属品加工时用切削机床进行金属品加工时,为了适当地调整为了适当地调整机床机床,需要测定刀具的磨损速度需要测定刀具的磨损速度.在一定的时间测量在一定的时间测量刀具的厚度刀具的厚度,得数据如表所示得数据如表所示:切削时间切削时间 t/h030.0129.1228.4328.1428.0527.7627.5727.2827.0刀具厚度刀具厚度 y/cm切削时间切削时间 t/h926.81026.51126.31226.11325.71425.31524.81624.0刀具厚度刀具厚度 y/cm拟合曲线为拟合曲线为:y=-0.3012t+29.3804例例4 一个一个15.4cm30.48cm的混凝土柱在加压实验中的应的混凝土柱在加压实验中的应力力-应变关系测试点的数据如表所示应变关系测试点的数据如表所示1.552.472.933.03已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述,即假设即假设式中式中,表示应力表示应力,单位是单位是 N/m2;表示应表示应变变.2.89已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述,即假设即假设式中式中,表示应力表示应力,单位是单位是 N/m2;表示应表示应变变.解解 选取选取指数函数指数函数作拟合时作拟合时,在拟合前需作变量代换在拟合前需作变量代换,化为化为 k1,k2 的线性函的线性函数数.于是于是,令令即即在命令窗口输入在命令窗口输入:x=500*1.0e-6 1000*1.0e-6 1500*1.0e-6 2000*1.0e-6 2375*1.0e-6y=3.103*1.0e+3 2.465*1.0e+3 1.953*1.0e+3 1.517*1.0e+3 1.219*1.0e+3z=log(y)a=polyfit(x,z,1)k1=exp(8.3009)w=1.55 2.47 2.93 3.03 2.89plot(x,w,*)y1=exp(8.3009)*x.*exp(-494.5209*x)plot(x,w,*,x,y1,r-)已知应力已知应力-应变关系可以用一条指数曲线来描述应变关系可以用一条指数曲线来描述,即假设即假设式中式中,表示应力表示应力,单位是单位是 N/m2;表示应表示应变变.拟合曲线为拟合曲线为:令令则则求得求得于是于是在实际应用中常见的拟合曲线有在实际应用中常见的拟合曲线有:直线直线多项式多项式一般一般 n=2,3,不宜过高不宜过高.双曲线双曲线(一支一支)指数曲线指数曲线二、人口预测线性模型二、人口预测线性模型对于开始提出的实验问题对于开始提出的实验问题,代如数据，计算得代如数据，计算得从而得到人口数与年份的函数关系为从而得到人口数与年份的函数关系为把把x=1999代如，估算出代如，估算出1999年的人口数为年的人口数为 y=1252.1（百万）（百万）12.52亿亿1999年实际人口数量为年实际人口数量为.亿。亿。线性预测模型线性预测模型 英国人口学家英国人口学家MalthusMalthus根据百余年的人口统计资根据百余年的人口统计资料，于料，于17981798年提出了著名的人口自然增长的年提出了著名的人口自然增长的指数增指数增长模型长模型。三、人口预测的三、人口预测的Malthus模型模型基本假设基本假设 :人口人口(相对相对)增长率增长率 r 是常数是常数设设x(t)为为t时刻的人口数量时刻的人口数量,t=0时人口数为时人口数为x0指数增长模型指数增长模型实际中，常用实际中，常用1.1.由前由前100100年的数据求出美国的人口增长年的数据求出美国的人口增长Malthus模型模型。2.2.预测后预测后100100年（每隔年（每隔1010年）的人口状况。年）的人口状况。3.3.根据预测的人口状况和实际的人口数量根据预测的人口状况和实际的人口数量,讨论人讨论人口模型的改进情况。口模型的改进情况。美国美国17901790年年19801980年每隔年每隔1010年的人口记录年的人口记录226.5204.0179.3150.7131.7123.2106.592.076.062.9人口人口(百万百万)1980197019601950194019301920191019001890年份年份50.238.631.423.217.112.99.67.25.33.9人口人口(百万百万)1880187018601850184018301820181018001790年份年份例例解：解：取得最小值取得最小值.其中其中,表示人口数量表示人口数量。表示年份表示年份,解方程组解方程组:即得参数即得参数的值的值.使得使得问题转化为求参数问题转化为求参数 prog41.m%prog41.m%This program is to predict the number of%This program is to predict the number of population%population%format longformat longt1=1790;1800;1810;1820;1830;1840;1850;1860;187t1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;0;1880;t2=1890;1900;1910;1920;1930;1940;1950;1960;197t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;0;1980;x1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.x1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;2;x2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.x2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;3;204.0;226.5;lnx1=log(x1);lnx2=log(x2);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);forfor i=1:10 i=1:10 xx1(i)=exp(ab(1)+ab(2)*t1(i);xx1(i)=exp(ab(1)+ab(2)*t1(i);endendforfor i=1:10 i=1:10 xx2(i)=exp(ab(1)+ab(2)*t2(i);xx2(i)=exp(ab(1)+ab(2)*t2(i);endendplot(t1,x1,r*-,t1,xx1,b+-,t2,x2,g*-plot(t1,x1,r*-,t1,xx1,b+-,t2,x2,g*-,t2,xx2,m+-);,t2,xx2,m+-);a=-49.79535457790735b=0.02859807120038仿真结果表明：仿真结果表明：人口增加的指人口增加的指数模型在短期数模型在短期内基本上能比内基本上能比较准确地反映较准确地反映人口自然增长人口自然增长的规律，但长的规律，但长期预测误差很期预测误差很大，需要修正大，需要修正预测模型。预测模型。拟合曲线拟合曲线原始数据曲线原始数据曲线四、人口预测的四、人口预测的Logistic模型模型人口增长到一定数量后，增长率下降的原因：人口增长到一定数量后，增长率下降的原因：资源、环境等因素对人口增长的阻滞作用资源、环境等因素对人口增长的阻滞作用且阻滞作用随人口数量增加而变大且阻滞作用随人口数量增加而变大假设假设r固有增长率固有增长率(x很小时很小时)k人口容量（资源、环境能容纳的最大数量）人口容量（资源、环境能容纳的最大数量）r是是x的减函数的减函数例的例的LogisticLogistic模型留给同学们练习模型留给同学们练习用多项式拟合人口模型用多项式拟合人口模型%This program is to predict the model of population by%This program is to predict the model of population by 4-degree polynomial%4-degree polynomial%prog42.m%prog42.m%format longformat longt1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;t1=1790;1800;1810;1820;1830;1840;1850;1860;1870;1880;t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;t2=1890;1900;1910;1920;1930;1940;1950;1960;1970;1980;t=t1;t2;t=t1;t2;P1=3.9;5.3;7.2;9.6;12.9;17.1;23.2;31.4;38.6;50.2;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;22P2=62.9;76.0;92.0;106.5;123.2;131.7;150.7;179.3;204.0;226.5;6.5;P=P1;P2;P=P1;P2;n=4;n=4;%The degree of the fitting polynomial%The degree of the fitting polynomial%a,s=polyfit(t1,P1,n);a,s=polyfit(t1,P1,n);y=polyval(a,t);y=polyval(a,t);%a a is the coefficients vector from n-degree to 0-degree%is the coefficients vector from n-degree to 0-degree%plot(t,P,r*-,t,y,b+-);plot(t,P,r*-,t,y,b+-);23a=1.0e+006*-0.00000000000014 0.00000000107892 -0.00000304878595 0.00381927346813 -1.79012132225427仿真结果表明仿真结果表明,人口增加的模型用多项式拟合能人口增加的模型用多项式拟合能比较准确地反映人口自然增长的规律，对长期预比较准确地反映人口自然增长的规律，对长期预测具有指导意义。测具有指导意义。例例2:2:海底光缆线长度预测模型海底光缆线长度预测模型某某一一通通信信公公司司在在一一次次施施工工 中中,需需要要在在水水面面宽宽为为2 20 0m m的的河河沟沟底底沿沿直直线线走走向向铺铺设设一一条条沟沟底底光光缆缆.在在铺铺设设光光缆缆之之前前需需要要对对沟沟 底底 的的 地地 形形 做做 初初B2468101214161820986420ADC探测到一组等分点位置的深度数据如下表所示探测到一组等分点位置的深度数据如下表所示.25步探测步探测,从而估计所需光缆的长度从而估计所需光缆的长度,为工程预算为工程预算提供依据提供依据.基本情况如图所示基本情况如图所示.10.9310.809.818.867.957.959.1510.2211.2912.6113.32201918171615141312111013.2812.2611.1810.139.058.027.967.968.969.01深度深度(m)9876543210分点分点2121个等分点处的深度个等分点处的深度(1)(1)预测通过这条河沟所需光缆长度的近似值预测通过这条河沟所需光缆长度的近似值.(2)(2)作出铺设沟底光缆的曲线图作出铺设沟底光缆的曲线图.解：解：用用12次多项式函数拟合光缆走势的曲线图如下次多项式函数拟合光缆走势的曲线图如下仿真结果表仿真结果表明明,拟合曲拟合曲线能较准确线能较准确地反映光缆地反映光缆的走势图的走势图.The length of the label is L=26.3809(m)假设所铺设的光缆足够柔软假设所铺设的光缆足够柔软,在铺设过程中光缆触地在铺设过程中光缆触地走势光滑走势光滑,紧贴地面紧贴地面,并且忽略水流对光缆的冲击并且忽略水流对光缆的冲击.%prog45.m This program is to fit the data by polynomial%prog45.m This program is to fit the data by polynomial%format longformat longt=linspace(0,20,21);t=linspace(0,20,21);x=linspace(0,20,100);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.2P=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.88,13.32,12.61,11.29,10.22,9.15,7.90,7.95,8.86,9.81,10.80,10.93;0,10.93;a,s=polyfit(t,P,12);a,s=polyfit(t,P,12);yy=polyval(a,x);yy=polyval(a,x);disp(yy=);disp(yy);disp(yy=);disp(yy);plot(x,yy,r*-,t,P,b+-);plot(x,yy,r*-,t,P,b+-);L=0;L=0;forfor i=2:100 i=2:100 L=L+sqrt(x(i)-x(i-1)2+(yy(i)-yy(i-1)2);L=L+sqrt(x(i)-x(i-1)2+(yy(i)-yy(i-1)2);endenddisp(The length of the label is L=);disp(L);disp(The length of the label is L=);disp(L);1.lsqcurvefit1.lsqcurvefit已知数据点数据点：xdata=xdata=（xdata1，xdata2，xdataxdatan n），），ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）用用MATLAB作非线性最小二乘拟合作非线性最小二乘拟合 Matlab Matlab的提供了两个求非线性最小二乘拟合的函数：的提供了两个求非线性最小二乘拟合的函数：lsqcurvefitlsqcurvefit和lsqnonlinlsqnonlin。两个命令都要先建立。两个命令都要先建立M-M-文件文件fun.mfun.m，在其中定义函数，在其中定义函数f(x)f(x)，但两者定义，但两者定义f(x)f(x)的方式是不同的的方式是不同的,可可参考例题参考例题.lsqcurvefitlsqcurvefit用以求含参量用以求含参量x x（向量）的向量值函数（向量）的向量值函数F(x,xdata)=F(x,xdata)=（F F（x x，xdataxdata1 1），），F F（x x，xdataxdatan n）T T中的参变量中的参变量x(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 x（向量）的向量值函数（向量）的向量值函数 f(x)f(x)=(f=(f1 1(x),f(x),f2 2(x),(x),f,fn n(x)(x)T T 中的参量中的参量x x，使得，使得 最小。最小。其中其中 fi（x）=f（x，xdatai，ydatai）=F(x,xdatai)-ydatai 2.lsqnonlin已知数据点：已知数据点：xdata=xdata=（xdata1，xdata2，xdataxdatan n）ydata=ydata=（ydataydata1 1，ydataydata2 2，ydataydatan n）输入格式为：输入格式为：1）x=lsqnonlin（fun，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=lsqnonlinlsqnonlin（fun，x0，options）；）；fun是一个事先建立的是一个事先建立的定义函数定义函数f(x)的的M-文件，文件，自变量为自变量为x迭代初值迭代初值选项见无选项见无约束优化约束优化例例6 已知观测数据点如表所示已知观测数据点如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.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保存为文件保存为文件 nihehanshu.m例例6 已知观测数据点如表所示已知观测数据点如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三个参数求三个参数 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=3.0022 4.0304 0.9404resnorm=0.0912例例6 已知观测数据点如表所示已知观测数据点如表所示xy03.10.13.270.23.810.34.50.45.180.560.67.050.78.560.89.690.911.25113.17求三个参数求三个参数 a,b,c的值的值,使得曲线使得曲线 f(x)=aex+bx2+cx3 与与已知数据点在最小二乘意义上充分接近已知数据点在最小二乘意义上充分接近.说明说明:最小二乘意义上的最佳拟合函数为最小二乘意义上的最佳拟合函数为f(x)=3ex+4.03x2+0.94 x3.此时的残差是此时的残差是:0.0912.f(x)=3ex+4.03x2+0.94 x3.拟合函数为拟合函数为:最小二乘拟合：最小二乘拟合：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)1 1）编写）编写M-M-文件文件 curvefun1.m curvefun1.m 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:1000tdata=100:100:1000cdata=cdata=1e-03*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,*4.54,4.99,5.35,5.65,5.90,6.10,6.26,6.39,6.50,6.59;6.50,6.59;x0=0.2,0.05,0.05;x0=0.2,0.05,0.05;x=lsqcurvefit(curvefun1,x0,tdata,cdata)x=lsqcurvefit(curvefun1,x0,tdata,cdata)f=f=curvefun1(x,tdata)F(x，tdata)=，x=(a，b，k)解法解法1 1.用命令用命令lsqcurvefitlsqcurvefit3 3）运算结果为）运算结果为：f=0.0043 0.0051 0.0056 0.0059 0.0061 f=0.0043 0.0051 0.0056 0.0059 0.0061 0.0062 0.0062 0.0063 0.0063 0.0063 0.0062 0.0062 0.0063 0.0063 0.0063 x=0.0063 -0.0034 0.2542 x=0.0063 -0.0034 0.25424）结论）结论：a=0.0063,b=-0.0034,k=0.2542MATLAB(fzxec2)解法解法 2 用命令用命令lsqnonlin f(x)=F(x,tdata,ctada)=x=（a，b，k）1）编写编写M-M-文件文件 curvefun2.m curvefun2.m 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)-cdata2）输入命令）输入命令:x0=0.2,0.05,0.05;x=lsqnonlin(curvefun2,x0)f=curvefun2(x)函数函数curvefun2的自变量是的自变量是x，cdata和和tdata是已知参数，故应是已知参数，故应将将cdata tdata的值写在的值写在curvefun2.m中中3 3）运算结果为）运算结果为 f=1.0e-003*(0.2322 -0.1243 -0.2495 -0.2413 f=1.0e-003*(0.2322 -0.1243 -0.2495 -0.2413-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792-0.1668 -0.0724 0.0241 0.1159 0.2030 0.2792 x=0.0063 -0.0034 0.2542 x=0.0063 -0.0034 0.2542可以看出可以看出,两个命令的计算结果是相同的两个命令的计算结果是相同的.4）结论）结论：即拟合得即拟合得a=0.0063 b=-0.0034 k=0.25420.0063 b=-0.0034 k=0.2542说明：拟合与统计回归区别与联系统计回归线性回归（regress命令）非线性回归非线性回非线性回 归归（1）确定回归系数的命令：beta，r，J=nlinfit（x，y，model,beta0）（2）非线性回归命令：nlintool（x，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得结果：