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

1、907求多元线形回归方程及预报一功能 x1，x2, x3,.xp为自变量，Y为随机变量，求线形回归方程Y=+ 其中。，为常数，是随机变量，且N（0，） 来描述Y与X的变化规律。并用T检验法检验线形回归是否显著。如果线性回归显著，可用经验回归平面方程对Y作出预报，并给出预报值的置信区间。二算法间介16，15(1) 求回归方程设x1,x2xp是确定变量，Y是随机变量，他们之间有关系Y=+ 其中。，为常数，是随机变量，且N（0，），这是P元线性回归模型，我们讨论P1的情形。作n次独立试验，得到n组数据 （1+,(k=1,2,n)记j=kj j=1,2,p,则（1）式可写为 Y=（x1-1）+p)+

2、其中 同于（1）中之，而=1+p 对上面得到的几组试验数据，便有 Y=-)+(-)+其中独立分布：N（0，）。 为了用最小乘法求（2）中，的估计值，我们引如下述符号 = Y= A= =A为准对角阵，子块L是P介实对称可逆阵。， 即A为L=B=由此得正规方程组=利用分块乘法得n及l=从而得估计计算公式:u=， = 经验回归方程为 （2）假设检验线性回归的显著性检验.在线性模型(2)中作假设 He2 b1=0,b2=0,bp=0 由 QI= 利正规方程组,可知右边第三项为0,从而Qi=记 Qe= Q2=j其中Bj(j=1,2,p)是正规方程组的右端项.称为剩余离差(平方和),是由试验或引起的误差;

3、(平方和),是由线性回归引起.由分解定理知: 服从自由度为n-p-1的分布 服从自由度为p的分布记 = ,=则 F=F(p,n-p-1)由该式计算出F值,对给出的显著性水平,F(p,n-p-1), F=F(p,n-p-1)则拒H0,即认为线性回归显著;则接受H0,认为线性回归不显著.用复相关系数R,且R=,也可说明Y与自变量的密切程度,01, 越大越密切回归系数的显著性检验若经检验线性回归显著，即说明回归系数b1,b2,bp不全为0,但是不能说明每个自变量对Y都是重要的，如果某个系数为0或接近于0相应的自变量对Y不起作用或作用很小，可以忽略，因而检验每个回归系数bj(1=j=p)是否为0，相当

