Matlab解决线性回归问题论文

Matlab解决线性回归问题2010级电子信息工程班【摘要】MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指°MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。回归分析，是对现有数据进行处理、从中发现有用信息的一种重要手段。而线性回归，特别是一元线性回归分析更是人们优先考虑采用的方式。基于此，本文就一元线性回归的MATLAB实现作了一番探讨，给出了多种实现方式，并通过一个实例加以具体展示，在数据处理时可根据自己的需要灵活地加以选用。【关键词】MATLAB程序、一元线性回归、多元线性回归问题、线性回归拟合问题、Subplot、三次样条插值函数。一、 提出问题在回归分析中，如果有两个或两个以上的自变量，就称为多元回归。事实上，一种现象常常是与多个因素相联系的，由多个自变量的最优组合共同来预测或估计因变量，比只用一个自变量进行预测或估计更有效，更符合实际。二、 简述原理（一）一元线性回归命令polyfit最小二乘多项式拟合[p，S]=polyfit（x，y，m）多项式y二a1xm+a2xm—1+…+amx+am+1其中x二（xl,x2,…，xm）xl・・・xm为（n*l）的矩阵；y为（n*l）的矩阵；p=（al,a2,…，am+1）是多项式y二a1xm+a2xm-1+・・・+amx+am+1的系数；S是一个矩阵，用来估计预测误差.命令polyval多项式函数的预测值Y二polyval（p,x）求polyfit所得的回归多项式在x处的预测值Y；p是polyfit函数的返回值；x和polyfit函数的x值相同。命令polyconf残差个案次序图[Y,DELTA]二polyconf（p,x,S,alpha）求polyfit所得的回归多项式在x处的预测值Y及预测值的显著性为1-alpha的置信区间DELTA；alpha缺省时为0.05。p是polyfit函数的返回值；x和polyfit函数的x值相同；S和polyfit函数的S值相同。命令polytool（x,y,m）—元多项式回归命令命令regress多元线性回归（可用于一元线性回归）b=regress（Y,X）[b,bint,r,rint,stats]=regress（Y,X,alpha）b回归系数bint回归系数的区间估计r残差rint残差置信区间stats用于检验回归模型的统计量，有三个数值：相关系数R2、F值、与F对应的概率p,相关系数R2越接近1,说明回归方程越显著；F>F1-a（k,n-k-1）时拒绝HO,F越大，说明回归方程越显著；与F对应的概率p时拒绝HO,回归模型成立。Y为n*1的矩阵；乂为（ones（n,1）,x1,…，xm）的矩阵；alpha显著性水平（缺省时为0.05）。（二）多元线性回归命令regress命令rstool多元二项式回归命令：rstool（x,y,'model',alpha）x为n*m矩阵y为n维列向量model由下列4个模型中选择1个（用字符串输入,缺省时为线性模型）：linear（线性）：purequadratic（纯二次）：interaction（交叉）：quadratic（完全二次）：alpha显著性水平(缺省时为0.05)返回值beta系数返回值rmse剩余标准差返回值residuals残差非线性回归1．命令nlinfit[beta,R,J]=nlinfit(X,Y,''model',beta0)X为n*m矩阵Y为n维列向量model为自定义函数betaO为估计的模型系数beta为回归系数R为残差2．命令nlintoolnlintool(X,Y,'model',betaO,alpha)X为n*m矩阵Y为n维列向量model为自定义函数beta0为估计的模型系数alpha显著性水平(缺省时为0.05)3．命令nlparcibetaci=nlparci(beta,R,J)beta为回归系数R为残差返回值为回归系数beta的置信区间4．命令nlpredci[Y,DELTA]=nlpredci(‘model',X,beta,R,J)Y为预测值DELTA为预测值的显著性为1-alpha的置信区间；alpha缺省时为0.05。X为n*m矩阵model为自定义函数beta为回归系数R为残差三、多元线性回归问题在实际经济问题中，一个变量往往受到多个变量的影响。例如，家庭消费支出，除了受家庭可支配收入的影响外，还受诸如家庭所有的财富、物价水平、金融机构存款利息等多种因素的影响，表现在线性回归模型中的解释变量有多个。这样的模型被称为多元线性回归模型。(multivariablelinearregressionmodel)1、多元线性回归模型的一般形式多元线性回归模型的一般形式为：

