求多元线形回归方程及预报

907求多元线形回归方程及预报一．功能x1，x2,x3,………..xp为自变量，Y为随机变量，求线形回归方程Y=b0+xb11+…+xbpp+其中,,10bb。。。，bp为常数，是随机变量，且—N（0，26）来描述Y与X的变化规律。并用T检验法检验线形回归是否显著。如果线性回归显著，可用经验回归平面方程对Y作出预报，并给出预报值的置信区间。二．算法间介[16]，[15](1)求回归方程设x1,x2…xp是确定变量，Y是随机变量，他们之间有关系Y=b0+xb11+…+xbpp+其中,,10bb。。。，bp为常数，是随机变量，且—N（0，26），这是P元线性回归模型，我们讨论P1的情形。作n次独立试验，得到n组数据（xk1+2xk+,…),Yxpkp(k=1,2,…,n)记Xj=nkXn11kjj=1,2,…,p,则（1）式可写为Y=b1（x1-X1）+…+(bpxpXp)+其中同于（1）中之，而=Xbb101+…+Xbpp对上面得到的几组试验数据，便有Y=Xbk11(-X1)+…+bp(Xkp-Xp)+k其中独立分布：k—N（0，26）。为了用最小乘法求（2）中，bbbp,...,,21的估计值，我们引如下述符号Y=n1niY11Y=XXXXXXXXXXXXXXXXXXpnpnnPPPP......,1......................,1....,1.22112222121212111,,....,2,1,,)(),()(`11pkjyXXLXXXXLinijijjykiknijijjkA=XXtippipppiiXXXXXXXXXXXXn211111111112111)(...)()(0...................................))(....)(000...0(=LLLLLLLLLpppPppn....0.......................0...00...00212222111211=Ln00A为准对角阵，子块L是P介实对称可逆阵。，即A为L=PPPPPPLLLLLLLLL212222111211B=YXT)()(111PPIIIIXXYXXYY=LPYYLYI1由此得正规方程组lon0=^^2^1^pbbbu=pyyyillly21利用分块乘法得niiyu^及l^^2^1pbbb=pyyylll21从而得估计计算公式:u=^y，^^^21...bbbp=1LLLLpyyy...21经验回归方程为^^1^1^,....xbxbppY（2）假设检验线性回归的显著性检验.在线性模型(2)中作假设He2b1=0,b2=0,…bp=0由QI=2)]ˆ(2)ˆ[(1YIYIYYIninininiYYIIYYIYIYIYYI111))(ˆ(22)ˆ(2)ˆ(利正规方程组,可知右边第三项为0,从而Qi=niniYIYIYYI112)ˆ(2)ˆ(记Qe=niiYYi12)ˆ(Q2=nipjBbYiY10ˆ2)ˆ(j其中Bj(j=1,2,…,p)是正规方程组的右端项.QE称为剩余离差(平方和),是由试验或2引起的误差;称为回归离差Q2(平方和),是由线性回归引起.由分解定理知:21QE服从自由度为n-p-1的X2分布21Q1服从自由度为p的X2分布记SE2=1pnQE,S22=pQ2则F=SSE222~F(p,n-p-1)由该式计算出F值,对给出的显著性水平,F(p,n-p-1),F=F(p,n-p-1)则拒H0,即认为线性回归显著;则接受H0,认为线性回归不显著.用复相关系数R,且R2=12QQ,也可说明Y与自变量的密切程度,02R1,2R越大越密切回归系数的显著性检验若经检验线性回归显著，即说明回归系数b1,b2,…,bp不全为0,但是不能说明每个自变量对Y都是重要的，如果某个系数为0或接近于0相应的自变量对Y不起作用或作用很小，可以忽略，因而检验每个回归系数bj(1=j=p)是否为0，相当于检验相应的Xj对Y的值是否起作用。在线性回归模型上作假设H0：bj=0，其中j固定，1=j=p,这里用b，作检验。在H0成立的前提下，可知=,j=1,2,3,…p服从自由度为n-p-1的t分布。C是C=A-1的主对角线上第j个元素。给定显著性水平a，一次抽样后，若=则拒绝假设H0，即认为b，显著地不等于0或说b于0有显著差别；否则接受H0，即认为b于0无显著差别。（3）预报所谓预报，就是给出X1=X01，X2=X02,…,Xp=Xop,对Y的值Y0作去件估计，即求出Y0，并指出它的置信区间（置信度1-a给出时）用=++…+对Y0作点估计。为了作区间估计，须讨论Y0-Y0的概率分布。可知E（-）=0D（-）=++-）（其中Cij是C=A-1第I行第j列元素，62可用下式计算记))((1101011120jjpjiipiXXXXCnd则()^00YY~N(0,2d20)从而U=0^00dYY~N(0,1)由此知T=)1(~0^00pndYY给出置信概率1-a，那么Y0的置信区间为（+ta/2（n-p-1）0d式中da可用（18）式计算后经开方得到。三、程序说明程序中分主程序与子程序，子程序有五个。子程序（1）求线性回归方程子程序（2）线性回归的显著性检验子程序（3）回归系数的显著检验子程序（4）预报子程序（5）计算正规方程组系数矩阵之逆阵注意，当检验出线性回归不显著时，就不再执行子程序（3）和（4）四、程序使用说明1、输入参数p实数，自变量个数N实数，实验次数R实数，给定的显著性水平，可取R=0.05FR实数FR=Fa(p,n-p-1),查表得到TR实数，TR=ta/2(n-p-1),由表给出A(N,P)二维实数组，P1=p+1，存放原始数据2、输入参数（1）对1中输入的参数全部打印输出（2）计算回归方程的一些中间运算结果及最后结果。X0(P1)一维实数组，存放各变量的均值L（p,p）正规方程组的系数矩阵LY（p）正规方程的右端项C（p,p）L（p,p）的逆矩阵B（p）回归系数回归方程（3）检验线性回归显著性的一些中间运算结果和检验结果方差分析表检验结果RR复相关系数，0RR1（4）检验回归系数显著的结果B1与0有显著差别（5）预报结果3、中间工作单元略10********************20907多元线莘回归及预报*30*******************40INPUT”P,N=”P,N:P1=P+150DIMA(N,P1),X)(P1),L(P,P),C(P,),LY(P),S(P,P1)60DIMB(P)70FORI=1TON80FORJ=1TOP190READA(I,J)100NEXTJ,I110LPRINTTAB(3);”例：“；“P=”；P；“N=”；N：LPRINT120LPRINT“序号”；TAB（9）；“X1”；TAB（17）；“X2”；TAB（25）；“X3”；130LPRINTTAB（35）；“Y”140FORI=1TOP1150LPRINTI；160FORJ=1TOP1170LPRINTTAB（8*J）；A（I，J）；180NEXTJ190LRINT200NEXTI210LRINT：LPRINT220GOSUB1：100230LPRINT240LPRINTTAB（3）：“回归系数和回归方程”250FORJ=1TOP260LPRINT“B（”；J；“）=”；USING“#。#####”；B（J）；270IFINT（I/4）*4=ITHENLPRINT280NEXTJ290LPRINT300LPRINT“Y=”；USING“#。#####”；B0；310FORI-1TOP320IFB（1）》0THENLPRINT“+”；330LPRINTUSING“#。####”；B（1）；340LPRINT“*X”；I；350IFINT（I/3）*3=ITHENLPRINT360NEXTI370LPRINT380****390INPUT“R，Fr，Tr=”；R，FR，TR400GOSUB1770410LPRINT420LPRINTTAB(3);”R=”;R;”Fr=”FR;”Tr=”TR:LPRINT430LPRINTTAB(15);”方差分析表：“440LPRINT“来源“;TAB(8);”离差”;TAB(16);”自由度”；450LPRINTTAB(26);”均方离差”;TAB(38);”F值”460LPRINT“Q1”;TAB(5);Q1;TAB(18);P;TAB(25);S1;TAB(35);F470LPRINT“Q2”;TAB(5);Q2;TAB(18);N-P1;TAB(25);S2480LPRINT“Qt”;TAB(5);QT;TAB(18);N-1490LPRINT:LPRINTTAB(3);”RR=”;RR:LPRINT500IFF=FRTHEN510ELSE520510LPRINT“线形回归显著”:GOTO530520LPRINT“线形回归不显著”530IFFFRTHEN1000540GOSUB1930550LPRINT:LPRINTTAB(3);”T(J):”560FORJ=1TOP570LPRINTT(J);:LPRINT”“;580NEXTJ590LPRINT:LPRINT600FORJ=1TOP610IFABS(T(J))=TRTHEN620ELSE630620LPRINT“B(“;J;”)与0有显著差异“:GOTO640630LPRINT“B(“;J;”)与0无显著差异”640NEXTJ650LPRINT660****670INPUT“WW=”:WW680IFWW=0THEN690ELSE700690LPRINTTAB(3);”WW=”WW:GOTO1000700FORJ=1TOP710INPUT“X(J)=”;X(J)720NEXTJ730GOSUB1980740LPRINT:LPRINTTAB(3);”预报：X(J):”750FORJ=1TOP760LPRINTX(J);:LPRINT;770IF(J/5)*5=JTHENLPRINT780NEXTJ790LPRINT800LPRINTTAB(12);”Y0=”;Y0810LPRINTTAB(3);Y1;”Y0”;Y2820LPRINT:GOTO670830DATA2,18,50,4.3302,7,9,40,3.6485,5,14,46,4.4830840DATA12,3,43,5.5468,1,20,64,5.4970,3,12,40,3.1125850DATA3,17,64,5.1182,6,5,39,3.8759,7,8,37,4.6700860DATA0,23,55,4.9536,3,16,60,5.0060,0,18,40,5.2701870DATA8,4,50,5.3772,6,14,51,5.4849,0,21,51,4.5960880DATA3,14,51,5.6645,7,12,56,6.0795,16,0,48,3.2194890DATA6,16,45,5.8076,0,15,52,4.7308,9,0,40,4.6805900DATA4,