回归系数的显著性检验

907求多元线形回归方程及预报
一．功能 x1，x2, x3,………..xp 为自变量，Y 为随机变量， 求线形回归方程
Y=b 0+x b 11+…+x b p p +ε
其中,,10b b 。

，b p 为常数，ε是随机变量，且ε—N （0，26） 来描述Y 与X 的变化规律。

并用T 检验法检验线形回归是否显著。

如果线性回归显著，可用经验回归平面方程对Y 作出预报，并给出预报值的置信区间。

二．算法间介[16]，[15] (1) 求回归方程
设x1,x2…xp 是确定变量，Y 是随机变量，他们之间有关系 Y=b 0+x b 11+…+x b p p +ε
其中,,10b b 。

，b p 为常数，ε是随机变量，且ε—N （0，26），这是P 元线性回归模型， 我们讨论P>1的情形。

作n 次独立试验，得到n 组数据 （x k 1+2x k +,…),Y x p k p (k=1,2,…,n)
记X j=∑=n
k X n 1
1kj j=1,2,…,p,则（1）式可写为
Y=+μb 1（x1-X 1）+…+(b p -x
p
X
p)+ε
其中 ε同于（1）中之ε，而μ=X b b 10+1+…+X b p p 对上面得到的几组试验数据，便有
Y=X b k 11(+μ-X 1)+…+b p (X kp -X p )+εk 其中ε独立分布：εk —N （0，26）。