TOC\o"1-5"\h\zY=0+lXli+2X2i++kXki+i,i=l,2,,n (1)其中k为解释变量的数目，j(jl,2,k)称为回归系(regressioncoefficient)。式也被称为总体回归函数的随机表达式。它的非随机表达式为：Y=0+1X1i+2X2i++kXki,i=1,2,n (2)j也被称为偏回归系数(partialregressioncoefficient)。2、多元线性回归计算模型Y=0+lxl+2x2+kxk+N(0,2) (3)多元性回归模型的参数估计，同一元线性回归方程一样，也是在要求误差平方和(Ze)为最小的前提下，用最小二乘法或最大似然估计法求解参数。设(xll, x12 , x1 p, y1 ),,( xn1, xn2,xnp,yn)是一个样本，用最大似然估计法估计参数：取b0,blbp，当bObO,b1b1,bb达到最小。例1输入以下程序:%数据y=[8SJ85.08,91,92.93,93,95.96,98,97.96,98,99,100,102]:3=[onesIlength(y)jl)js']%数据[b,bintjr：rint,stata]=reEre53(Y,X); %制图输出结果如图1所示:

图1图1在Matlab命令窗口输入z=b(l)+b(2)*x%运算%输出运算plot(x,Y,‘k+，,x,z,，r，)%输出运算得到预测比较图如图2、3所示:

例2：（多元线性回归）水泥凝固时放出的热量y与水泥中的四种化学成分x1,x2,x3,x4有关,今测得一组数据如下，试确定多元线性模型

