Bootstrap法による半減期の信頼区間の推定と検定

先ほどの記事では半減期の推定を行いましたが、統計的なばらつきを考慮するためにBootstrap法という手法を使って信頼区間を推定し検定します。

前回に引き続き仙台市青葉区のデータを解析しています。


Bootstrap法は分布を仮定しなくてもノンパラメトリックに推定値の分布を推定できる優れた方法です。

簡単に説明しますと、使っているデータをリサンプリングすることによって擬似的にデータを作り、そのデータでパラメータの推定をします。

これを何回か繰り返してパラメータの分布を得るのです（推奨されているのは2,000回）。


前回の半減期を求める関数を使ってBootstrap法でデータをリサンプリングして、半減期を2,000回求めるとこのような分布になりました。


f:id:isseing333:20110331093836j:image


このヒストグラムのパーセンタイルから、95%信頼区間と中央値を求めると次のようになります。


95%下側限界中央値95%上側限界
（2.5パーセンタイル）（50パーセンタイル）（97.5パーセンタイル）
6.04日7.11日9.26日


中央値の7.11日は前回のブログで推定した7.09日とほぼ一致しているので、Bootstrap法のバイアスは無いようです。

この表から95%信頼区間は約6日から約9.3日となります。

I133の半減期である約8日はこの中に入っています。


ということは、この8日を帰無仮説にしてBootstrap法によって検定すると棄却されない。

つまり「統計的には放射物質がI133であることが示唆される」という事になります。



もし棄却された場合は次のようなものが考えられます。

  1. I133以外の物質が混ざっている
  2. 半減期を早くする理由が別にある

1の場合はモデルを少し修正して、物質が混合しているモデルにして推定する必要があります。

2の場合は風による拡散などの影響でしょう。これも風の効果などをモデル化することは可能です。


今回の検定の結果では、1も2も示唆されてません（もちろん最初のモデルが妥当であれば、です）。

繰り返しになりますが私は物理屋さんではありませんので、今回の記事は「統計学をどうやって使えばいいのか」という1つの見本になればと思ってます。



ブートストラップ法はとても汎用性の高い方法ですので、使えるようになると重宝します。


ブートストラップ法による信頼区間の構成を応用したABC法という方法がさらに精度が高かったはずです。

これはまだプログラムにしていませんので、また更新致します。


※ベースライン補正の方法が適切でないとのご指摘を受け以下のブログに修正版の記事を書きました。

http://d.hatena.ne.jp/isseing333/20110331/1301575593


HoushaSen <- read.csv("放射線.csv", as.is = T)


OptimHangenki <- function(Tvec, Hvec){
	T <- Tvec
	H <- Hvec
	Hangenki <- function(Beta){
		a  <- -1.1
		T1 <- 14
		T2 <- 24
		a1 <- Beta[1]
		a2 <- Beta[2]
		b1 <- Beta[3]
		t  <- T
		Pred  <- ifelse(t  <= T1, a, 
			 ifelse(t  <= T2, a + a1 + b1*(t - T1),
			 ifelse(T2 <  t, a + a1 + a2 + b1*(t - T1), NA)))
		Resid <- log10(H) - Pred
		Resid%*%Resid
	}
	optim(c(0, 0, 0), Hangenki)
}


BootHangenki <- function(nboot, Housha, T1){
	#------bootstrap
	library(bootstrap)
	
	#---make data frame using T and D
	Data  <- data.frame(Housha, T1)
	
	#---produce NULL data
	ResultBoot      <- matrix(NA, nboot, 3)
	
	#---perform bootstrap SSI
	for(i in 1:nboot){
		set.seed(i)
		BootID          <- sample(1:nrow(Data), nrow(Data), replace = "T")
		BootSample      <- Data[BootID, ]
		BootSample      <- BootSample[order(as.numeric(rownames(BootSample))), ]
		
		T <- BootSample[, 2]
		H <- BootSample[, 1]
		PredHangenki    <- OptimHangenki(T, H)
		ResultBoot[i, ] <- PredHangenki$par
		print(i)
	}
	
	a1 <- ResultBoot[, 1]
	a2 <- ResultBoot[, 2]
	b1 <- ResultBoot[, 3]
	hangen <- -(log10(20) - log10(10)) / b1
	
	print("median, 95% CI of Bootstrap hangenki")
	print(round(quantile(hangen, probs = c(0.025, 0.5, 0.975), na.rm = T), digits =2))
	
	ResultBoot
}

Boot <- BootHangenki(2000, HoushaSen[, 2], HoushaSen[, 3])

library(ggplot2)
Temp <- data.frame(半減期 = -(log10(20) - log10(10)) / Boot[, 3])
ggplot(Temp, aes(半減期)) + geom_histogram()

ページTOPへ