Matlab解决线性回归问题论文.docx

1、Matlab解决线性回归问题论文Matlab解决线性回归问题2010级电子信息工程班【摘要】MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。回归分析，是对现有数据进行处理、从中发现有用信息的一种重要手段。而线性回归，特别是一元线性回归分析更是人们优先考虑采用的方式。基于此，本文就一元线性回归的MATLAB实现作了一番探讨，给出了多种实现方式，

2、并通过一个实例加以具体展示，在数据处理时可根据自己的需要灵活地加以选用。【关键词】MATLAB程序、一元线性回归、多元线性回归问题、线性回归拟合问题、Subplot、三次样条插值函数。一、提出问题在回归分析中，如果有两个或两个以上的自变量，就称为多元回归。事实上，一种现象常常是与多个因素相联系的，由多个自变量的最优组合共同来预测或估计因变量，比只用一个自变量进行预测或估计更有效，更符合实际。二、简述原理（一）一元线性回归1命令 polyfit最小二乘多项式拟合 p，S=polyfit（x，y，m）多项式y=a1xm+a2xm-1+amx+am+1其中x=（x1，x2，xm）x1xm为（n*1）

3、的矩阵;y为（n*1）的矩阵；p=（a1，a2，am+1）是多项式y=a1xm+a2xm-1+amx+am+1的系数；S是一个矩阵，用来估计预测误差.2命令 polyval多项式函数的预测值Y=polyval（p，x）求polyfit所得的回归多项式在x处的预测值Y；p是polyfit函数的返回值；x和polyfit函数的x值相同。3命令 polyconf 残差个案次序图Y，DELTA=polyconf（p，x，S，alpha）求polyfit所得的回归多项式在x处的预测值Y及预测值的显著性为1-alpha的置信区间DELTA；alpha缺省时为0.05。p是polyfit函数的返回值；x和p

4、olyfit函数的x值相同；S和polyfit函数的S值相同。4. 命令 polytool（x，y，m）一元多项式回归命令5命令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-（k，n-k-1）时拒绝H0，F越大，说明回归方程越显著；与F对应的概率p 时拒绝H0，回归模型成立。

5、Y为n*1的矩阵；X为（ones(n,1),x1,xm）的矩阵；alpha显著性水平（缺省时为0.05）。（二）多元线性回归1命令 regress2命令 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命令 n

6、linfitbeta,R,J=nlinfit(X,Y,model,beta0)X 为n*m矩阵Y为 n维列向量model为自定义函数beta0为估计的模型系数beta为回归系数R为残差2命令 nlintoolnlintool(X,Y,model,beta0,alpha)X 为n*m矩阵Y为 n维列向量model为自定义函数beta0为估计的模型系数alpha显著性水平（缺省时为0.05）3命令 nlparcibetaci=nlparci(beta,R,J)beta为回归系数R为残差返回值为回归系数beta的置信区间4命令 nlpredciY,DELTA=nlpredci(model,X,bet

7、a,R,J)Y为预测值DELTA为预测值的显著性为1-alpha的置信区间；alpha缺省时为0.05。X 为n*m矩阵model为自定义函数beta为回归系数R为残差三、多元线性回归问题在实际经济问题中，一个变量往往受到多个变量的影响。例如，家庭消费支出，除了受家庭可支配收入的影响外，还受诸如家庭所有的财富、物价水平、金融机构存款利息等多种因素的影响，表现在线性回归模型中的解释变量有多个。这样的模型被称为多元线性回归模型。（multivariable linear regression model ）1、多元线性回归模型的一般形式多元线性回归模型的一般形式为：Y =0 +1X1i+2X2i+

8、kXki +i, i=1,2,n （1）其中 k 为解释变量的数目,j (j1,2,k ) 称为回归系（regressioncoefficient)。式也被称为总体回归函数的随机表达式。它的非随机表达式为：Y =0+1X1i +2X2i + k X ki , i=1,2,n (2) j 也被称为偏回归系数（partial regression coefficient)。2、 多元线性回归计算模型Y=0+1x1+2x2+kxk +N (0,2 ) (3) 多元性回归模型的参数估计，同一元线性回归方程一样，也是在要求误差平方和（e)为最小的前提下，用最小二乘法或最大似然估计法求解参数。 设 ( x

9、11, x12 , x1 p , y1 ), ( xn1, xn 2 , xnp , yn ) 是一个样本，用最大似然估计法估计参数：取 b0 , b1 bp ,当 b0b0 , b1b1,b b 达到最小。例1输入以下程序： 输出结果如图1所示： 在 Matlab 命令窗口输入z=b(1)+b(2)*x %运算plot(x,Y,k+,x,z,r) %输出运算 得到预测比较图如图2、3所示：图2图2例 2：（多元线性回归）水泥凝固时放出的热量 y 与水泥中的四种化学成分 x1, x2, x3, x4 有关, 今测得一组数据如下, 试确定多元线性模型%数据%数据%数据%运行方程输出结果如图4所示

