簡単なBootstrap信頼区間の計算（bootstrap confidence interval）

分布が良く分からない統計量の信頼区間を求める場合は、Bootstrap sampleを使って信頼区間を構成するのが簡単です。

Efron先生が論文で良く例に出すのが、主成分分析の固有値の信頼区間。

確かに何の分布に従ってるのか良く分からないですね。


パッケージはbootstrapがあります。

プログラム例はこんな感じ↓

library(bootstrao)
x <- rnorm(20)
theta <- function(x){mean(x)}
results <- bootstrap(x,100,theta)

でもこれだとパラメータが1次元なので、汎用性は低いです（うまいやり方があるのかもしれないけど）。


そこで、sample関数を使ってbootstrap sampleを作る方法を紹介します。

プログラムはこちら↓

#---make null data
ResultBoot      <- matrix(NA, nboot, npar)

#---perform bootstrap sample
for(i in 1:nboot){
	BootID          <- sample(1:nrow(Data), nrow(Data), replace = "T")
	#---make bootstrap sample
	BootSample      <- Data[BootID, ]
	#---calculate bootstrap "star" estimate
	ResultBoot[i, ] <- SomeFunction(BootSample, parameter)
}

(MeanBoot <- mean(ResultBoot[, 1]))
(SDBoot   <- sd(ResultBoot[, 1]))

#---95% Bootstrap Confidence Interval（95% ブートストラップ信頼区間）
MeanBoot - 1.96*SDBoot
MeanBoot + 1.96*SDBoot

残念ながらコピペで動くプログラムじゃありませんが、適切にパラメータ設定できればどんな関数で推定したものもブートストラップ信頼区間を作れます。

いじる必要のあるパラメータは、

で、信頼区間を付けたい統計量が1つしかない場合は、ResultBoot[i, ]をResultBoot[i]にすれば問題ないはず（未確認汗）。

こんなわけでブートストラップ信頼区間も簡単に作れます。

是非お試し下さい。

ページTOPへ