Matlab解决线性回归问题论文.docx
《Matlab解决线性回归问题论文.docx》由会员分享，可在线阅读，更多相关《Matlab解决线性回归问题论文.docx（27页珍藏版）》请在冰豆网上搜索。
Matlab解决线性回归问题论文
Matlab解决线性回归问题
2010级电子信息工程班
【摘要】MATLAB和Mathematica、Maple并称为三大数学软件。
它在数学类科技应用软件中在数值计算方面首屈一指。
MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。
回归分析，是对现有数据进行处理、从中发现有用信息的一种重要手段。
而线性回归，特别是一元线性回归分析更是人们优先考虑采用的方式。
基于此，本文就一元线性回归的MATLAB实现作了一番探讨，给出了多种实现方式，并通过一个实例加以具体展示，在数据处理时可根据自己的需要灵活地加以选用。
【关键词】MATLAB程序、一元线性回归、多元线性回归问题、线性回归拟合问题、Subplot、三次样条插值函数。
一、提出问题
在回归分析中，如果有两个或两个以上的自变量，就称为多元回归。
事实上，一种现象常常是与多个因素相联系的，由多个自变量的最优组合共同来预测或估计因变量，比只用一个自变量进行预测或估计更有效，更符合实际。
二、简述原理
（一）一元线性回归
1．命令polyfit最小二乘多项式拟合
[p，S]=polyfit（x，y，m）
多项式y=a1xm+a2xm-1+…+amx+am+1
其中x=（x1，x2，…，xm）x1…xm为（n*1）的矩阵;
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和polyfit函数的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，回归模型成立。
Y为n*1的矩阵；
X为（ones（n,1）,x1,…,xm）的矩阵；
alpha显著性水平（缺省时为0.05）。
（二）多元线性回归
1．命令regress
2．命令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为自定义函数
beta0为估计的模型系数
beta为回归系数
R为残差
2．命令nlintool
nlintool（X,Y,’model’,beta0,alpha）
X为n*m矩阵
Y为n维列向量
model为自定义函数
beta0为估计的模型系数
alpha显著性水平（缺省时为0.05）
3．命令nlparci
betaci=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、多元线性回归模型的一般形式
多元线性回归模型的一般形式为：
Y=0+1X1i+2X2i++kXki+i,i=1,2,,n
（1）
其中k为解释变量的数目,j（j1,2,k）称为回归系（regressioncoefficient）。
式也被称为总体回归函数的随机表达式。
它的非随机表达式为：
Y=0+1X1i+2X2i++kXki,i=1,2,n
（2）
®j也被称为偏回归系数（partialregressioncoefficient）。
2、多元线性回归计算模型
Y=0+1x1+2x2+kxk+∼N（0,2）（3）
多元性回归模型的参数估计，同一元线性回归方程一样，也是在要求误差平方和（Σe）为最小的前提下，用最小二乘法或最大似然估计法求解参数。
设（x11,x12,x1p,y1）,,（xn1,xn2,xnp,yn）是一个样本，用最大似然估计法估计参数：
取b0,b1bp,当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所示：
图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）给出回归系数
的估计值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,yy）
结果如图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）
Beta估计出的回归系数
r残差
JJacobian矩阵
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,c,d,e。
（1）对回归模型建立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）;%运算出图
（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,beta0）%运算出图
输出结果如图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（1,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所示：
图9
lsqnonlin
非线性最小二乘（非线性数据拟合）的标准形式为
其中：
L为常数
在MATLAB5.x中，用函数leastsq解决这类问题，在6.0版中使用函数lsqnonlin。
设
则目标函数可表达为
其中：
x为向量，F（x）为函数向量。
函数lsqnonlin
格式x=lsqnonlin（fun,x0）%x0为初始解向量；fun为
，i=1,2,…,m，fun返回向量值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]=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矩阵。
例求下面非线性最小二乘问题
初始解向量为x0=[0.3,0.4]。
解：
先建立函数文件，并保存为myfun.m，由于lsqnonlin中的fun为向量形式而不是平方和形式，因此，myfun函数应由
建立：
k=1,2,…,10
functionF=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.TolX
x=
0.25780.2578
resnorm=%求目标函数值
lsqcurvefit
非线性曲线拟合是已知输入向量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（…）%运算方法
[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，ydata为满足关系ydata=F（x,xdata）的数据；
lb、ub为解向量的下界和上界
，若没有指定界，则lb=[]，ub=[]；
options为指定的优化参数；fun为拟合函数，其定义方式为：
x=lsqcurvefit（@myfun,x0,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，拟合函数为
即目标函数为
其中：
初始解向量为x0=[0.3,0.4,0.1]。
解：
先建立拟合函数文件，并保存为myfun.m
functionF=myfun（x,xdata）
F=x
（1）*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.TolFun
x=0.22690.33850.3021
resnorm=6.2950
对于上述这些可化为线性模型的回归问题，一般先将其化为线性模型，然后再用最小二乘法求出参数的估计值，最后再经过适当的变换，得到所求回归曲线。
在熟练掌握最小二乘法的情况下，解决上述问题的关键是确定曲线类型和怎样将其转化为线性模型。
如图10所示
图10
确定曲线类型一般从两个方面考虑:
一是根据专业知识，从理论上推导或凭经验推测、二是在专业知识无能为力的情况下，通过绘制和观测散点图确定曲线大体类型。
六、Subplot问题
程序：
x1=[111112222333344445556666788881010101011111212131314151616161720]';%数据
x2=[1010010000111010000101011011000111001010110000]';%数据
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=[0001011001010000110000110100001011010011011010]';%数据
y=[13876116081870111283117672087211772105351219512313149752137119800114172026313231128841324513677159651236621352138392288416978148031740422184135481446715942231742378025410148611688224170159902633017949256852783718838174831920719346];%数据
size（x5）
x=[ones（46,1）,x1,x2,x3,x4]
[b,bint,r,rint,stats]=regress（y,x,0.05）%运算方法
plot（x1,r,'r+'）
holdon
fori=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;
end
end
k;
plot（k,r,'g+'）
输出如图11所示
图11
第二个
x1=[111112222333344445556666788881010101011111212131314151616161720]';%数据
x2=[101001000011101
升级会员 