10、：图4结果如图5所示：图5在matlab命令框中输入得到如图6所示：图6四、线性回归拟合对于多元线性回归模型：设变量的n组观测值为记 ，则 的估计值为 (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)得到的维列向量b即为(11.2)式给出的回归系数的估计值 b,bint,r,rint,stats=regress(y,x)给出回归系数的估计

11、值b，的95置信区间（向量）bint，残差r以及每个残差的95置信区间（向量）rint；向量stats给出回归的R2统计量和F以及临界概率p的值 如果的置信区间（bint的第行）不包含0，则在显著水平为时拒绝的假设，认为变量是显著的 b, bint,r,rint,stats=regress(y,x,alpha) 给出了bint和rint的100(1-alpha)%的置信区间三次样条插值函数的MATLAB程序matlab的spline x = 0:10; y = sin(x); %插值点xx = 0:.25:10; %绘图点yy = spline(x,y,xx); plot(x,y,o,xx,y

12、y)结果如图7所示图7五、非线性拟合非线性拟合可以用以下命令(同样适用于线形回归分析)：1.beta = nlinfit(X,y,fun,beta0) X给定的自变量数据,Y给定的因变量数据,fun要拟合的函数模型(句柄函数或 者内联函数形式), beta0函数模型中系数估计初值,beta返回拟合后的系数 2.x = lsqcurvefit(fun,x0,xdata,ydata) fun要拟合的目标函数,x0目标函数中的系数估计初值,xdata自变量数据,ydata函数值数据 X拟合返回的系数(拟合结果) nlinfit格式：beta，r，J=nlinfit（x，y，model, beta0）

13、Beta估计出的回归系数r 残差J Jacobian矩阵x，y 输入数据x、y分别为n*m矩阵和n维列向量，对一元非线性回归，x为n维列向量。model是事先用m-文件定义的非线性函数beta0 回归系数的初值例1已知数据： x1=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 ,

14、c , d , e 。(1)对回归模型建立M文件model.m如下: function yy=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); %运算出图(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;%

15、数据beta,r,j = nlinfit(x,y,myfun,beta0)%运算出图输出结果如图8所示图8例题2：混凝土的抗压强度随养护时间的延长而增加，现将一批混凝土作成12个试块，记录了养护日期（日）及抗压强度y（kg/cm2）的数据： 养护时间：x =2 3 4 5 7 9 12 14 17 21 28 56 抗压强度：y =35+r 42+r 47+r 53+r 59+r 65+r 68+r 73+r 76+r 82+r 86+r 99+r 建立非线性回归模型，对得到的模型和系数进行检验。 注明：此题中的+r代表加上一个-0.5,0.5之间的随机数 模型为：y=a+k1*exp(m*x

16、)+k2*exp(-m*x); Matlab程序：x=2 3 4 5 7 9 12 14 17 21 28 56; %数据r=rand(1,12)-0.5; y1=35 42 47 53 59 65 68 73 76 82 86 99; %数据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.5 0.5 0.5 0.5);%运算方法 a=beta(1),k1=beta(2),k2=beta(3),m=

17、beta(4) %数据%test the model%运算方法xx=min(x):max(x); %运算方法yy=a+k1*exp(m*xx)+k2*exp(-m*xx);%运算方法 plot(x,y,o,xx,yy,r) 结果如图9所示：图9lsqnonlin非线性最小二乘（非线性数据拟合）的标准形式为其中：L为常数在MATLAB5.x中，用函数leastsq解决这类问题，在6.0版中使用函数lsqnonlin。设则目标函数可表达为其中：x为向量，F(x)为函数向量。函数 lsqnonlin格式x = lsqnonlin(fun,x0) %x0为初始解向量；fun为，i=1,2,m，fun返