4、于检验相应的Xj对Y的值是否起作用。在线性回归模型上作假设 H0：bj=0，其中j固定，1=j= 则拒绝假设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可用下式计算记则(N(0,d)从而U=N(0,1)由此知T=给出置信概率1-a，那么

5、Y0的置信区间为（ +ta/2（n-p-1）式中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、

6、 输入参数（1） 对1中输入的参数全部打印输出（2） 计算回归方程的一些中间运算结果及最后结果。X0(P1)一维实数组，存放各变量的均值L（p,p）正规方程组的系数矩阵LY（p）正规方程的右端项C（p,p）L（p,p）的逆矩阵B（p）回归系数回归方程（3） 检验线性回归显著性的一些中间运算结果和检验结果方差分析表检验结果RR复相关系数，0RR=FR THEN 510 ELSE 520 510 LPRINT “线形回归显著”:GOTO 530 520 LPRINT “线形回归不显著” 530 IF F=TR THEN 620 ELSE 630 620 LPRINT “B(“;J;”)与0有显著差

7、异 “:GOTO 640 630 LPRINT “B(“;J;”)与0无显著差异” 640 NEXT J 650 LPRINT 660 * 670 INPUT “WW=”:WW 680 IF WW=0 THEN 690 ELSE 700 690 LPRINT TAB(3);”WW=”WW:GOTO 1000 700 FOR J=1 TO P 710 INPUT “X(J)=”;X(J) 720 NEXT J 730 GOSUB 1980 740 LPRINT :LPRINT TAB(3);”预报：X(J):” 750 FOR J=1 TO P 760 LPRINT X(J);:LPRINT;

8、770 IF (J/5)*5=J THEN LPRINT 780 NEXT J 790 LPRINT 800 LPRINT TAB(12);”Y0=”;Y0 810 LPRINT TAB(3);Y1;”Y0”;Y2 820 LPRINT :GOTO 670 830 DATA 2,18,50,4.3302,7,9,40,3.6485,5,14,46,4.4830 840 DATA 12,3,43,5.5468,1,20,64,5.4970,3,12,40,3.1125 850 DATA 3,17,64,5.1182,6,5,39,3.8759,7,8,37,4.6700 860 DATA 0,2

9、3,55,4.9536,3,16,60,5.0060,0,18,40,5.2701 870 DATA 8,4,50,5.3772,6,14,51,5.4849,0,21,51,4.5960 880 DATA 3,14,51,5.6645,7,12,56,6.0795,16,0,48,3.2194 890 DATA 6,16,45,5.8076,0,15,52,4.7308,9,0,40,4.6805 900 DATA 4,6,32,3.1272,0,17,47,3.6104,9,0,44,3.7174 910 DATA 2,16,39,3.8946,9,639,2.7066,12,5,51,5

10、.6314 920 DATA 6,13,41,5.8152,12,7,47,5.1302,0,24,61,5.3910 930 DATA 5,12,37,4.4533,4,15,49,4.6569,0,20,45,4.5212 940DATA 6,16,42,4.8650,4,17,48,5.3566,10,4,48,4.6098 950 DATA 4,14,36,2.3815,5,13,36,3.8746,9,8,51,4.5919 960 DATA 6,1354,5.1588,5,8,100,5.4373,5,11,44,3.9960 970 DATA 8,6,63,4.3970,2,13

11、,50,4.0622,7,8,50,2.2905 980DATA 4,10,45,4.7115,10,5,40,4.5310,3,17,64,5.3637 990 DATA 4,155,72,6.0771 1000 END 1100 子程序（1） 1110 FOR J=1 TO P1 1120 D=0 1130 FOR I=1 TO N 1140 D=D+A(I,J) 1150 NEXT I 1160 X0(J)=D/N 1170 NEXT J 1180 LPRINT TAB(3);”各变量匀值：” 1190 J=1 TO P1 1200 LPRINT X0(J);:LPRINT” “; 12

12、10 NEXT J 1220 LPRINT 1230 * 1240 FOR I=1 TO P 1250 FOR J=1 TO P 1260 D=0 1270 FOR K=1 TO N 1280 D= D+(A(K,I)-X0(I)*(A(K,J)-X0(J) 1290 NEXT K 1300 L(I,J)=D 1310 NEXT J 1320 NEXT J 1330 FOR I=1 TO P 1340 D=0 1350 FOR K=1 TO N 1360 D =D+(A(K,I)-X0(I)*(A(K,P1)-X0(P1) 1370 NEXT K 1380 LY(I)=D 1390 NEXT

13、I 1400 LPRINT 1410 LPRINT TAB(3);”正规方程组的系数矩阵” 1420 FOR I=1 TO P 1430 FOR J=1 TO P 1440 LPRINT USING #.#”:LY (I);:LPRINT” “; 1450 NEXT J 1460 LPRINT 1470 NEXT I 1480 LPRINT 1490 LPRINT TAB(3):”正规方程组的右端项” 1500 FOR I=1 TO P 1510 LPRINT USING#.#”LY(I);:LPRINT” ”; 1520 NEXT I 1530 LPRINT 1540 * 1550 FOR

14、I=1 TO P 1560 FOR J=1 TO P 1570 S(I,J)=L(I,J) 1580 NEXT J 1590 S(I,P1)=LY(I) 1600 NEXT I 1610 GOSUB 2130 1620 LPRINT 1630 LPRINT TAB(3);”系数矩阵L(P,P)的逆阵” 1640 FIR I=1 TO P 1650 FOR J=1 TO P 1660 LPRINT USING”#.#”;C(I,J):;LPRINT” “ 1670 NEXT J 1680 LPRINT 1690 NEXT I 1691*1692 BO=XO(P1)1693 FOR J=1 TO

15、P1694 B0=B0-B(J)*X0(J)1695 NEXT J1696 LPRINT1697 RETURN1698 子程序21699 Q1=0:QT=01700 FOR K=1 TO P1701 Q1=Q1/B(K)*LY(K)1702 NEXT K1703 FOR I=1 TO N1704 QT=QT/(A(I,P1)-X0(P1)21705 NEXT I1706 Q2=QT-Q11707 “*1708 P2=N-P11709 S1=Q1/P:S2=Q2/P21710 *1711 F=S1/S21712 RR=SQR(Q1/QT)1713 RETURN1714 子程序31715 FOR

16、J=1 TO P1716 T(J)=B(J)/(SQR(C(J,J)*S2)1717 NEXT J1718 子程序41719 Y0=0:H=01720 FOR J=1 TO P1721 Y0=Y0/B(J)*X(J)1722 NEXT J1723 Y0=Y0+B01724 *1725 FOR I=1 TO P1726 FOR J=1 TO P1727 H=H+C(I,J)*(X(I)-X0(I)*(X(J)-X0(J)1728 NEXT J,I1729 D0=SQR(1+1/N+H)1730 YY=TR*D0*SQR(S2)1731 Y1=Y0-YY:Y2=Y0+YY1732 RETURN17

17、33 子程序51734 FOR K=1 TO P1735 A=1/S(K,K)1736 FOR I=1 TO P1737 FOR J=1 TO P11738 IF IK AND JK THEN S(I,J)*-S(I,J)-S(I,K)*S(K,J)*A1739 NEXT J1740 NEXT I1741 FOR J=1 TO P11742 S(K,J)=S(K,J)*A1743 IF JP1 THEN S(,K)=-S(J,K)*A1744 NEXT J1745 FOR I=1 TO P1746 FOR J=1TO P1747 C(I,J)=S(I,J)1748 NEXT J1749 NEX

18、T I1750 *1751 FOR I=1 TO P1752 B(I)=S(I,P)1753 NEXT I1754 RETURN例:P=3 N=49序号 X1 X2 X3 Y1 2 18 50 4.33022 7 9 40 3.64853 5 14 46 4.4834 12 3 43 5.54685 1 20 64 5.11286 3 12 40 3.87597 3 17 64 5.11288 6 5 39 3.87599 7 8 37 4.6710 0 23 55 4.953611 3 16 60 5.00612 0 18 40 5.271013 8 4 50 5.377214 6 14 5

19、1 5.484915 0 21 51 4.59616 3 14 51 5.664517 7 12 56 6.079518 16 0 48 3.219419 6 16 45 5.807620 0 15 52 4.730621 9 0 40 4.680522 4 6 32 3.127223 0 17 47 3.610424 9 0 44 3.717425 2 16 39 3.894626 9 6 39 2.706627 12 5 51 5.631428 6 13 41 5.815229 12 7 47 5.130230 0 24 61 5.39131 5 12 37 4.453332 4 15 49 4.556933 0 20 45 4.521234 6 16 42 4.86535 4 17 48 5.256336 10 4 48 4.6008537 4 14 50 2.381538 5 12 50 3.37483