《新编MATLAB&Simulink自学一本通》课件第19章 回归分析.ppt

1、2023-5-5 一元线性回归一元线性回归 一元非线性回归一元非线性回归 多元线性和广义线性回归多元线性和广义线性回归 多元非线性回归多元非线性回归 数据拟合的其他数据拟合的其他MATLAB函数函数主要内容：主要内容：2023-5-5第一节第一节 一元线性回归一元线性回归2023-5-5【例例19.1-1】现有全国现有全国31个主要城市个主要城市2007年的气候情况观测数据，年的气候情况观测数据，如下表所示。试根据表中如下表所示。试根据表中31组观测数据研究年平均气温和全年组观测数据研究年平均气温和全年日照时数之间的关系。日照时数之间的关系。城城 市市年平均气温年平均气温单位：单位：年极端最高

2、气温年极端最高气温单位：单位：年极端最低气温年极端最低气温单位：单位：年均相对湿度年均相对湿度单位：单位：%全年日照时数全年日照时数单位：小时单位：小时全年降水量全年降水量单位：毫米单位：毫米序序号号北京14.037.3-11.7542351.1483.91天津13.638.5-10.6612165.4389.72石家庄14.939.7-7.4592167.7430.43太原11.435.8-13.2552174.6535.44呼和浩特9.035.6-17.6472647.8261.25沈阳9.033.9-23.1682360.9672.36兰州11.134.3-11.9532214.1407

3、.928西宁6.130.7-21.8572364.7523.129银川10.435.0-15.4522529.8214.730乌鲁木齐8.537.6-24.0562853.4419.5312023-5-5 令令x表示年平均气温，表示年平均气温，y表示全年日照时数。由于表示全年日照时数。由于 x和和 y均为一维变量，可以先从均为一维变量，可以先从 x和和 y的散点图上直观的观察它们之的散点图上直观的观察它们之间的关系，然后再作进一步的分析。间的关系，然后再作进一步的分析。ClimateData=xlsread(examp19_1_1.xls);x=ClimateData(:,1);y=Clima

