bootstrap 事始め

この本の第一章例3をRで。

データの定義

x<- c(5,10,15,20,25,30,35,40,45,50)
y<- c(8.6,10.4,21.6,22.8,19.5,30.0,29.5,34.3,38.7,47.4)

線形回帰

re<-lm(y~x)

その残差

resi<-summary(re)$residual

箱を作って

ncount<-100000
Res<-matrix(0,ncount,length(x))
res.boot<-matrix(0,ncount,length(x))
result<-matrix(0,ncount,2)

一応時間なんかも測ったりして

time<-proc.time()

R初心者です（ループ）

for(i in 1:ncount){
Res[i,]<-sample(resi,replace=T)
res.boot[i,]<-re$coefficients[1]*array(1,length(x))+re$coefficients[2]*x+Res[i,]
temp<-lm(res.boot[i,]~x)
result[i,1]<-temp$coefficients[1]
result[i,2]<-temp$coefficients[2]
}
時間見て
proc.time()-time

ユーザ	システム	経過(秒)
1.16	0.00	1.16

1000回で1.5秒か・・・。

summary(result[,1])

Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
-0.8874	3.6070	4.9510	4.9550	6.3630	11.0900

sd(result[,1])
[1] 2.070965

ちなみにただの線形回帰なら

re$coefficients[1]
(Intercept)
4.9

点推定値はいいけど結構散らばっている。まあサンプル10個だからな。

summary(result[,2])

Min.	1st Qu.	Median	Mean	3rd Qu.	Max.
0.6035	0.7282	0.7741	0.7757	0.8175	0.9668

sd(result[,2])
[1] 0.06594192

回帰係数は

re$coefficients[2]
x
0.7774545

こっちはほとんど動かない。
パッケージがあるかとかは調べてない。

一応これでブートストラップ。
事始めでもないか。一応回帰だし。