回归分析MATLAB实现

1/24MATLAB数据分析方法（机械工业出版社）第3章回归分析普通高等院校计算机课程规划教材MATLAB数据分析方法李柏年吴礼斌主编张孔生丁华参编2/24MATLAB数据分析方法（机械工业出版社）第3章回归分析回归分析是最常用的数据分析方法之一。它是根据已得的试验结果以及以往的经验来建立统计模型，并研究变量间的相关关系，建立起变量之间关系的近似表达式即经验公式，并由此对相应的变量进行预测和控制等.3.1一元回归模型3.1.1一元线性回归模型1.一元线性回归的基本概念3/24MATLAB数据分析方法（机械工业出版社）第3章回归分析通常，我们对总体(x,Y)进行n次的独立观测，获得n组数据（称为样本观测值）(x1,y1),(x2,y2),…,(xn,yn)利用最小二乘法可以得到回归模型参数0，1的最小二乘估计设Y是一个可观测的随机变量，它受到一个非随机变量因素x和随机误差的影响。若Y与x有如下线性关系：,10xY(3.1.1)且E=0,D=2,则称(3.1.1)为一元线性回归模型.其中0，1为回归系数，x为自变量，Y为因变量.10ˆ,ˆ4/24MATLAB数据分析方法（机械工业出版社）第3章回归分析(3.1.2)其中于是建立经验公式模型：(3.1.3)一元线性回归分析的主要任务：一是利用样本观测值对回归系数0，1和作点估计；二是对方程的线性关系即1作显著性检验；三是在x=x0处对Y作预测等.以下举例说明建立经验公式（3.1.3）的方法。.ˆ,ˆˆ110xxxyLLxy，niixnx11，niiyny11，niixxxxL12)()()(1yyxxLiniixy01ˆˆˆyx5/24MATLAB数据分析方法（机械工业出版社）第3章回归分析例3.1.1近10年来，某市社会商品零售总额与职工工资总额（单位：亿元）数据如下表3.1。表3.1商品零售总额与职工工资表（单位：亿元）建立社会商品零售总额与职工工资总额数据的回归模型工资总额23.827.631.632.433.734.943.252.863.873.4零售总额41.451.861.767.968.777.595.9137.4155.0175.0解：%首先输入数据x=[23.80,27.60,31.60,32.40,33.70,34.90,43.20,52.80,63.80,73.40];y=[41.4,51.8,61.70,67.90,68.70,77.50,95.90,137.40,155.0,175.0];6/24MATLAB数据分析方法（机械工业出版社）第3章回归分析%然后作散点图plot(x,y,'*')%作散点图xlabel('x(职工工资总额)')%横坐标名ylabel('y(商品零售总额)')%纵坐标名20304050607080406080100120140160180x(职工工资总额)y(商品零售总额)图3.1商品零售总额与职工工资总额数据散点图7/24MATLAB数据分析方法（机械工业出版社）第3章回归分析%计算最佳参数Lxx=sum((x-mean(x)).^2);Lxy=sum((x-mean(x)).*(y-mean(y)));b1=Lxy/Lxx;b0=mean(y)-b1*mean(x);运行后得到：b1=2.7991，b0=-23.5493所以，回归模型为5493.237991.2ˆxy问题1：当x=0，得到y=-23.5493亿元如何理解？问题2：如何检验E=0?D=2?8/24MATLAB数据分析方法（机械工业出版社）第3章回归分析2.一元多项式回归模型1110...nnnnyaxaxaxa在一元回归模型中，如果变量y与x的关系是n次多项式，即其中是随机误差，服从正态分布N(0,2)a0,a1,…,an为回归系数，则称(3.1.4)为多项式回归模型.(3.1.4)(1)多项式曲线拟合在MATLAB7的统计工具箱中，有多项式曲线拟合的命令polyfit，其调用格式有以下三种：9/24MATLAB数据分析方法（机械工业出版社）第3章回归分析p=polyfit(x,y,n)[p,S]=polyfit(x,y,n)[p,S,mu]=polyfit(x,y,n)其中，输入x,y分别为自变量与因变量的样本观测数据向量；n是多项式的阶数，对于一元线性回归则取n=1;输出p是按照降幂排列的多项式的系数向量，S是一个矩阵，用于估计预测误差或供MATLAB的其它函数的调用。例3.1.2某种合金中的主要成分为A,B两种金属，经过试验发现：这两种金属成分之和x与合金的膨胀系数y有如下关系，建立描述这种关系的数学表达式.10/24MATLAB数据分析方法（机械工业出版社）第3章回归分析表3.2合金的膨胀系数表解：%首先输入数据x=37:0.5:43;y=[3.4,3,3,2.27,2.1,1.83,1.53,1.7,1.8,1.9,2.35,2.54,2.9];%其次做散点图plot(x,y,‘*’)xlabel('x(两种合金之和)')%横坐标名ylabel(‘y(合金膨胀系数)’)%纵坐标名%然后根据散点图猜测曲线类别（2.1.7）x3737.53838.53939.54040.54141.54242.543y3.4332.272.11.831.531.71.81.92.352.542.911/24MATLAB数据分析方法（机械工业出版社）第3章回归分析由于散点图呈抛物线，故选择二次函数曲线进行拟合.p=polyfit(x,y,2)%注意取n=2运行得到回归系数：p=0.1660-13.3866271.6231即二次回归模型为：2ˆ0.16613.3866271.6231yxx12/24MATLAB数据分析方法（机械工业出版社）第3章回归分析多项式曲线拟合预测的命令polyval，其调用格式有以下两种：Y=polyval(p,x0)[Y,Delta]=polyconf(p,x0,S,alpha)其中，输入p,S是由多项式拟合命[p,S]=polyfit(x,y,n)的输出，x0是要预测的自变量的值.输出Y是polyfit所得的回归多项式在x处的预测值。(2)多项式回归的预测与置信区间如果输入数据的误差相互独立，且方差为常数，则Y±Delta至少包含95%的预测值；alpha缺省时为0.05。(Y-Delta,Y+Delta)即95%的置信区间13/24MATLAB数据分析方法（机械工业出版社）第3章回归分析(3)多项式回归的GUI界面命令多项式回归的GUI界面命令polytool，其典型调用格式polytool(x,y,n,alpha)其中，输入x,y分别为自变量与因变量的样本观测数据向量；n是多项式的阶数；置信度为(1-alpha)%，alpha缺省时为0.05。该命令可以绘出总体拟合图形以及(1-alpha)上、下置信区间的直线（屏幕上显示为红色）.此外，用鼠标拖动图中纵向虚线，就可以显示出对于不同的自变量数值所对应的预测状况，与此同时图形左端数值框中会随着自变量的变化而得到的预报数值以及(1-alpha)置信区间长度一半的数值。14/24MATLAB数据分析方法（机械工业出版社）第3章回归分析15/24MATLAB数据分析方法（机械工业出版社）第3章回归分析例3.1.3为了分析X射线的杀菌作用，用200千伏的X射线来照射细菌，每次照射6分钟用平板计数法估计尚存活的细菌数，照射次数记为t，照射后的细菌数y如表3.3所示。t123456789101112131415y3522111971601421061046056383632211915表3.3X射线照射次数与残留细菌数试求：①给出y与t的二次函数回归模型；②在同一坐标系内做出原始数据与拟合结果的散点图③预测t=16时残留的细菌数；④根据问题实际意义选择多项式函数是否合适？数据来源：http//~hadi/RABE16/24MATLAB数据分析方法（机械工业出版社）第3章回归分析解：%输入原始数据t=1:15;y=[352,211,197,160,142,106,104,60,56,38,36,32,21,19,15];p=polyfit(t,y,2);%作二次多项式回归y1=polyval(p,t);%模型估计与作图plot(t,y,'-*',t,y1,'-o');legend('原始数据','二次函数')xlabel('t(照射次数)')ylabel('y(残留细菌数)')t0=16;yc1=polyconf(p,t0)%预测t0=16时残留的细菌数运行结果为p=1.9897-51.1394347.8967，yc1=39.0396即二次回归模型为211.9897t-51.1394t+347.8967y17/24MATLAB数据分析方法（机械工业出版社）第3章回归分析yc1=39.0396,表明照射16次后，用二次函数计算出细菌残留数为39.0396，显然与实际不相符合。调用多项式回归的GUI界面命令polytool,如图3.4原始数据与拟合结果的散点图如图3.3所示，从图形可知拟合效果较好.图3.3原始数据与拟合结果的散点图18/24MATLAB数据分析方法（机械工业出版社）第3章回归分析根据实际问题的意义可知：尽管二次多项式拟合效果较好，但是用于预测并不理想。因此如何根据原始数据散点图的规律，选择适当的回归曲线是非常重要的，因此有必要研究非线性回归分析.图3.4二次函数预测交互图19/24MATLAB数据分析方法（机械工业出版社）第3章回归分析3.1.2一元非线性回归模型为了便于正确地选择合适的函数进行回归分析建模，我们给出通常选择的六类曲线如下所示：1.非线性曲线选择（1）双曲线1/y=a+b/x(见图3.5)。图3.5双曲线图3.5双曲线20/24MATLAB数据分析方法（机械工业出版社）第3章回归分析(2)幂函数曲线y=axb,其中x0,a0(图3.6)。图3.6幂函数曲线（3）指数曲线y=aebx，其中参数a0(见图3.7)。图3.7指数曲线21/24MATLAB数据分析方法（机械工业出版社）第3章回归分析xbaey/（4）倒指数曲线,其中a0(图3.8)。图3.8倒指数曲线（5）y=a+blnx(见图3.9)。图3.9对数曲线22/24MATLAB数据分析方法（机械工业出版社）第3章回归分析xbeay1（6）S型曲线(见图3.10)。图3.10S型曲线对于非线性回归建模通常有两种方法：一是通过适当的变换转化为线性回归模型，例如双曲线模型(图3.5)。如果无法实现线性化，可以利用最小二乘法直接建立非线性回归模型，求解最佳参数。23/24MATLAB数据分析方法（机械工业出版社）第3章回归分析2.非线性回归的MATLAB命令MATLAB统计工具箱中实现非线性回归的命令有nlinfit、nlparci、lpredci和nlintool。下面逐一介绍调用格式。非线性拟合命令nlinfit，调用格式：[beta,r,J]=nlinfit(x,y,'model',beta0)其中，输人数据x，y分别为n×m矩阵和n维列向量，对一元非线性回归，x为n维列向量，model是事先用M文件定义的非线性函数，beta0是回归系数的初值(需要通过解方程组得到)，beta是估计出的最佳回归系数，r是残差，J是Jacobian矩阵，它们是估计预测误差需要的数据。24/24MATLAB数据分析方法（机械工业出版社）第3章回归分析非线性回归预测命令nlpredci，调用格式：ypred=nlpredci(FUN,inputs,beta,r,J)其中，输入参数beta,r,J是非线性回归命令nlinfit的输出结果,FUN是拟合函数，inputs是需要预测的自变量；输出量ypred是inputs的预测值。非线性回归置