为了用最小乘法求（2）中μ，b b b p ,...,,21的估计值，我们引如下述符号
Y =n
1∑=n
i Y
1
1
Y=⎥
⎥
⎥
⎥⎥
⎥⎦⎤⎢⎢⎢⎢⎢
⎢⎣
⎡--
-
----
--X X X
X
X
X
X X X X X X X
X
X X X X p np n n P P P
P
......,
1...............
.......,1
....
,1
.2
2
1122
221212
12
1
11
⎪⎪⎪
⎩
⎪
⎪⎪⎨⎧
=-=--=∑∑==,,....,2,1,,
)(),
()(`11
p k j y X X L X X X X L i n i j ij jy k ik n
i j ij jk
A==X X t
⎥⎥⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢
⎢⎢⎣
⎡---
---∑∑∑∑i p p i
p
p p i
i X X X X X X
X X X X X X
n 211111
111112
1
11
)(...)()
(0.......................
............))(....)(00
0 0
(
=⎥⎥⎥⎥⎥⎥
⎦
⎤⎢⎢⎢⎢⎢⎢⎣
⎡L L
L
L L L L L L pp p P p p n ....
0...........
....
0
...00 0
2
1
222
21
11211=⎥⎦⎤
⎢⎣⎡L n 00 A 为准对角阵，子块L 是P 介实对称可逆阵。

，
即A 为L=⎥
⎥
⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡PP P P P P L L L L L L
L L L 2
1
22221
11211
B=Y X T ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--∑∑∑)()(111P P I I I I
X X Y X X Y Y =⎥⎥⎥⎥⎦
⎤
⎢⎢⎢⎢⎣⎡∑LPY Y L Y I 1
由此得正规方程组
⎥⎦⎤⎢⎣⎡l o n 0=⎥⎥⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡^^2^1^p b b b u =⎥⎥⎥
⎥⎥⎥
⎦⎤⎢⎢⎢⎢⎢
⎢⎣⎡
∑py y y i l l l y 21 利用分块乘法得 n ∑=i
i y u ^
及
l ⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎣⎡^^2^1p b b b =⎥⎥⎥⎥⎥
⎦⎤
⎢⎢⎢⎢⎢⎣⎡py y y l l l 21 从而得估计计算公式:
u=^y ，⎥⎥⎥⎥⎥⎥
⎥⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡^^^21...b b b p =1-L ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤
⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣
⎡L L L py y y ...21 经验回归方程为
^
^
1^
1^
,....μ+++=x b x b p p Y （2）假设检验
线性回归的显著性检验.在线性模型(2)中作假设 He2 b1=0,b2=0,…bp=0 由
QI=
=-+-∑
=2)]ˆ(2)ˆ[(1
Y I Y I Y YI n
i ∑∑∑===--+-+-n
i n
i n
i Y YI I Y YI Y I Y I Y YI 1
1
1
))(ˆ(22)ˆ(2)ˆ( 利正规方程组,可知右边第三项为0,从而 Qi=∑∑==-+-n
i n
i Y I Y I Y YI 1
1
2)ˆ(2)ˆ(
记 Qe=∑=-n
i i Y Yi 1
2)ˆ(
Q2=∑∑===-n i p
j B b Y i
Y 1
ˆ2)ˆ(j 其中Bj(j=1,2,…,p)是正规方程组的右端项.Q E 称为剩余离差(平方和),是由试验或
σ
2
引起的误差;称为回归离差Q 2(平方和),是由线性回归引起.
由分解定理知:
σ
21
Q
E
服从自由度为n-p-1的X 2
分布 σ
2
1
Q 1
服从自由度为p 的X 2分布
记
S E 2
=
1
--p n Q
E
,
S 2
2=
p
Q
2
则 F=S
S E
22
2
~F(p,n-p-1)
由该式计算出F 值,对给出的显著性水平α,F α(p,n-p-1),
F>=F α(p,n-p-1)
则拒H0,即认为线性回归显著;则接受H0,认为线性回归不显著.用复相关系数R,且R 2=1
2
Q Q ,也可说明Y 与自变量的密切程度,0<2R <1, 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 （
-
）=0
D （
-
）=
+
+
-
）（
其中Cij 是C=A-1第I 行第j 列元素，62可用下式计算
记))((1101
011
1
20
j j p
j i i p
i X X X X C
n d --+=∑∑
==
则()^
00Y Y -~N(0,2σd 2
0)
从而U=
^
0d Y Y σ-~N(0,1) 由此知
T=
)1(~0
^
0---p n d Y Y σ
给出置信概率1-a ，那么Y0的置信区间为
（
+ta/2（n-p-1）0d
式中da 可用（18）式计算后经开方得到。

三、程序说明
程序中分主程序与子程序，子程序有五个。

子程序（1）求线性回归方程 子程序（2）线性回归的显著性检验 子程序（3）回归系数的显著检验 子程序（4）预报
子程序（5）计算正规方程组系数矩阵之逆阵
注意，当检验出线性回归不显著时，就不再执行子程序（3）和（4） 四、程序使用说明 1、 输入参数 p 实数，自变量个数 N 实数，实验次数
R 实数，给定的显著性水平，可取R=0.05
FR 实数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）正规方程组的系数矩阵
L Y（p）正规方程的右端项
C（p,p）L（p,p）的逆矩阵
B（p）回归系数
回归方程
（3）检验线性回归显著性的一些中间运算结果和检验结果方差分析表
检验结果
RR 复相关系数，0<RR<1
（4）检验回归系数显著的结果
B1与0有显著差别
（5）预报结果
3、中间工作单元略
10********************
20 907 多元线莘回归及预报*
30 *******************
40 INPUT”P,N=”P,N:P1=P+1
50 DIM A(N,P1),X)(P1),L(P,P),C(P,),L Y(P),S(P,P1) 60 DIM B(P)
70 FOR I=1 TO N
80 FOR J=1 TO P1
90 READ A(I,J)
100 NEXT J,I
110 LPRINT TAB(3);”例：“；“P=”；P；“N=”；N：LPRINT
120 LPRINT“序号”；TAB（9）；“X1”；TAB（17）；“X2”；TAB（25）；“X3”；
130 LPRINT TAB（35）；“Y”
140 FOR I=1 TO P1
150 LPRINT I；
160FOR J=1 TO P1
170LPRINT TAB（8*J）；A（I，J）；
180NEXT J
190LRINT
200 NEXT I
210 LRINT：LPRINT
220 GOSUB 1：100
230 LPRINT
240 LPRINT TAB（3）：“回归系数和回归方程”
250 FOR J=1 TO P
260 LPRINT“B（”；J；“）=”；USING“#。

#####”；B（J）；270 IF INT（I/4）*4=I THEN LPRINT
280NEXT J
290 LPRINT
300 LPRINT“Y=”；USING“#。

#####”；B0；
310 FOR I-1 TO P
320IF B（1）》0 THEN LPRINT“+”；
330LPRINT USING“#。

####”；B（1）；
340 LPRINT“*X”；I；
350 IF INT（I/3）*3=I THEN LPRINT
360 NEXT I
370 LPRINT
380 ****
390 INPUT“R，Fr，Tr=”；R，FR，TR
400 GOSUB 1770
410 LPRINT
420 LPRINT TAB(3);”R=”;R;”Fr=”FR;”Tr=”TR:LPRINT
430 LPRINT TAB(15);”方差分析表：“
440 LPRINT “来源“;TAB(8);”离差”;TAB(16);”自由度”；
450 LPRINT TAB(26);”均方离差”;TAB(38);”F值”
460 LPRINT “Q1”;TAB(5);Q1;TAB(18);P;TAB(25);S1;TAB(35);F 470 LPRINT “Q2”;TAB(5);Q2;TAB(18);N-P1;TAB(25);S2
480 LPRI NT “Qt”;TAB(5);QT;TAB(18);N-1
490 LPRINT:LPRINT TAB(3);”RR=”;RR:LPRINT
500 IF F>=FR THEN 510 ELSE 520
510 LPRINT “线形回归显著”:GOTO 530
520 LPRINT “线形回归不显著”
530 IF F<FR THEN 1000
540 GOSUB 1930
550 LPRINT:LPRINT TAB(3);”T(J):”
560 FOR J=1 TO P
570 LPRI NT T(J);:LPRINT” “;
580 NEXT J
590 LPRINT:LPRINT
600 FOR J=1TO P
610 IF ABS(T(J))>=TR THEN 620 ELSE 630
620 LPRINT “B(“;J;”)与0有显著差异“: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;
770 IF (J/5)*5=J THEN LPRINT
780 NEXT J
790 LPRINT
800 L PRINT TAB(12);”Y0=”;Y0
810 LPRINT TAB(3);Y1;”<Y0<”;Y2
820 LPRINT :GOTO 670
830 DA TA 2,18,50,4.3302,7,9,40,3.6485,5,14,46,4.4830 840 DA TA 12,3,43,5.5468,1,20,64,5.4970,3,12,40,3.1125
850 DA TA 3,17,64,5.1182,6,5,39,3.8759,7,8,37,4.6700 860 DA TA 0,23,55,4.9536,3,16,60,5.0060,0,18,40,5.2701 870 DA TA 8,4,50,5.3772,6,14,51,5.4849,0,21,51,4.5960 880 DA TA 3,14,51,5.6645,7,12,56,6.0795,16,0,48,3.2194 890 DA TA 6,16,45,5.8076,0,15,52,4.7308,9,0,40,4.6805 900 DA TA 4,6,32,3.1272,0,17,47,3.6104,9,0,44,3.7174 910 DA TA 2,16,,39,3.8946,9,639,2.7066,12,5,51,5.6314 920 DA TA 6,13,41,5.8152,12,7,47,5.1302,0,24,61,5.3910 930 DA TA 5,12,37,4.4533,4,15,49,,4.6569,0,20,45,4.5212 940DA TA 6,16,42,4.8650,4,17,48,5.3566,10,4,48,4.6098 950 DA TA 4,14,36,2.3815,5,13,36,3.8746,,9,8,51,4.5919 960 DA TA 6,1354,5.1588,5,8,100,5.4373,5,11,44,3.9960 970 DA TA 8,6,63,4.3970,2,13,50,4.0622,7,8,50,2.2905 980DA TA 4,10,45,4.7115,10,5,40,4.5310,3,17,64,5.3637 990 DA TA 4,155,72,6.0771
1000 END
1100 子程序（1）
1110 FOR J=1 TO P1
1120 D=0
1130FOR I=1 TO N
1140D=D+A(I,J)
1150NEXT I
1160X0(J)=D/N
1170 NEXT J
1180 LPRINT TAB(3);”各变量匀值：” 1190 J=1 TO P1
1200LPRINT X0(J);:LPRINT” “; 1210 NEXT J
1220 LPRINT
1230 ****
1240 FOR I=1 TO P
1260D=0
1270FOR K=1 TO N
1280D= D+(A(K,I)-X0(I))*(A(K,J)-X0(J)) 1290 NEXT K
1300L(I,J)=D
1310NEXT J
1320 NEXT J
1330 FOR I=1 TO P
1340 D=0
1350FOR K=1 TO N
1360 D =D+(A(K,I)-X0(I))*(A(K,P1)-X0(P1)) 1370NEXT K
1380LY(I)=D
1390 NEXT I
1400 LPRINT
1410 LPRINT TAB(3);”正规方程组的系数矩阵” 1420 FOR I=1 TO P
1440LPRINT USING #####.###”:LY (I);:LPRINT” “; 1450NEXT J
1460LPRINT
1470 NEXT I
1480 LPRINT
1490 LPRINT TAB(3):”正规方程组的右端项”
1500FOR I=1 TO P
1510 LPRINT USING#####.###”L Y(I);:LPRINT” ”; 1520 NEXT I
1530 LPRINT
1540 ****
1550 FOR I=1 TO P
1560FOR J=1 TO P
1570S(I,J)=L(I,J)
1580NEXT J
1590S(I,P1)=L Y(I)
1600 NEXT I
1610 GOSUB 2130
1620 LPRINT
1630 LPRINT TAB(3);”系数矩阵L(P,P)的逆阵” 1640 FIR I=1 TO P
1650FOR J=1 TO P
1660LPRINT USING”##.#####”;C(I,J):;LPRINT” “ 1670NEXT J
1680LPRINT
1690 NEXT I
1691…****
1692 BO=XO(P1)
1693 FOR J=1 TO P
1694 B0=B0-B(J)*X0(J)
1695 NEXT J
1696 LPRINT
1697 RETURN
1698 …子程序2
1699 Q1=0:QT=0
1700 FOR K=1 TO P
1701 Q1=Q1/B(K)*L Y(K)
1702 NEXT K
1703 FOR I=1 TO N
1704 QT=QT/(A(I,P1)-X0(P1))^2 1705 NEXT I
1706 Q2=QT-Q1
1707 “****
1708 P2=N-P1
1709 S1=Q1/P:S2=Q2/P2
1710 …***
1711 F=S1/S2
1712 RR=SQR(Q1/QT)
1713 RETURN
1714 …子程序3
1715 FOR J=1 TO P
1716 T(J)=B(J)/(SQR(C(J,J)*S2)) 1717 NEXT J
1718 …子程序4
1719 Y0=0:H=0
1720 FOR J=1 TO P
1721 Y0=Y0/B(J)*X(J)
1722 NEXT J
1723 Y0=Y0+B0
1724 …****
1725 FOR I=1 TO P
1726 FOR J=1 TO P
1727 H=H+C(I,J)*(X(I)-X0(I))*(X(J)-X0(J))
1728 NEXT J,I
1729 D0=SQR(1+1/N+H)
1730 YY=TR*D0*SQR(S2)
1731 Y1=Y0-YY:Y2=Y0+YY
1732 RETURN
1733 …子程序5
1734 FOR K=1 TO P
1735 A=1/S(K,K)
1736 FOR I=1 TO P
1737 FOR J=1 TO P1
1738 IF I<>K AND J<>K THEN S(I,J)*-S(I,J)-S(I,K)*S(K,J)*A 1739 NEXT J
1740 NEXT I
1741 FOR J=1 TO P1
1742 S(K,J)=S(K,J)*A
1743 IF J<>P1 THEN S(,K)=-S(J,K)*A
1744 NEXT J
1745 FOR I=1 TO P
1746 FOR J=1TO P
1747 C(I,J)=S(I,J)
1748 NEXT J
1749 NEXT I
1750 …****
1751 FOR I=1 TO P
1752 B(I)=S(I,P)
1753 NEXT I
1754 RETURN
例:P=3 N=49
序号X1 X2 X3 Y
1 2 18 50 4.3302
2 7 9 40 3.6485
3 5 1
4 46 4.483
4 12 3 43 5.5468
5 1 20 64 5.1128
6 3 12 40 3.8759
7 3 17 64 5.1128
8 6 5 39 3.8759
9 7 8 37 4.67
10 0 23 55 4.9536
11 3 16 60 5.006
12 0 18 40 5.2710
13 8 4 50 5.3772
14 6 14 51 5.4849
15 0 21 51 4.596
16 3 14 51 5.6645
17 7 12 56 6.0795
18 16 0 48 3.2194
19 6 16 45 5.8076
20 0 15 52 4.7306
21 9 0 40 4.6805
22 4 6 32 3.1272
23 0 17 47 3.6104
24 9 0 44 3.7174
25 2 16 39 3.8946
26 9 6 39 2.7066
27 12 5 51 5.6314
28 6 13 41 5.8152
29 12 7 47 5.1302
30 0 24 61 5.391
31 5 12 37 4.4533
32 4 15 49 4.5569
33 0 20 45 4.5212
34 6 16 42 4.865
35 4 17 48 5.2563
36 10 4 48 4.60085
37 4 14 50 2.3815
38 5 12 50 3.3748
39 9 8 51 4.5919
40 6 13 54 5.1538
41 5 8 100 5.4373
42 5 11 44 3.990
43 8 6 63 4.397
44 2 13 50 4.8622
45 7 8 50 2.2805
46 4 10 45 4.7115
47 10 5 40 4.531
48 3 17 64 5.3637
49 4 15 72 6.0771 各变量均值:
5.285714 11.79592 48.91873 4.002284 662.000 -918.143 -324.857
918.143 1753.960 714.184
324.857 714.184 6295.674
正规方程组的系数矩阵
-11.717 74.334 240.118
正规方程组的右端顶
0.005525 0.002911 -0.000045
0.002911 0.002131 -0.000092
0.000045 -0.000092 0.000157
回归系数和回归方程
B(1)=0.140798
B(2)=0.102322
B(3)=0.033798
Y=0.99775/0.1408*X1/0.1023*X2/0.0338*X3
R=.05 FR=2.84 TR=2.0141
方差分析表
来源离差自由度均方离差F值
Q1 14.07168 3 4.690559 7.573205 Q2 27.87098 45 .6193551
QT 41.94265 48
RR=.5792219
线性回归显著
T(J):
2.406815 2.816406
3.324249
B(1)与0有显著差异
B(2)
B(3)
预报:X(J);
12
20
60
Y0=6.761631
4.676331<Y0<8.846931
预报:X(J):
15
55]
80
Y0=11.44124
6.999876<Y0<15.8826
预报:X(J):
3
15
55
Y0=4.813854
3.203319<U0<6.424349
VW=0。