4、teData(:,5);plot(x,y,k.,Markersize,15);xlabel(年平均气温年平均气温(x);ylabel(全年日照时数全年日照时数(y);%计算计算x和和y的线性相关系数矩阵的线性相关系数矩阵R R=corrcoef(x,y)一、一、数据的散点图数据的散点图2023-5-568101214161820222426500100015002000250030003500年 平 均 气 温(x)全年日照时数(y)0.7095xy 2023-5-5二、二、调用调用regress函数作一元线性回归函数作一元线性回归1.p元广义线性回归模型元广义线性回归模型1112121011

5、121222212211221()()()1()()()1()()()ppppnnpnppnnyXf xfxfxyf xfxfxyf xfxfxy矩阵形式矩阵形式01 11222()()(),1,iiippipiyf xfxfxin2023-5-52.regress函数的用法函数的用法nYYYY.21pb.10 b =regress (Y,X )1112121121222211221()().()1()().().1()().()ppppnnpnpf xfxfxf xfxfxXf xfxfx2023-5-5回归系数的区间估计回归系数的区间估计残差残差用于检验回归模型的统计量，用于检验回归模型的

6、统计量，有有4个数值：判定系数个数值：判定系数r2、F 统计量观测值、检验的统计量观测值、检验的p值值、误差方差的估计误差方差的估计置信区间置信区间 显著性水平显著性水平（缺省时为（缺省时为0.05）b,bint,r,rint,stats=regress(Y,X,alpha)2.regress函数的用法（续）函数的用法（续）2023-5-53.调用调用regress函数作一元线性回归函数作一元线性回归 ClimateData=xlsread(examp19_1_1.xls);x=ClimateData(:,1);y=ClimateData(:,5);xdata=ones(size(x,1),1

7、),x;%调用调用regress函数作一元线性回归函数作一元线性回归 b,bint,r,rint,s=regress(y,xdata);yhat=xdata*b;plot(x,y,k.,Markersize,15);hold on;plot(x,yhat,linewidth,3);xlabel(年平均气温年平均气温(x);ylabel(全年日照时数全年日照时数(y);legend(原始散点原始散点,回归直线回归直线);3115.476.962yx2023-5-54.残差分析残差分析 figure rcoplot(r,rint)通过对残差和残差的置信区间进行分析，可以看出原始通过对残差和残差的置

8、信区间进行分析，可以看出原始数据中是否存在异常点，若残差的置信区间不包括数据中是否存在异常点，若残差的置信区间不包括0点，可点，可认为该组观测为异常点。认为该组观测为异常点。51015202530-1500-1000-500050010001500Residual Case Order PlotResidualsCase Number68101214161820222426500100015002000250030003500年平均气温(x)全年日照时数(y)原始数据散点原始数据回归直线剔除异常数据后回归直线2023-5-55.剔除异常值后重新回归剔除异常值后重新回归 xt=x(y1250);

9、yt=y(y1250);xtdata=ones(size(xt,1),1),xt;b,bint,r,rint,s=regress(yt,xtdata)从残差图可以看出有从残差图可以看出有4条线段（红色虚线）与水平线条线段（红色虚线）与水平线 没有没有交点，它们对应的观测序号分别为交点，它们对应的观测序号分别为22、23、24和和26，也就是，也就是说这说这4组观测对应的残差的置信区间不包含组观测对应的残差的置信区间不包含0点，可认为这四点，可认为这四组观测数据为异常数据。剔除异常数据后重新作回归。组观测数据为异常数据。剔除异常数据后重新作回归。2983.863.628yx2023-5-5三、三

10、、调用自编调用自编reglm函数作一元线性回归函数作一元线性回归1.reglm函数的用法函数的用法 reglm(y,X)stats=reglm(y,X)stats=reglm(y,X,model)stats=reglm(y,X,model,varnames)参数说明请参考课本。参数说明请参考课本。2023-5-52.调用调用reglm函数作一元线性回归函数作一元线性回归 ClimateData=xlsread(examp19_1_1.xls);x=ClimateData(:,1);y=ClimateData(:,5);varname=x;reglm(y,x,varname);2023-5-53

11、.结果结果-方差分析表方差分析表-方差来源方差来源 自由度自由度 平方和平方和 均方均方 F值值 p值值回归回归 1.0000 4316961.1182 4316961.1182 29.3884 0.0000残差残差 29.0000 4259914.3811 146893.5993总计总计 30.0000 8576875.4994 均方根误差均方根误差(Root MSE)383.2670 判定系数判定系数(R-Square)0.5033 因变量均值因变量均值(Dependent Mean)1965.1742 调整的判定系数调整的判定系数(Adj R-Sq)0.4862-参数估计参数估计-变量变

12、量 估计值估计值 标准误标准误 t值值 p值值常数项常数项 3115.3773 223.0588 13.9666 0.0000 x -76.9616 14.1967 -5.4211 0.00002023-5-5四、四、调用调用robustfit函数作稳健回归函数作稳健回归1.robustfit函数的用法函数的用法 b=robustfit(X,Y)b,stats=robustfit(X,Y)b,stats=robustfit(X ,Y,wfun,tune,const)加权函数加权函数12nyyy调整常数调整常数是否显示常是否显示常数项的标示，数项的标示，取值为取值为on或或 off2023-5-

13、52.调用调用robustfit函数作一元线性回归函数作一元线性回归 ClimateData=xlsread(examp19_1_1.xls);x=ClimateData(:,1);y=ClimateData(:,5);%利用利用robustfit函数作稳健回归函数作稳健回归 b,stats=robustfit(x,y)b=1.0e+003*3.0348 -0.06833034.868.3yx68101214161820222426500100015002000250030003500年平均气温(x)全年日照时数(y)原始数据散点regress函数对应的回归直线robustfit函数对应的回归

14、直线2023-5-53.regress和和robustfit函数结果对比函数结果对比 xsort,id=sort(x);ysort=y(id);xdata=ones(size(xsort,1),1),xsort;b1=regress(ysort,xdata);yhat1=xdata*b1;b2=robustfit(xsort,ysort);yhat2=xdata*b2;plot(x,y,ko)hold on plot(xsort,yhat1,r-,linewidth,3)plot(xsort,yhat2,linewidth,3)legend(原始数据散点原始数据散点,.regress函数对应的

15、回归直线函数对应的回归直线,.robustfit函数对应的回归直线函数对应的回归直线);xlabel(年平均气温年平均气温(x)ylabel(全年日照时数全年日照时数(y)2023-5-5第二节第二节 一元非线性回归一元非线性回归2023-5-5【例例19.2-1】头围头围(head circumference)是反映婴幼儿大脑和颅骨是反映婴幼儿大脑和颅骨发育程度的重要指标之一，对头围的研究具有非常重要的意义。发育程度的重要指标之一，对头围的研究具有非常重要的意义。笔者研究了天津地区笔者研究了天津地区1281位儿童（位儿童（700个男孩，个男孩，581个女孩）的个女孩）的颅脑发育情况，测量了年

16、龄、头宽、头长、头宽颅脑发育情况，测量了年龄、头宽、头长、头宽/头长、头围和头长、头围和颅围等指标，测得颅围等指标，测得1281组数据，年龄跨度从组数据，年龄跨度从7个星期到个星期到16周岁，周岁，试根据这试根据这1281组数据建立头围关于年龄的回归方程。组数据建立头围关于年龄的回归方程。2023-5-5 令令x表示年龄，表示年龄，y表示头围。由于表示头围。由于 x和和 y均为一维变量，均为一维变量，可以先从可以先从 x和和 y的散点图上直观的观察它们之间的关系，然后的散点图上直观的观察它们之间的关系，然后再作进一步的分析。再作进一步的分析。HeadData=xlsread(examp19_2

17、_1.xls);x=HeadData(:,4);y=HeadData(:,9);plot(x,y,k.)xlabel(年龄年龄(x)ylabel(头围头围(y)一、一、数据的散点图及备选方程数据的散点图及备选方程1.散点图散点图2023-5-50246810121416354045505560年龄(x)头围(y)2023-5-5 负指数函数：负指数函数：2.备选方程备选方程231xye 双曲线函数：双曲线函数：123xyx 幂函数：幂函数：312()yx Logistic曲线函数：曲线函数：31()21xye 对数函数：对数函数：123ln()yx2023-5-5二、二、调用调用nlinfit

18、函数作一元非线性回归函数作一元非线性回归1.nlinfit函数的用法函数的用法 beta,r,J,COVB,mse=nlinfit(X,y,fun,b0,options)1212,;,1,2,ikiiipiyfxxxin 未知参数未知参数事先用事先用m-文件文件定义的非线性定义的非线性函数函数111212122212ppnnnpxxxxxxXxxx12nyyyy12.k回归系数初值回归系数初值优化属性设置优化属性设置雅可比矩阵雅可比矩阵残差残差2023-5-52.选择合适的理论回归方程选择合适的理论回归方程 负指数函数：负指数函数：231xye2023-5-53.调用调用nlinfit函数作一

19、元非线性回归函数作一元非线性回归 HeadData=xlsread(examp19_2_1.xls);HeadData=sortrows(HeadData,4);x=HeadData(:,4);y=HeadData(:,9);%理论回归方程对应的匿名函数理论回归方程对应的匿名函数 HeadCir2=(beta,x)beta(1)*exp(beta(2)./(x+beta(3);options=statset;options.Robust=on;%调用调用nlinfit函数作非线性回归函数作非线性回归 beta,r,J,COVB,mse=nlinfit(x,y,HeadCir2,53,-0.26

20、04,0.6276,options);0.25950.760452.3766xye2023-5-54.绘制一元非线性回归曲线绘制一元非线性回归曲线 yhat=HeadCir2(beta,x);figure plot(x,y,k.)hold on plot(x,yhat,linewidth,3)xlabel(年龄年龄(x)ylabel(头围头围(y)legend(原始数据散点原始数据散点,.非线性回归曲线非线性回归曲线)0246810121416354045505560年 龄(x)头围(y)原 始 数 据 散 点非 线 性 回 归 曲 线2023-5-55.参数估计值的置信区间参数估计值的置信区

21、间 在调用在调用nlinfit函数得到未知参数的估计值函数得到未知参数的估计值beta、残差向量、残差向量r、雅可比矩阵、雅可比矩阵J、未知参数的协方差矩阵的估计、未知参数的协方差矩阵的估计COVB和误和误差方差差方差 的估计的估计mse后，还可以调用后，还可以调用nlparci函数函数计算参数估计计算参数估计值的置信区间。值的置信区间。beta,resid,J,Sigma=nlinfit(X,y,fun,b0)ci=nlparci(beta,resid,covar,Sigma)或或 ci=nlparci(beta,resid,jacobian,J)%求参数估计值的求参数估计值的95%置信区间

22、的第置信区间的第1种方式种方式 ci1=nlparci(beta,r,covar,COVB,alpha,0.05)%求参数估计值的求参数估计值的95%置信区间的第置信区间的第2种方式种方式 ci1=nlparci(beta,r,jacobian,J,alpha,0.05)2023-5-56.头围平均值的置信区间和观测值的预测区间头围平均值的置信区间和观测值的预测区间 对于对于 x（年龄）的一个给定值（年龄）的一个给定值 x0，相应的，相应的 y（头围）是一（头围）是一个随机变量，具有一定的分布。个随机变量，具有一定的分布。x 给定时给定时y 的总体均值的区的总体均值的区间估计称为平均值（或预测

23、值）的置信区间，间估计称为平均值（或预测值）的置信区间，y的观测值的的观测值的区间估计称为观测值的预测区间。区间估计称为观测值的预测区间。求出头围关于年龄的回归曲线后，对于给定的年龄，可求出头围关于年龄的回归曲线后，对于给定的年龄，可以调用以调用nlpredci函数函数求出头围的预测值（即求出头围的预测值（即x 给定时给定时y 的总体的总体均值）、预测值的置信区间和观测值的预测区间。均值）、预测值的置信区间和观测值的预测区间。beta,resid,J,Sigma=nlinfit(X,y,fun,b0)ypred,delta=nlpredci(fun,x,beta,resid,covar,Sig

24、ma)或或 ypred,delta=nlpredci(fun,x,beta,resid,jacobian,J)2023-5-57.预测区间图预测区间图 xx=0:0.1:16;%计算给定年龄处头围预测值和预测区间计算给定年龄处头围预测值和预测区间 ypred,delta=nlpredci(HeadCir2,xx,beta,r,covar,.COVB,mse,mse,predopt,observation);yup=ypred+delta;ydown=ypred-delta;figure hold on h1=fill(xx;flipud(xx),yup;flipud(ydown),0.5,0.

25、5,0.5);set(h1,EdgeColor,none,FaceAlpha,0.5);plot(xx,yup,r-,LineWidth,2);plot(xx,ydown,b-.,LineWidth,2);plot(xx,ypred,k,linewidth,2)grid on ylim(32,57)xlabel(年龄年龄(x)ylabel(头围头围(y)h2=legend(预测区间预测区间,预测区间上限预测区间上限,预测区间下限预测区间下限,回归曲线回归曲线);set(h2,Location,SouthEast)2023-5-52023-5-5三、三、利用曲线拟合工具利用曲线拟合工具cftoo

26、l作一元非线性拟合作一元非线性拟合1.cftool函数的用法函数的用法 cftool cftool(xdata,ydata)cftool(xdata,ydata,w)12nxxxx12nyyyy权重向量权重向量2023-5-52.导入数据导入数据 HeadData=xlsread(examp19_2_1.xls);HeadData=sortrows(HeadData,4);x=HeadData(:,4);y=HeadData(:,9);2023-5-53.数据的平滑处理数据的平滑处理2023-5-54.数据筛选数据筛选数据筛选界面数据筛选界面手动选点去除数据界面手动选点去除数据界面2023-5

27、-55.数据拟合数据拟合 导入头围和年龄的导入头围和年龄的数据之后，曲线拟合数据之后，曲线拟合主界面的坐标系里出主界面的坐标系里出现了相应的散点图。现了相应的散点图。单击曲线拟合主界面单击曲线拟合主界面上的上的“Fitting”按钮，按钮，在弹出的数据拟合界在弹出的数据拟合界面上单击面上单击“New fit”按钮，将创建一个新按钮，将创建一个新的拟合的拟合2023-5-56.拟合结果拟合结果 对于前面给出的对于前面给出的1281组头围和年龄的观测数据，至少可组头围和年龄的观测数据，至少可以用以用5种函数进行拟合，得到的非线性回归方程分别为种函数进行拟合，得到的非线性回归方程分别为 负指数曲线：

28、负指数曲线：0.26760.790652.43xye 双曲线：双曲线：0.66440.019070.01779xyx 幂函数曲线：幂函数曲线：0.0590745.22(0.05)yx Logistic曲线：曲线：(4.634)50.31 32.96xye 对数曲线：对数曲线：45.182.859ln(0.05)yx拟合较好拟合较好2023-5-57.拟合效果图拟合效果图2023-5-5第三节第三节 多元线性和广义线性回归多元线性和广义线性回归2023-5-5【例例19.3-1】在有氧锻炼中，人的耗氧能力在有氧锻炼中，人的耗氧能力 y 是衡量身体状况的是衡量身体状况的重要指标，它可能与以下因素有

29、关：年龄重要指标，它可能与以下因素有关：年龄 x1（岁），体重（岁），体重 x2（kg），），1500米跑所用的时间米跑所用的时间 x3（min），静止时心速），静止时心速 x4（次（次/min），跑步后心速），跑步后心速 x5（次（次/min）.对对24名名40至至57岁的志愿者进岁的志愿者进行了测试，结果如表行了测试，结果如表19.3-1所示。表所示。表19.3-1中的数据保存在文件中的数据保存在文件examp19_3_1.xls中，试根据这些数据建立耗氧能力中，试根据这些数据建立耗氧能力y 与诸因素与诸因素之间的回归模型。之间的回归模型。2023-5-5 xydata=xlsread(e

30、xamp19_3_1.xls);y=xydata(:,2);X=xydata(:,3:7);xdata=ones(size(X,1),1),X;b,bint,r,rint,s=regress(y,xdata);一、一、调用调用regress函数作多元线性回归函数作多元线性回归1.理论回归方程理论回归方程01 122334455 ybb xb xb xb xb x2.参数估计及显著性检验参数估计及显著性检验3.经验回归方程经验回归方程12345121.16550.34710.01674.29030.03990.1587yxxxxx2023-5-5 reglm(y,X)二、二、调用自编调用自编re

31、glm函数作多元线性回归函数作多元线性回归 理论回归方程理论回归方程01 122334455 ybb xb xb xb xb x 参数估计及显著性检验参数估计及显著性检验1.五元线性回归五元线性回归-方差分析表方差分析表-方差来源方差来源 自由度自由度 平方和平方和 均方均方 F值值 p值值回归回归 5.0000 625.3110 125.0622 16.0069 0.0000残差残差 18.0000 140.6340 7.8130总计总计 23.0000 765.9450 均方根误差均方根误差(Root MSE)2.7952 判定系数判定系数(R-Square)0.8164 因变量均值因变量

32、均值(Dependent Mean)47.6750 调整的判定系数调整的判定系数(Adj R-Sq)0.7654-参数估计参数估计-变量变量 估计值估计值 标准误标准误 t值值 p值值常数项常数项 121.1655 17.4064 6.9610 0.0000X1 -0.3471 0.1435 -2.4185 0.0264X2 -0.0167 0.0874 -0.1914 0.8504X3 -4.2903 1.0268 -4.1784 0.0006X4 -0.0399 0.0942 -0.4236 0.6769X5 -0.1587 0.0788 -2.0122 0.05942023-5-5 X1

33、35=X(:,1 3 5);varnames=X1,X3,X5;reglm(y,X135,varnames);理论回归方程理论回归方程01 13355 ybb xb xb x 参数估计及显著性检验参数估计及显著性检验2.三元线性回归（剔除不显著项重新回归）三元线性回归（剔除不显著项重新回归）-方差分析表方差分析表-方差来源方差来源 自由度自由度 平方和平方和 均方均方 F值值 p值值回归回归 3.0000 623.7206 207.9069 29.2364 0.0000残差残差 20.0000 142.2244 7.1112总计总计 23.0000 765.9450 均方根误差均方根误差(Ro

34、ot MSE)2.6667 判定系数判定系数(R-Square)0.8143 因变量均值因变量均值(Dependent Mean)47.6750 调整的判定系数调整的判定系数(Adj R-Sq)0.7865-参数估计参数估计-变量变量 估计值估计值 标准误标准误 t值值 p值值常数项常数项 118.0135 14.3399 8.2297 0.0000X1 -0.3254 0.1288 -2.5274 0.0200X3 -4.5694 0.7741 -5.9026 0.0000X5 -0.1561 0.0750 -2.0809 0.05052023-5-5 reglm(y,X,quadratic

35、)三、三、调用自编调用自编reglm函数作二次回归函数作二次回归1.理论回归方程理论回归方程5455201111iiijijiiiiij iiybb xb x xb x 2.参数估计及显著性检验参数估计及显著性检验-方差分析表方差分析表-方差来源方差来源 自由度自由度 平方和平方和 均方均方 F值值 p值值回归回归 20.0000 765.0129 38.2506 123.1051 0.0010残差残差 3.0000 0.9321 0.3107总计总计 23.0000 765.9450 均方根误差均方根误差(Root MSE)0.5574 判定系数判定系数(R-Square)0.9988 因变

36、量均值因变量均值(Dependent Mean)47.6750 调整的判定系数调整的判定系数(Adj R-Sq)0.99072023-5-5四、拟合效果图四、拟合效果图 s1=reglm(y,X);X135=X(:,1 3 5);varnames=X1,X3,X5;s2=reglm(y,X135,varnames);s3=reglm(y,X,quadratic);%绘制拟合效果图绘制拟合效果图 figure;plot(y,ko);hold on plot(s1.yhat,:,linewidth,2);plot(s2.yhat,r-.,linewidth,2);plot(s3.yhat,k,li

37、newidth,2);legend(y的原始散点的原始散点,5元线性回归拟合元线性回归拟合,.3元线性回归拟合元线性回归拟合,完全二次回归拟合完全二次回归拟合);xlabel(y的观测序号的观测序号);%为为X轴加标签轴加标签 ylabel(y);%为为Y轴加标签轴加标签2023-5-5051015202535404550556065y的观测序号y y的原始散点5元线性回归拟合3元线性回归拟合完全二次回归拟合2023-5-5五、五、调用调用stepwise函数作逐步回归函数作逐步回归1.stepwise函数的用法函数的用法111212122212ppnnnpxxxxxxxxx stepwise

38、(X,y,inmodel,penter,premove)初始模型初始模型中所包含中所包含变量的索变量的索引引,p个元个元素的向量素的向量剔除变量剔除变量的最小的最小p，默认为默认为0.112nyyy引入变量引入变量的最大的最大p，默认为默认为0.05 函数运行后出现一交互式界面，通过该界面进行引入和函数运行后出现一交互式界面，通过该界面进行引入和剔除变量的操作，还可以导出相关结果。剔除变量的操作，还可以导出相关结果。2023-5-52.调用调用stepwise函数作逐步回归函数作逐步回归 xydata=xlsread(examp19_3_1.xls);y=xydata(:,2);X=xydat

39、a(:,3:7);inmodel=1:5;stepwise(X,y,inmodel);135118.0130.32544.56940.1561yxxx 初始理论回归方程初始理论回归方程01 122334455 ybb xb xb xb xb x2023-5-52.调用调用stepwise函数作逐步回归（续）函数作逐步回归（续）xydata=xlsread(examp19_3_1.xls);y=xydata(:,2);x1=xydata(:,3);x2=xydata(:,4);x3=xydata(:,5);x4=xydata(:,6);x5=xydata(:,7);X=x1,x2,x3,x4,x

40、5,x1.*x2,x2.*x3,x2.*x5,x4.2,log(x2),sqrt(x5);inmodel=1:size(X,2);stepwise(X,y,inmodel);初始理论回归方程初始理论回归方程520612723825941021151lniiiybb xb x xb x xb x xb xbxbx2023-5-5第四节第四节 多元非线性回归多元非线性回归 地震震中位置的确定地震震中位置的确定2023-5-5【例例19.4-1】2011年年4月月1日某时在某一地点发生了一次地震，图日某时在某一地点发生了一次地震，图19.4-1中中10个地震观测站点均接收到了地震波，观测数据如表个地

41、震观测站点均接收到了地震波，观测数据如表19.4-1所列。假定地震波在各种介质和各个方向的传播速度均相所列。假定地震波在各种介质和各个方向的传播速度均相等，并且在传播过程中保持不变。请根据表等，并且在传播过程中保持不变。请根据表19.4-1中的数据确定中的数据确定这次地震的震中位置、震源深度以及地震发生的时间（不考虑这次地震的震中位置、震源深度以及地震发生的时间（不考虑时区因素，建议时间以分为单位）。时区因素，建议时间以分为单位）。一、问题描述一、问题描述2023-5-5图图19.4-1 地震观测站点示意图地震观测站点示意图05001000150020002500300035000500100

42、015002000250030003500X（单 位：千 米）Y（单 位：千 米）ABCDEFGHIJ2023-5-5表表19.4-1 地震观测站坐标及接收地震波时间地震观测站坐标及接收地震波时间地震观测站地震观测站横坐标横坐标x（千米）（千米）纵坐标纵坐标y（千米）（千米）接收地震波时间接收地震波时间A50033004月1日9时21分9秒B3002004月1日9时19分29秒C80016004月1日9时14分51秒D140022004月1日9时13分17秒E17007004月1日9时11分46秒F230028004月1日9时14分47秒G250019004月1日9时10分14秒H290090

43、04月1日9时11分46秒I320031004月1日9时17分57秒J34001004月1日9时16分49秒2023-5-5二、模型建立二、模型建立22200000()(),1,2,1060iiiixxyyzTtiv其中其中 为随机误差，为随机误差，为模型参数。为模型参数。i00000,xyzvt2023-5-5%理论回归方程所对应的理论回归方程所对应的M函数（需要在程序编辑窗口编写）函数（需要在程序编辑窗口编写）function T=modelfun(beta,xy)x=xy(:,1);%x坐标坐标y=xy(:,2);%y坐标坐标T=sqrt(x-beta(1).2+.(y-beta(2).

44、2+beta(3).2)/(60*beta(4)+beta(5);1.理论回归方程所对应的理论回归方程所对应的M函数函数三、调用三、调用nlinfit函数作多元非线性回归函数作多元非线性回归2023-5-5%定义地震观测站位置坐标及接收地震波时间数据矩阵定义地震观测站位置坐标及接收地震波时间数据矩阵x,y,Minutes,Secondsxyt=500 3300 21 9 300 200 19 29 800 1600 14 51 1400 2200 13 17 1700 700 11 46 2300 2800 14 47 2500 1900 10 14 2900 900 11 46 3200 3

45、100 17 57 3400 100 16 49;xy=xyt(:,1:2);Minutes=xyt(:,3);Seconds=xyt(:,4);T=Minutes+Seconds/60;beta0=1000 100 1 1 1;b=nlinfit(xy,T,modelfun,beta0)b=1.0e+003*2.2005 1.3999 0.0351 0.0030 0.00702.参数估计参数估计2023-5-53.结果结果000002200.51399.935.137xyzvt 也就是说地震发生的时间为也就是说地震发生的时间为2011年年4月月1日日09时时07分，震分，震中位于中位于 处，

46、震源深度处，震源深度35.1公里。公里。002200.5,1399.9xy2023-5-5第五节第五节 补充内容补充内容 数据拟合的其他数据拟合的其他MATLAB函数函数2023-5-5一、多项式回归函数：一、多项式回归函数：polyfit，polytool12.nxxXx1.确定回归系数的点估计值确定回归系数的点估计值 p,s =polyfit (X,Y,m )1121.mmmmya xa xa xanYYYY.21121.maapa次数次数2023-5-52.预测和预测误差估计预测和预测误差估计 Yhat =polyval (p,xdat )预测值预测值系数估计值系数估计值指定指定x 值值

47、Yhat,DELTA =polyconf(p,xdat,s,alpha)(1)Yhat=polyval(p,xdat)求求polyfit所得的回归多项式在所得的回归多项式在xdat 处的预测值处的预测值Y；(2)Yhat,DELTA=polyconf(p,xdat,s,alpha)求求polyfit所得的回归多项式在所得的回归多项式在x 处的预测值处的预测值Yhat及预测及预测值的置信水平为值的置信水平为1-alpha的的置信区间置信区间Yhat DELTA；alpha缺省时为缺省时为0.05.2023-5-5【例例19.5-1】出钢时所用的盛钢水的钢包，由于钢水对耐火出钢时所用的盛钢水的钢包

48、，由于钢水对耐火材料的浸蚀，容积不断增大。我们希望找到使用次数与增大材料的浸蚀，容积不断增大。我们希望找到使用次数与增大的容积之间的关系。对一钢包做试验，测得数据列于下表：的容积之间的关系。对一钢包做试验，测得数据列于下表：（1）作出散点图；）作出散点图；（2）求）求y 关于关于x 的经验回归方程；的经验回归方程；2023-5-5q 原始数据散点与折线图原始数据散点与折线图24681012141666.577.588.599.51010.511使 用 次 数增大容积2023-5-5xy=2 3 4 5 6 7 8 9 10 11 12 13 14 15 166.42 8.2 9.58 9.5

49、9.7 10 9.93 9.99 10.49 10.59 10.6 10.8 10.6 10.9 10.76;x=xy(1,:);y=xy(2,:);figure(1)plot(x,y,bo);grid on;xlabel(使用次数使用次数);ylabel(增大容积增大容积)set(gca,color,none)p,s=polyfit(x,y,2);yhat,delta=polyconf(p,x,s);py yhat y-yhat yhat-delta yhat+deltafigure(2)plot(x,y,bo,x,yhat,r,x,yhat-delta,c,x,yhat+delta,c);

50、grid on;xlabel(使用次数使用次数);ylabel(增大容积增大容积)set(gca,color,none)q 调用调用polyfit函数作多项式拟合的函数作多项式拟合的Matlab程序程序12023-5-5 p=-0.029 0.7408 6.0927q 部分结果部分结果Y YYhatYhatr rYhat-deltaYhat-deltaYhat+deltaYhat+delta6.426.427.45847.4584-1.0384-1.03846.19216.19218.72468.72468.28.28.05428.05420.14580.14586.8766.8769.232