18、回向量值F，而不是平方和值，平方和隐含在算法中，fun的定义与前面相同。x = lsqnonlin(fun,x0,lb,ub)%lb、ub定义x的下界和上界：。x = lsqnonlin(fun,x0,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 = lsqno

19、nlin()%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矩阵。例 求下面非线性最小二乘问题初始解向量为x0=0.3, 0.4。解：先建立函数文件，并保存为myfun.m，由于lsq

20、nonlin中的fun为向量形式而不是平方和形式，因此，myfun函数应由建立： k=1,2,10function F = myfun(x)k = 1:10;F = 2 + 2*k-exp(k*x(1)-exp(k*x(2);然后调用优化程序：x0 = 0.3 0.4；x,resnorm = lsqnonlin(myfun,x0) %主运算方法结果为：Optimization terminated successfully:Norm of the current step is less than OPTIONS.TolXx = 0.2578 0.2578resnorm =%求目标函数值lsq

21、curvefit非线性曲线拟合是已知输入向量xdata和输出向量ydata，并且知道输入与输出的函数关系为ydata=F(x, xdata)，但不知道系数向量x。今进行曲线拟合，求x使得下式成立：在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

22、() %运算方法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() %运算方法参数说明：x0为初始解向量；xdata

23、，ydata为满足关系ydata=F(x, xdata)的数据；lb、ub为解向量的下界和上界，若没有指定界，则lb= ，ub= ；options为指定的优化参数；fun为拟合函数，其定义方式为：x = lsqcurvefit(myfun,x0,xdata,ydata)，其中myfun已定义为 function F = myfun(x,xdata)F = % 计算x处拟合函数值fun的用法与前面相同；resnorm=sum (fun(x,xdata)-ydata).2)，即在x处残差的平方和；residual=fun(x,xdata)-ydata，即在x处的残差；exitflag为终止迭代的条

24、件；output为输出的优化信息；lambda为解x处的Lagrange乘子；jacobian为解x处拟合函数fun的jacobian矩阵。例 求解如下最小二乘非线性拟合问题已知输入向量xdata和输出向量ydata，且长度都是n，拟合函数为即目标函数为其中：初始解向量为x0=0.3, 0.4, 0.1。解：先建立拟合函数文件，并保存为myfun.mfunction F = myfun(x,xdata)F = x(1)*xdata.2 + x(2)*sin(xdata) + x(3)*xdata.3;然后给出数据xdata和ydataxdata = 3.6 7.7 9.3 4.1 8.6 2.

25、8 1.3 7.9 10.0 5.4;ydata = 16.5 150.6 263.1 24.7 208.5 9.9 2.7 163.9 325.0 54.3;x0 = 10, 10, 10; %初始估计值x,resnorm = lsqcurvefit(myfun,x0,xdata,ydata)结果为：Optimization terminated successfully:Relative function value changing by less than OPTIONS.TolFunx = 0.2269 0.3385 0.3021resnorm = 6.2950对于上述这些可化为线性

26、模型的回归问题，一般先将其化为线性模型，然后再用最小二乘法求出参数的估计值，最后再经过适当的变换，得到所求回归曲线。 在熟练掌握最小二乘法的情况下，解决上述问题的关键是确定曲线类型和怎样将其转化为线性模型。 如图10所示图10确定曲线类型一般从两个方面考虑:一是根据专业知识，从理论上推导或凭经验推测、二是在专业知识无能为力的情况下，通过绘制和观测散点图确定曲线大体类型。六、Subplot问题程序：x1=1 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 6 6 6 6 7 8 8 8 8 10 10 10 10 11 11 12 12 13 13 14 15 16

27、16 16 17 20;%数据x2=1 0 1 0 0 1 0 0 0 0 1 1 1 0 1 0 0 0 0 1 0 1 0 1 1 0 1 1 0 0 0 1 1 1 0 0 1 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=1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 0 1 0 1 0 1 1 0 0 0 0 1 0 0 1 0 0 0 0 0

28、1 0 1;%数据x4=0 0 0 1 0 1 1 0 0 1 0 1 0 0 0 0 1 1 0 0 0 0 1 1 0 1 0 0 0 0 1 0 1 1 0 1 0 0 1 1 0 1 1 0 1 0;%数据y=13876 11608 18701 11283 11767 20872 11772 10535 12195 12313 14975 21371 19800 11417 20263 13231 12884 13245 13677 15965 12366 21352 13839 22884 16978 14803 17404 22184 13548 14467 15942 23174

29、 23780 25410 14861 16882 24170 15990 26330 17949 25685 27837 18838 17483 19207 19346;%数据size(x5)x = ones(46,1),x1,x2,x3,x4b,bint,r,rint,stats=regress(y,x,0.05)%运算方法plot(x1,r,r+)hold onfor i=1:1:46 if(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; %交替运算 else k(i)=6; endendk;plot(k,r,g+)输出如图11所示图11第二个x1=1 1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 6 6 6 6 7 8 8 8 8 10 10 10 10 11 11 12 12 13 13 14 15 16 16 16 17 20;%数据x2=1 0 1 0 0 1 0 0 0 0 1 1 1 0 1