多元回归程序MATLAB程序.doc

1、matlab回归（拟合）总结前言1、学三条命令polyfit(x,y,n)-拟合成一元幂函数（一元多次）regress(y,x)-可以多元， nlinfit(x,y,fun,beta0) (可用于任何类型的函数，任意多元函数，应用范围最广，最万能的)2、同一个问题，这三条命令都可以使用，但结果肯定是不同的，因为拟合的近似结果，没有唯一的标准的答案。相当于咨询多个专家。3、回归的操作步骤：根据图形（实际点），选配一条恰当的函数形式（类型）-需要数学理论与基础和经验。（并写出该函数表达式的一般形式，含待定系数）-选用某条回归命令求出所有的待定系数。所以可以说，回归就是求待定系数的过程（需确定函数的

2、形式）一、多元回归分析对于多元线性回归模型(其实可以是非线性，它通用性极高)：设变量的n组观测值为记 ，则 的估计值为排列方式与线性代数中的线性方程组相同（），拟合成多元函数-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:

3、残差 ，rint: 残差的置信区间，stats: 用于检验回归模型的统计量，有四个数值：相关系数r2、F值、与F对应的概率p和残差的方差（前两个越大越好，后两个越小越好），alpha: 显著性水平（缺省时为0.05，即置信水平为95%），（alpha不影响b,只影响bint(区间估计)。它越小，即置信度越高，则bint范围越大。显著水平越高，则区间就越小）（返回五个结果）-如有n个自变量-有误（n个待定系数），则b 中就有n+1个系数（含常数项，-第一项为常数项）(b-b的范围/置信区间-残差r-r的置信区间rint-点估计-区间估计 如果的置信区间（bint的第行）不包含0，则在显著水平为时

4、拒绝的假设，认为变量是显著的*(而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

5、+c*x2+d*x3+e*x12+f*x2*x3+g*x12,一定有一个常数项，且必须放在最前面（即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名成年女子的身

6、高与腿长所得数据如下：身高143145146147149150153154155156157158159160162164腿长8885889192939395969897969899100102配成y=a+b*x形式 x=143 145 146 147 149 150 153 154 155 156 157 158 159 160 162 164; y=88 85 88 91 92 93 93 95 96 98 97 96 98 99 100 102; plot(x,y,r+)z=x;x=ones(16,1),x;-常数项b,bint,r,rint,stats=regress(y,x);-处结

7、果与polyfit(x,y,1)相同b,bint,stats得结果：b = bint = -16.0730 -33.7071 1.5612-每一行为一个区间 0.7194 0.6047 0.8340stats = 0.9282 180.9531 0.0000即；的置信区间为-33.7017，1.5612, 的置信区间为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

8、=regress(Y,X,0.03);polyfit(x,y,1)-当为一元时（也只有一组数），则结果与regress是相同的，只是命令中x,y要交换顺序，结果的系数排列顺序完全相反，x中不需要全1列。ans =0.7194 -16.0730-此题也可用polyfit求解，杀鸡用牛刀，脖子被切断。3、残差分析，作残差图：rcoplot(r,rint)从残差图可以看出，除第二个数据外，其余数据的残差离零点均较近，且残差的置信区间均包含零点，这说明回归模型 y=-16.073+0.7194x能较好的符合原始数据，而第二个数据可视为异常点(而剔除)4、预测及作图：plot(x,y,r+) hold

9、on a=140:165; b=b(1)+b(2)*a;plot(a,b,g)例2 观测物体降落的距离s与时间t的关系，得到数据如下表，求s关于t的回归方程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.86 15.67 20.60 26.69 33.71 4

10、1.93 51.13 61.49 72.90 85.44 99.08 113.77 129.54 146.48;p,S=polyfit(t,s,2)p =489.2946 65.8896 9.1329得回归模型为 ：方法二-化为多元线性回归：t=1/30:1/30:14/30;s=11.86 15.67 20.60 26.69 33.71 41.93 51.13 61.49 72.90 85.44 99.08 113.77 129.54 146.48;T=ones(14,1), t, (t.2) %？是否可行？等验证.-因为有三个待定系数，所以有三列，始于常数项b,bint,r,rint,st

11、ats=regress(s,T);b,statsb = 9.1329 65.8896 489.2946stats =1.0e+007 * 0.0000 1.0378 0 0.0000得回归模型为 ： %结果与方法1相同polyfit-一元多次regress-多元一次-其实通过技巧也可以多元多次regress最通用的，万能的，表面上是多元一次，其实可以变为多元多次且任意函数，如x有n列（不含全1列），则表达式中就有n+1列（第一个为常数项，其他每项与x的列序相对应）。例3 设某商品的需求量与消费者的平均收入、商品价格的统计数据如下，建立回归模型，预测平均收入为1000、价格为6时的商品需求量.需

12、求量10075807050659010011060收入1000600 1200500300400130011001300300价格5766875439选择纯二次模型，即 -用户可以任意设计函数x1=1000 600 1200 500 300 400 1300 1100 1300 300;x2=5 7 6 6 8 7 5 4 3 9;y=100 75 80 70 50 65 90 100 110 60;X=ones(10,1) x1 x2 (x1.2) (x2.2); %10指有10组数据，x1 x2 (x1.2) (x2.2)时方程的自变量b,bint,r,rint,stats=regress

13、(y,X)b,stats b = 110.5313 0.1464 -26.5709 -0.0001 1.8475stats = 0.9702 40.6656 0.0005 20.5771故回归模型为：剩余标准差为4.5362, 说明此回归模型的显著性较好.三、非线性回归（拟合）使用格式：beta = nlinfit(x,y, 程序名,beta0) beta,r,J = nlinfit(X,y,fun,beta0)X给定的自变量数据,Y给定的因变量数据,fun要拟合的函数模型(句柄函数或者内联函数形式), beta0函数模型中待定系数估计初值（即程序的初始实参）beta返回拟合后的待定系数其中b

14、eta为估计出的回归系数；r为残差；J为Jacobian矩阵输入数据x、y分别为n*m矩阵和n维列向量，对一元非线性回归，x为n维列向量。可以拟合成任意函数。最通用的，万能的命令x,y顺序，x不需要任何加工，直接用原始数据。-所编的程序一定是两个形参（待定系数/向量，自变量/矩阵：每一列为一个自变量）结果要看残差的大小和是否有警告信息，如有警告则换一个b0初始向量再重新计算。本程序中也可能要用.* ./ .如结果中有警告信息，则必须多次换初值来试算.难点是编程序与初值存在的问题：不同的beta0,则会产生不同的结果，如何给待定系数的初值以及如何分析结果的好坏，如出现警告信息，则换一个待定系数试

15、一试。因为拟合本来就是近似的，可能有多个结果。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

16、.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如下:function yy=myfun(beta,x) %一定是两个参数：系数和自变量-一个向量/一个矩阵a=beta(1)b=beta(2)c=

17、beta(3)d=beta(4)e=beta(5)x1=x(:,1); %系数是数组，b(1),b(2),b(n)依次代表系数1, 系数2, 系数nx2=x(:,2); %自变量x是一个矩阵,它的每一列分别代表一个变量,有n列就可以最多nx3=x(:,3);yy=beta(1)+beta(2)*x2+beta(3)*x3+beta(4)*(x2.2)+beta(5)*(x3.2);(b(i)与待定系数的顺序关系可以任意排列，并不是一定常数项在最前，只是结果与自己指定的相对应)（x一定是一列对应一个变量，不能x1=x(1),x2=x(2),x3=x(3)）(2)主程序如下:x=0.5,0.4,0

18、.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; %每一列为一个变量,如果是倒入数据矩阵,只能把x的数据倒进去，不能把全部数据都倒进去，然后选某几列y=0.785,0.703,0.583,0.571,0.126;beta0=1,1, 1,1, 1,1; %有多少个待定系数，就给多少个初始值。beta,r,j = nlinfit(x,y,myfun,beta0)beta = -0.4420 5.5111 0.3837 -8.1734 -0.1340例2混凝土的抗压强度随养护时间的延长而增加，现将一批混凝土作成12个试块，记录了养护日期（日）及

19、抗压强度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)+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 6

20、5 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); %初值为0.2也可以，如为1则不行，则试着换系数初值-此处为一元，x,y行/列向量都可以a=beta(1),k1=beta(2),k2=beta(3),m=beta(4)%test the modelxx=min(x):max(x);2：56yy=a+k1*exp(m*xx)+k2*exp(-m*xx

21、); plot(x,y,o,xx,yy,r) %,xx,yy,r是画曲线，相当于拟合结果： a = 87.5244 k1 = 0.0269 k2 = -63.4591 m = 0.1083 图形：例3 出钢时所用的盛钢水的钢包，由于钢水对耐火材料的侵蚀，容积不断增大.我们希望知道使用次数与增大的容积之间的关系.对一钢包作试验，测得的数据列于下表：使用次数增大容积使用次数增大容积234567896.428.209.589.509.7010.009.939.991011121314151610.4910.5910.6010.8010.6010.9010.76对将要拟合的非线性模型y= aeb/x，

