多元回归程序MATLAB程序

matlab回归（拟合）总结前言1、学三条命令polyfit(x,y,n)---拟合成一元幂函数（一元多次）regress(y,x)----可以多元，nlinfit(x,y,’fun’,beta0)(可用于任何类型的函数，任意多元函数，应用范围最广，最万能的)2、同一个问题，这三条命令都可以使用，但结果肯定是不同的，因为拟合的近似结果，没有唯一的标准的答案。相当于咨询多个专家。3、回归的操作步骤：根据图形（实际点），选配一条恰当的函数形式（类型）---需要数学理论与基础和经验。（并写出该函数表达式的一般形式，含待定系数）------选用某条回归命令求出所有的待定系数。所以可以说，回归就是求待定系数的过程（需确定函数的形式）一、多元回归分析对于多元线性回归模型(其实可以是非线性，它通用性极高)：exxypp110设变量12,,,pxxxy的n组观测值为12(,,,)1,2,,iiipixxxyin记npnnppxxxxxxxxxx212222111211111，nyyyy21，则p10的估计值为排列方式与线性代数中的线性方程组相同（），拟合成多元函数---regress使用格式：左边用b=[b,bint,r,rint,stats]右边用=regress(y,x)或regress(y,x,alpha)---命令中是先y后x,---须构造好矩阵x(x中的每列与目标函数的一项对应)---并且x要在最前面额外添加全1列/对应于常数项---y必须是列向量---结果是从常数项开始---与polyfit的不同。）其中：b为回归系数，的估计值(第一个为常数项)，bint为回归系数的区间估计，r:残差，rint:残差的置信区间，stats:用于检验回归模型的统计量，有四个数值：相关系数r2、F值、与F对应的概率p和残差的方差（前两个越大越好，后两个越小越好），alpha:显著性水平（缺省时为0.05，即置信水平为95%），（alpha不影响b,只影响bint(区间估计)。它越小，即置信度越高，则bint范围越大。显著水平越高，则区间就越小）（返回五个结果）---如有n个自变量-有误（n个待定系数），则b中就有n+1个系数（含常数项，---第一项为常数项）(b---b的范围/置信区间---残差r---r的置信区间rint-----点估计----区间估计如果i的置信区间（bint的第1i行）不包含0，则在显著水平为时拒绝0i的假设，认为变量ix是显著的．*******(而rint残差的区间应包含0则更好)。b,y等均为列向量,x为矩阵(表示了一组实际的数据)必须在x第一列添加一个全1列。----对应于常数项。相关系数r2越接近1，说明回归方程越显著；(r2越大越接近1越好)F越大，说明回归方程越显著；（F越大越好）与F对应的概率p越小越好，一定要Pa时拒绝H0而接受H1，即回归模型成立。乘余（残差）标准差（RMSE）越小越好(此处是残差的方差，还没有开方)（前两个越大越好，后两个越小越好）重点：regress(y,x)重点与难点是如何加工处理矩阵x。y是函数值，一定是只有一列。也即目标函数的形式是由矩阵X来确定如s=a+b*x1+c*x2+d*x3+e*x1^2+f*x2*x3+g*x1^2,一定有一个常数项，且必须放在最前面（即x的第一列为全1列）X中的每一列对应于目标函数中的一项（目标函数有多少项则x中就有多少列）X=[ones,x1,x2,x3,x1.^2,x2.*x3，x1.ˆ2](剔除待定系数的形式)regress:y/x顺序，矩阵X需要加工处理nlinfit:x/y顺序，X/Y就是原始的数据，不要做任何的加工。(即regress靠矩阵X来确定目标函数的类型形式（所以X很复杂，要作很多处理）而nlinfit是靠程序来确定目标函数的类型形式（所以X就是原始数据，不要做任何处理）例1测16名成年女子的身高与腿长所得数据如下：身高143145146147149150153154155156157158159160162164腿长88858891929393959698979698991001021401451501551601658486889092949698100102配成y=a+b*x形式x=[143145146147149150153154155156157158159160162164]';y=[8885889192939395969897969899100102]';plot(x,y,'r+')z=x;x=[ones(16,1),x];----常数项[b,bint,r,rint,stats]=regress(y,x);---处结果与polyfit(x,y,1)相同b,bint,stats得结果：b=bint=-16.0730-33.70711.5612------每一行为一个区间0.71940.60470.8340stats=0.9282180.95310.0000即7194.0ˆ,073.16ˆ10；0ˆ的置信区间为[-33.7017，1.5612],1ˆ的置信区间为[0.6047,0.834];r2=0.9282,F=180.9531,p=0.0。p0.05,可知回归模型y=-16.073+0.7194x成立.[b,bint,r,rint,stats]=regress(Y,X,0.05);-----结果相同[b,bint,r,rint,stats]=regress(Y,X,0.03);polyfit(x,y,1)-----当为一元时（也只有一组数），则结果与regress是相同的，只是命令中x,y要交换顺序，结果的系数排列顺序完全相反，x中不需要全1列。ans=0.7194-16.0730--此题也可用polyfit求解，杀鸡用牛刀，脖子被切断。3、残差分析，作残差图：rcoplot(r,rint)246810121416-5-4-3-2-101234ResidualCaseOrderPlotResidualsCaseNumber从残差图可以看出，除第二个数据外，其余数据的残差离零点均较近，且残差的置信区间均包含零点，这说明回归模型y=-16.073+0.7194x能较好的符合原始数据，而第二个数据可视为异常点(而剔除)4、预测及作图：plot(x,y,'r+')holdona=140:165;b=b(1)+b(2)*a;plot(a,b,'g')1401451501551601658486889092949698100102例2观测物体降落的距离s与时间t的关系，得到数据如下表，求s关于t的回归方程2ˆctbtast(s)1/302/303/304/305/306/307/30s(cm)11.8615.6720.6026.6933.7141.9351.13t(s)8/309/3010/3011/3012/3013/3014/30s(cm)61.4972.9085.4499.08113.77129.54146.48法一：直接作二次多项式回归t=1/30:1/30:14/30;s=[11.8615.6720.6026.6933.7141.9351.1361.4972.9085.4499.08113.77129.54146.48];[p,S]=polyfit(t,s,2)p=489.294665.88969.1329得回归模型为：1329.98896.652946.489ˆ2tts方法二----化为多元线性回归：2ˆctbtast=1/30:1/30:14/30;s=[11.8615.6720.6026.6933.7141.9351.1361.4972.9085.4499.08113.77129.54146.48];T=[ones(14,1),t',(t.^2)']%？？？是否可行？？？等验证...----因为有三个待定系数，所以有三列，始于常数项[b,bint,r,rint,stats]=regress(s',T);b,statsb=9.132965.8896489.2946stats=1.0e+007*0.00001.037800.0000得回归模型为：22946.4898896.651329.9ˆtts%结果与方法1相同polyfit------一元多次regress----多元一次---其实通过技巧也可以多元多次regress最通用的，万能的，表面上是多元一次，其实可以变为多元多次且任意函数，如x有n列（不含全1列），则表达式中就有n+1列（第一个为常数项，其他每项与x的列序相对应）。例3设某商品的需求量与消费者的平均收入、商品价格的统计数据如下，建立回归模型，预测平均收入为1000、价格为6时的商品需求量.需求量10075807050659010011060收入10006001200500300400130011001300300价格5766875439选择纯二次模型，即2222211122110xxxxy----用户可以任意设计函数x1=[10006001200500300400130011001300300];x2=[5766875439];y=[10075807050659010011060]';X=[ones(10,1)x1'x2'(x1.^2)'(x2.^2)'];%10指有10组数据，x1'x2'(x1.^2)'(x2.^2)'时方程的自变量[b,bint,r,rint,stats]=regress(y,X)b,statsb=110.53130.1464-26.5709-0.00011.8475stats=0.970240.66560.000520.5771故回归模型为：2221218475.10001.05709.261464.05313.110xxxxy剩余标准差为4.5362,说明此回归模型的显著性较好.三、非线性回归（拟合）使用格式：beta=nlinfit(x,y,‘程序名’,beta0)[beta,r,J]=nlinfit(X,y,fun,beta0)X给定的自变量数据,Y给定的因变量数据,fun要拟合的函数模型(句柄函数或者内联函数形式),beta0函数模型中待定系数估计初值（即程序的初始实参）beta返回拟合后的待定系数其中beta为估计出的回归系数；r为残差；J为Jacobian矩阵输入数据x、y分别为n*m矩阵和n维列向量，对一元非线性回归，x为n维列向量。可以拟合成任意函数。最通用的，万能的命令x,y顺序，x不需要任何加工，直接用原始数据。---所编的程序一定是两个形参（待定系数/向量，自变量/矩阵：每一列为一个自变量）结果要看残差的大小和是否有警告信息，如有警告则换一个b0初始向量再重新计算。本程序中也可能要用.*./.^如结果中有警告信息，则必须多次换初值来试算.难点是编程序与初值存在的问题：不同的beta0,则会产生不同的结果，如何给待定系数的初值以及如何分析结果的好坏，如出现警告信息，则换一个待定系数试一试。因为拟合本来就是近似的，可能有多个结果。1:重点（难点）是预先编程序（即确定目标函数的形式，而regress的目标函数由x矩阵来确定，其重难点为构造矩阵a）2:x/y顺序—列向量----x/y是原始数据，不要做任何修改3：编程：一定两个形参（beta,x）a=beta(1);b=beta(2);c=beta(3);…x1=x(:,1);x2=x(:,2);x3=x(:,3);即每一列为一个自变量4：regress/nlinfit都是列向量5：regress:有n项（n个待定系数），x就有n列；nlinfit:有m个变量则x就有m列例1已知数据：x1=[0.5,0.4,0.3,0.2,0.1];x2=[0.3,0.5,0.2,0.4,0.6];x