setwd("E:")
mydata<-read.table(file="布莱泽帕根与怀特检验.txt",header=TRUE)
#mydata
names(mydata)
#"id"    "EXP"   "POP"   "PCAID" "PCINC"
#要求用POP AID INC解释EXP
attach(mydata)
mylm<-lm(EXP~POP+PCAID+PCINC)
summary(mylm)
# b<-c()
# for(i in 1:4){
#   b[i]<-mylm$coefficients[i]
# }
# b
# EXP_predict<-b[1]+b[2]*POP+b[3]*PCAID+b[4]*PCINC
# res<-EXP-EXP_predict
# res
#或者这样
res<-as.numeric(mylm$residuals)
res
plot(id,res,col="red")
#做怀特检验
#认为残差由POP引起
lm_w<-lm(res~POP)
lm_w_sum<-summary(lm_w)
r2<-lm_w_sum$r.squared
n<-length(res)
if(n*r2>qchisq(p=1-0.001,df=1)){
  cat("在置信水平0.001下,存在异方差：","因为",
    "n*r2",n*r2,">","qchisq(p=1-0.01,df=1",qchisq(p=1-0.001,df=1))
}else{
  cat("在置信水平0.001下,不存在异方差：","因为",
    "n*r2",n*r2,"<","qchisq(p=1-0.01,df=1",qchisq(p=1-0.001,df=1))
}
#做布莱泽-帕根检验
#写出估计的sigma2
sigma2<-sum(res^2)/n
res2_sigma2<-res^2/sigma2
#认为异方差由POP引起
lm_breusch<-lm(res2_sigma2~POP)
lm_b_sum<-summary(lm_breusch)
ESS<-sum(as.numeric(lm_b_sum$residuals)^2)
TSS<-sum((res2_sigma2-mean(res2_sigma2))^2)
RSS<-TSS-ESS
if((RSS/2)>qchisq(p=1-0.001,df=1)){
  cat("在置信水平0.001下，存在异方差：","因为",
    "RSS/2",RSS/2,">","qchisq(p=1-0.001,df=1)",qchisq(p=1-0.001,df=1))
}else{
  cat("在置信水平0.001下，存在异方差：","因为",
    "RSS/2",RSS/2,"<","qchisq(p=1-0.001,df=1)",qchisq(p=1-0.001,df=1))
}

运行结果：囧
怀特检验：
在置信水平0.01下,不存在异方差： 因为 n*r2 2.364729e-31 < qchisq(p=1-0.01,df=1 6.634897
布莱泽-帕根检验：
在置信水平0.01下，存在异方差： 因为 RSS/2 91.97504 > qchisq(p=1-0.01,df=1) 6.634897
太奇怪了！！！！！


附录：原始数据是这样~~~~~~

id EXP POP PCAID PCINC
1 704 1026 186 3664
2 526 774 123 4279
3 411 460 235 3703
4 5166 5796 190 4825
5 699 969 184 4513
6 2546 3080 145 5414
7 22750 18367 240 5275
8 5911 7349 141 5379
9 8840 11905 136 4545
10 6867 10722 112 4572
11 3457 5286 103 4364
12 8935 11244 156 5162
13 7799 9013 147 4982
14 3757 4526 116 4279
15 3528 3877 163 4343
16 2108 2884 113 4316
17 3156 4747 151 4307
18 475 634 201 4128
19 521 680 195 3766
20 1052 1528 134 4451
21 1551 2268 132 4535
22 571 571 170 5222
23 3392 4048 135 5017
24 3037 4756 131 4396
25 1250 1795 252 3624
26 2938 5221 141 3868
27 1512 2688 153 3500
28 3197 4733 178 3956
29 4771 7347 114 4450
30 2063 3306 181 3634
31 2446 4072 175 3708
32 2104 3521 193 3476
33 1427 2256 255 3188
34 1014 2008 199 3345
35 2691 3738 196 3565
36 1767 2633 190 3837
37 7246 11604 141 4085
38 587 716 252 4083
39 512 755 180 3711
40 368 346 369 4269
41 1920 2364 183 4600
42 823 1076 277 3512
43 1532 1963 150 4273
44 821 1127 196 3741
45 543 533 180 5209
46 3070 3418 184 4601
47 1766 2185 201 4339
48 20052 20411 200 5087
49 698 325 570 5222
50 940 816 202 5153