仙台市の放射線値の解析
この方が行っている解析を参考にして、宮城県のホームページからデータを取って放射線値の推移を解析しました。
数値が17日分しかなかったので、1つずつページを開いてcsvに記録しました。
今回は仙台市青葉区のデータです。
まずは例によってスプラインの結果です。
横軸が日で、縦軸が自然対数のマイクロシーベルト/時の値です。
この解析の目的は、放射線を発している物質(?)が何なのかつきとめること(?)なので、半減期を推定する必要があります。
(物理が専門ではないので?が多くて申し訳ありません)
なのでスプラインで傾向を掴むだけでは役に立ちません。
そこで、この対数値に対して次のような簡単なモデルを当てはめます。
- if t <= T1: 観測値 = a
- if t <= T2: 観測値 = a + a1 + b1*(t - T1)
- if T2 < t : 観測値 = a + a1 + a2 + b1*(t - T1)
a、a1、a2はベースラインの飛散値と何かによってまき散らされて増えた飛散値。
b1は減衰速度(log10(0.5)/b1が半減期)。
T1、T2はまき散らされたタイミングです(推定パラメータを減らすため、今回の解析では14日と24日に決め打ちしました)。
このT1とT2の次の日に飛散量が増えるというモデルになっています。
このモデルを適応して最小二乗法でパラメータを推定し、グラフにすると次のようになりました。
半減期は7.09日と推定されました。
参考にした記事では5.23日と推定されていたのでそれよりも長いですが、I131の半減期は約8日らしいのでそれなりの値に思えます。
またこのモデルより、15日と25日に増加した放射線量は0.255と0.0738マイクロシーベルト/時と推定されました。
物理が専門ではないので、物理学的にはこの結果の良し悪しは判断できないのですが、、、
このようにRを使えば適当な式を立ててoptim関数で最適化することで、様々なパラメータを推定することが出来ます。
例によって利用したデータはダウンロードのページにアップしておきます。
※ベースライン補正の方法が適切でないとのご指摘を受け以下のブログに修正版の記事を書きました。
http://d.hatena.ne.jp/isseing333/20110331/1301575593
HoushaSen <- read.csv("放射線.csv", as.is = T) plot(HoushaSen[, 3], HoushaSen[, 2], lwd = 2) plot(HoushaSen[, 3], log10(HoushaSen[, 2]), lwd = 2) Spline <- smooth.spline(HoushaSen[, 3], HoushaSen[, 2]) PredSpline <- predict(Spline, seq(10, 30, by = 0.1)) plot(HoushaSen[, 3], log10(HoushaSen[, 2]), lwd = 2, xlab = "日", ylab = "log10(放射線量)", main = "スプライン") lines(PredSpline$x, log10(PredSpline$y), lwd = 2) Hangenki <- function(Beta){ a <- -1.1 T1 <- 14 T2 <- 24 a1 <- Beta[1] a2 <- Beta[2] b1 <- Beta[3] t <- HoushaSen[, 3] 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(HoushaSen[, 2]) - Pred Resid%*%Resid } (PredHangenki <- optim(c(0, 0, 0), Hangenki)) a <- -1.1 T1 <- 14 T2 <- 24 Beta <- PredHangenki$par a1 <- Beta[1] a2 <- Beta[2] b1 <- Beta[3] t <- HoushaSen[, 3] 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(HoushaSen[, 2]) - Pred Resid%*%Resid Half <- -(log10(20) - log10(10)) / b1 plot(HoushaSen[, 3], log10(HoushaSen[, 2]), lwd = 2, xlab = "日", ylab = "log10(放射線量)", main = "仙台市青葉区") lines(t, Pred, lwd = 2) text(25, -0.5, paste("推定半減期:", round(Half, digits = 2), "日", sep = "")) #放射線の増加量 10^(-1.1+a1)-10^(-1.1) t=25 aa1 <- a + a1 + a2 + b1*(t - T1) aa2 <- a + a1 + b1*(t - T1) 10^aa1 - 10^aa2