22、（如再加y= c*sin(x)+aeb/x）建立m-文件volum.m如下：function yhat=volum(beta,x) yhat=beta(1)*exp(beta(2)./x);或function f= volum (beta,x)a=beta(1);b=beta(2);f=a*exp(b./x);2、输入数据：主程序：x=2:16;y=6.42 8.20 9.58 9.5 9.7 10 9.93 9.99 10.49 10.59 10.60 10.80 10.60 10.90 10.76;plot(x,y,*-) beta0=8 2; beta,r ,J=nlinfit(x,y,

23、zhang1,beta0)3、求回归系数：beta = 11.6037 -1.0641即得回归模型为：4、预测及作图： plot(x,y,ro) hold onxx=2:0.05:16; yy=zhang1(beta,xx);%-通过调用用户自编的函数plot(xx,yy,g) % 拟合成线或者YY,delta=nlpredci(zhang1,x,beta,r ,J)plot(x,y,k+,x,YY,r)或plot(x,y,ro)hold onxx=2:0.05:16;yy=beta(1)*exp(beta(2)./xx); plot(xx,yy,g)例4 财政收入预测问题：财政收入与国民收入

24、、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。下表列出了1952-1981年的原始数据，试构造预测模型。财政收入预测问题：财政收入与国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资等因素有关。下表列出了1952-1981年的原始数据，试构造预测模型。年份国民收入（亿元）工业总产值(亿元)农业总产值（亿元）总人口（万人）就业人口（万人）固定资产投资（亿元）财政收入(亿元)19525983494615748220729441841953586455475587962136489216195470752049160266218329724819557375585

25、296146522328982541956825715556628282301815026819578377985756465323711139286195810281235598659942660025635719591114168150967207261733384441960107918704446620725880380506196175711564346585925590138271196267796446167295251106623019637791046514691722664085266196494312505847049927736129323196511521581632

26、725382867017539319661322191168774542298052124661967124916476977636830814156352196811871565680785343191512730319691372210168880671332252074471970163827477678299234432312564197117803156790852293562035563819721833336578987177358543546581973197836848558921136652374691197419933696891908593736939365519752

27、1214254932924213816846269219762052430995593717388344436571977218949259719497439377454723197824755590105896259398565509221979270260651150975424058156489019802791659211949870541896568826198129276862127310007273280496810解 设国民收入、工业总产值、农业总产值、总人口、就业人口、固定资产投资分别为x1、x2、x3、x4、x5、x6，财政收入为y，设变量之间的关系为：y= ax1+bx2+cx3+dx4+ex5+fx6使用非线性回归方法求解。1 对回归模型建立M文件model.m如下: function yy=model(beta0,X) %一定是两个参数，第一个为系数数组，b(1),b(2),b(n) %分别代表每个系数，而第二个参数代表所有的自变量， %是一个矩阵，它的每一列分别代表一个自变量。