序号卢和；冠8'、fip初；11^'P；1V'愛：也n2闘2沖亦5V71Px3^6^15P8r鬧存阶込47V：33p6^7«,5^74.3^104.3^87.6^953109.2^1睦.斑序号户扣伽他12*^T口公「21^侗11+J伽3佃珈4衍斗W盹x3^耳1腳4+1湖■:22卩26^珈12^72.5P93.143115.^83.8^109.4+1□分折：诅U413,1坯4,11J0J4J竝=[2伉29,托」1.52,557131,54,47,40^,15盯汁^=[615^8 1了卫,1S,4罟兌町屮材=[150,52^,4733^2,6,44^2益』4,12,12]屮尸卩8.5,743,104.3,87.6,95.9,109：2,102.77^->,93.1,115.9,8^.8,113.3,109.4]；^由武?应可得X=[@TRT,J«T,0T护T](旳为单位列向星1屮Y=jT..Matlab程序为;:，:*■'d输入勿]區命令：3小卩m?；厂门£ %数据必呼护㈠“沐壬门“沖尹汽％数据J、：、•，•:.-■■-■二：■..■；'■ %数据■--■■：.二-'二■■-、…:-■=■- %数据y=[78.5.74.3,104.3.87.6.95.9J0^2.10^.7.72.5.93.U1S.9.&3.&.113.3.1014]>%数据X=[ones(length(y)J):x1,.x2\xJ\x4l：%把行向蛍转铁为列向鱼屮丫二丫：％把行向HfeO列向量屮[tirtiintnr:rintzstatsj=regress“肉屮 %运算公式b:bint;stats+J在Matlab图示所示：心'. .'. ' ％数据X率[施29.55..31>52,55,71.丸51,47,迪66,68];灶&尽&血&9,17,22,也嗨网;x4=[60,52,绚卸33,22禺咼22,绳轴12,121.注厲尽他3谢餌肌缶恫嗾叽型疏人隆氐0$】昇遽瘀閱禺1注迢呃4];X=[onesQwifthfyh1）注1’注2*,讶注4’];%数据%数据%数据%运行方程输出结果如图4所示:=62.40S4:p 的置信区间(t99.1786,2235893)b-.=1.5511^£的置信区间(10.1663,3.26E5)/因此我们可得矶=05102」bnz的置信区间(11.1589,2.1792)^^=0.1019>阴的置信区间<n.63S5:1.S423)鋭=口441彳获的置信区间017791：1.4910>口=0朋24』=111.4792^=0.0000.^p<W回归模型」162.4054+1.5511x1+0.5102x2+0.1019^3-0.1441^立一结果如图5所示:

图5在matlab命令框中输入卫二feCHb僅门1+b(3)*x2+b(4}^3+t)译片4屮plotfX.Y.-k+rXz.T>得到如图6所示:四、线性回归拟合对于多元线性回归模型y=B+BxH Px+e011pp设变量Xi,xj…Xp,y盼组观测值为(x,x，•••兀,y)i=1,2,•…,n・i1 i2 ipi、1

记x=••、1

记x=••1x…xry)rp)121p10x22•••…x2p••••••,y=y2•，则p=p1•xn2…xnp<y丿njp;x11x21•••xn1的估计值为b=P=(x'x)-1x'y(11.2)在Matlab中，用regress函数进行多元线性回归分析，应用方法如下:语法：b=regress(y,x)[b,bint,r,rint,stats]=regress(y,x)[b,bint,r,rint,stats]=regress(y,x,alpha)b=regress(y,x)得到的p+1维列向量b即为(11.2)式给出的回归系数P的估计值・[b,bint,r,rint,stats]=regress(y,x)给出回归系数P的估计值b，p的95%置信区间((p+1)x2向量)bint,残差r以及每个残差的95%置信区间(nx2向量)rint；向量stats给出回归的R2统计量和F以及临界概率p的值.如果卩的置信区间(bint的第i+1行)不包含0,则在显著水平为a时拒绝ip=0的假设，认为变量x是显著的.ii[b,bint,r,rint,stats]=regress(y,x,alpha)给出了bint和rint的100(l-alpha)%的置信区间.三次样条插值函数的MATLAB程序matlab的splinex=0:10;y=sin(x); %插值点xx=0:.25:10; %绘图点yy=spline(x,y,xx);plot(x,y,'o',xx,yy)结果如图7所示图7五、非线性拟合非线性拟合可以用以下命令（同样适用于线形回归分析）：beta=nlinfit（X,y,fun,betaO）X给定的自变量数据,Y给定的因变量数据,fun要拟合的函数模型（句柄函数或者内联函数形式），beta0函数模型中系数估计初值,beta返回拟合后的系数x=lsqcurvefit（fun,x0,xdata,ydata）fun要拟合的目标函数,x0目标函数中的系数估计初值,xdata自变量数据,ydata函数值数据X拟合返回的系数（拟合结果）nlinfit格式：[beta,r，J]=nlinfit（x,y,'model',beta0）Beta估计出的回归系数r残差JJacobian矩阵x，y输入数据x、y分别为n*m矩阵和n维列向量，对一元非线性回归，x为n维列向量。model是事先用m-文件定义的非线性函数beta0回归系数的初值例1已知数据：xl二[0.5,0.4,0.3,0.2,0.1];x2二[0.3,0.5,0.2,0.4,0.6];x3=[1.8,1.4,1.0,1.4,1.8];y二[0.785,0.703,0.583,0.571,0.126]';且y与x1，x2,x3关系为多元非线性关系（只与x2,x3相关）为：y=a+b*x2+c*x3+d*（x2「2）+e*（x3「2）求非线性回归系数a,b,c,d,e。对回归模型建立M文件model.m如下：functionyy二myfun(beta,x)x1=x(:,1);%数据x2=x(:,2);%数据x3=x(:,3);%数据yy=beta(1)+beta(2)*x2+beta(3)*x3+beta(4)*(x2."2)+beta(5)*(x3."2);%运算出图主程序如下：x二[0.5,0.4,0.3,0.2,0.1;0.3,0.5,0.2,0.4,0.6;1.8,1.4,1.0,1.4,1.8]; %数据y=[0.785,0.703,0.583,0.571,0.126];%数据beta0=[1,1,1,1,1];%数据[beta,r,j]=nlinfit(x,y,@myfun,betaO)%运算出图输出结果如图8所示图8例题2：混凝土的抗压强度随养护时间的延长而增加，现将一批混凝土作成12个试块，记录了养护日期(日)及抗压强度y(kg/cm2)的数据：养护时间：x=[234579121417212856]抗压强度：y=[35+r42+r47+r53+r59+r65+r68+r73+r76+r82+r86+r99+r]建立非线性回归模型，对得到的模型和系数进行检验。注明：此题中的+r代表加上一个［-0.5,0.5］之间的随机数模型为：y=a+k1*exp(m*x)+k2*exp(-m*x);Matlab程序：x=[234579121417212856];%数据r二rand(l,12)-0.5;y1=[354247535965687376828699];%数据y=y1+r;%运算路径myfunc二inline('beta(1)+beta(2)*exp(beta(4)*x)+beta(3)*exp(-beta(4)*x)','beta','x');%运算路径beta=nlinfit(x,y,myfunc,[0.50.50.50.5]);%运算方法a二beta(1),k1=beta(2),k2二beta(3),m=beta(4)%数据%testthemodel%运算方法xx=min(x):max(x);%运算方法yy=a+k1*exp(m*xx)+k2*exp(—m*xx);%运算方法plot(x,y,'o',xx,yy,'r')结果如图9所示:图9lsqnonlin非线性最小二乘(非线性数据拟合)的标准形式为minf(x)=f］(x)2+f2(x)2++fm(x)2+Lx其中：L为常数在MATLAB5.X中，用函数leastsq解决这类问题，在6.0版中使用函数lsqnonlin。

设F(x)设F(x)=f1(x)f2(x)fm(x)m则目标函数可表达为min1||F(x)||2=1工f(x)2x2 2 2ii其中：X为向量，F(x)为函数向量。函数lsqnonlin格式x=lsqnonlin(fun,xO)%xO为初始解向量；fun为fi(x),i=l,2,…，m,fun返回向量值F,而不是平方和值，平方和隐含在算法中，fun的定义与前面相同。x=lsqnonlin(fun,xO,lb,ub)%lb、ub定义x的下界和上界：lb<x<ub。x=lsqnonlin(fun,xO,lb,ub,options)%options为扌旨定优化参数，若x没有界，则lb=[]，ub=[]。[x,resnorm]=lsqnonlin(…)％resnorm二sum(fun(x)."2)，即解x处目标函数值。[x,resnorm,residual]=lsqnonlin(…)％residual二fun(x)，即解x处fun的值。[x,resnorm,residual,exitflag]=lsqnonlin(…)％exitflag为终止迭代条件。[x,resnorm,residual,exitflag,output]二lsqnonlin(…)％output输出优化信息。[x,resnorm,residual,exitflag,output,lambda]=lsqnonlin(…)%lambda为Lagrage乘子。[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(…)%fun在解x处的Jacobian矩阵。例求下面非线性最小二乘问题f(2+2k-ekx1-ekx2)2初始解向量为k=1x0=[0.3,0.4]。解：先建立函数文件，并保存为myfun.m，由于lsqnonlin中的fun为向量形式而不是平方和形式，因此，myfun函数应由fi(x)建立：fk(x)=2+2k一ekx1一ekx2 k=l,2,…，10functionF=myfun(x)k=1:10;F=2+2*k-exp(k*x(1))-exp(k*x(2));然后调用优化程序：x0=[0.30.4]；[x,resnorm]=lsqnonlin(@myfun,x0)%主运算方法结果为：Optimizationterminatedsuccessfully:NormofthecurrentstepislessthanOPTIONS.TolXx=0.25780.2578resnorm=%求目标函数值lsqcurvefit非线性曲线拟合是已知输入向量xdata和输出向量ydata,并且知道输入与输出的函数关系为ydata二F(x,xdata)，但不知道系数向量x。今进行曲线拟合,求x使得下式成立：min

x(F(x,xdatai)—ydatai)2min

x(F(x,xdatai)—ydatai)2i在MATLAB5.X中，使用函数curvefit解决这类问题。函数lsqcurvefit格式x=lsqcurvefit(fun,x0,xdata,ydata)%运算方法x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options)[x,resnorm]=lsqcurvefit(…)％运算方法[x,resnorm,residual]=lsqcurvefit(…)％运算方法[x,resnorm,residual,exitflag]=lsqcurvefit(…)％运算方法[x,resnorm,residual,exitflag,output]二lsqcurvefit(…)[x,resnorm,residual,exitflag,output,lambda]=lsqcurvefit(…)%运算方法[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit(…)％运算方法参数说明：xO为初始解向量；xdata,ydata为满足关系ydata二F(x,xdata)的数据;lb、ub为解向量的下界和上界lb<x<ub，若没有指定界，则lb=[]，ub=[]；options为指定的优化参数；fun为拟合函数，其定义方式为：x=lsqcurvefit(@myfun,xO,xdata,ydata),其中myfun已定义为 functionF二myfun(x,xdata)F=…％计算x处拟合函数值fun的用法与前面相同；resnorm二sum((fun(x,xdata)-ydata)."2),即在x处残差的平方和；residual二fun(x,xdata)-ydata,即在x处的残差；exitflag为终止迭代的条件；output为输出的优化信息；lambda为解x处的Lagrange乘子；jacobian为解x处拟合函数fun的jacobian矩阵。例求解如下最小二乘非线性拟合问题已知输入向量xdata和输出向量ydata,且长度都是n,拟合函数为ydata(i)=x(1)-xdata(i)2+x(2)-sin(xdata(i))+x(3)-xdata(i)3即目标函数为 xdal)—ydata.)2x 2i=1其中：F(x,xdata)=x(1)-xdata2+x(2)-sin(xdata)+x(3)-xdata3初始解向量为x0=[0.3,0.4,0.1]。解：先建立拟合函数文件，并保存为myfun.mfunctionF=myfun(x,xdata)

F=x(l)*xdata.2+x(2)*sin(xdata)+x(3)*xdata.3;然后给出数据xdata和ydata>>xdata=[3.67.79.34.18.62.81.37.910.05.4];>>ydata=[16.5150.6263.124.7208.59.92.7163.9325.054.3];>>x0=[10,10,10]; %初始估计值>>[x,resnorm]=lsqcurvefit(@myfun,x0,xdata,ydata)结果为：Optimizationterminatedsuccessfully:RelativefunctionvaluechangingbylessthanOPTIONS.TolFunx=0.2269 0.3385 0.3021resnorm=6.2950对于上述这些可化为线性模型的回归问题，一般先将其化为线性模型，然后再用最小二乘法求出参数的估计值，最后再经过适当的变换，得到所求回归曲线。在熟练掌握最小二乘法的情况下，解决上述问题的关键是确定曲线类型和怎样将其转化为线性模型。如图10所示图10亡3图10亡3整握聲SBW确定曲线类型一般从两个方面考虑：一是根据专业知识，从理论上推导或凭经验推测、二是在专业知识无能为力的情况下，通过绘制和观测散点图确定曲线大体类型。六、Subplot问题程序：x1=[1 1 1 1 1 2 2 2 2 3 333 4 4 4 4 5 5 5 6 6 6 6 7 8 8 8 8 1010 101011 1112 1213 1314 1516 1616 1720]';%数据x2=[1 0 1 0 0 1 0 0 0 0 111 0 1 0 0 0 0 1 0 1 0 1 1 0 1 1 0 001 11001 0 1 0 1 1 0 0 0 0]';%数据x5二[1,3,3,2,3,2,2,1,3,2,1,2,3,1,3,3,2,2,3,1,1,3,2,2,1,2,1,3,1,1,2,3,2,2,1,2,3,1,2,2,3,2,2,1,2,1]';%数据x3=[1000000100100100000110001010110000100100000101]';%数据x4=[0001 011 001010000 110 000 110 100001011010011 011 010]';%数据y=[13876116081870111283117672087211772 10535 121951231314975213711980011417 20263 132311288413245 13677 15965 12366 2135213839228841697814803 17404 221841354814467 15942 23174 23780 2541014861168822417015990 26330 179492568527837 18838 17483 19207 19346];%数据size(x5)x=[ones(46,1),x1,x2,x3,x4][b,bint,r,rint,stats]二regress(y,x,0.05)%运算方法plot(x1,r,'r+')holdonfori=1:1:46if(x2(i,1)==0&x5(i,1)==1)k(i)=1;%交替运算elseif(x2(i,1)==1&x5(i,1)==1)k(i)=2;%交替运算elseif(x2(i,1)==0&x5(i,1)==2)k(i)=3;%交替运算elseif(x2(i,1)==1&x5(i,1)==2)k(i)=4;%交替运算elseif(x2(i,1)==0&x5(i,1)==3)k(i)=5;%交替运算elsek(i)=6;endendk;plot(k,r,'g+')输出如图11所示图11第二个TOC\o"1-5"\h\zx1=[11 1 1 1 2 2 2 2 3 333 4 4 4 4 5 5 5 6 6 6 6 7 8 8 8 8 1010 1010111112 1213 1314 1516 16 16 1720]';%数据x2=[10 1 0 0 1 0 0 0 0 111 0 1 0 0 0 0 1 0 1 0 1 1 0 1 1 0 001 110010 1 0 1 1 0 0 0 0]';%数据x5二[1,3,3,2,3,2,2,1,3,2,1,2,3,1,3,3,2,2,3,1,1,3,2,2,1,2,1,3,1,1,2,3,2,2,1,2,3,1,2,2,3,2,2,1,2,1]';%数据x3=[10 0 0 0 0 0 1 0 01001 0 0 0 0 0 1 1 0 0 0 1 0 1 0 1 1 0 0 0 01001 0 0 0 0 0 1 0 1]';%数据x4=[00 0 1 0 1 1 0 0 10100 0 0 1 1 0 0 0 0 1 1 0 1 0 0 0 0 1 0 1 10100 1 1 0 1 1 0 1 0]';%数据y二[1387611608187011128311