東日本大震災後5年間の地震発生データ分析
Mr.Unadon(Twitter @MrUnadon)
2017/11/26
はじめに
本稿は2017年11月26日に開催された「Hijiyama.R The Final」での発表から、分析部分のみを抜き出した資料となります。
テーマは東日本大震災以降の地震発生に関するモデリングです。 ポアソン分布と二項分布の状態空間モデルを実装します。
メインの解析は二項分布の状態空間モデルとなります。
以下の構成で進めていきます。
データの取得
データの可視化と確認
メカニズムの想像
Stanモデル(ポアソン分布と二項分布の状態空間モデル)
解析結果
library(tidyverse) #データ成形、可視化まわり色々パッケージ群
library(dygraphs) #時系列データ可視化(java)
library(formattable) #html table
library(xts) #Create or Test For An xts Time-Series Object
library(rstan) #HMC ベイズ推定
library(shinystan) #収束診断等
require(RCurl)
#stan parallel
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
#set.seed
set.seed(1234)1.データの取得
東日本大震災以降、5年間のデータを取得しました。日毎のデータで、その日の地震発生回数を震度別に記録したものです。
データの取得元はこちらです。
tublaを利用してpdfからデータを抜き出してきました。csvに変換したものを利用できるようにしておきます。以下のコードでデータを成形します。
#データの読込
dat<-read.csv('https://raw.githubusercontent.com/MrUnadon/mrunadon.github.io/master/images/tabula-yukan1.csv',
fileEncoding = "utf-8",
header = F)
#最初の行だけ日付に時刻が付いているため変換 i.e. "3月14日14時"
dat$V1[1]<-"3月11日"
#列名を付与
names(dat)<-c("date_jap", #日本語表記の日時列
"Lv_1","Lv_2","Lv_3","Lv_4","Lv_5w","Lv_5s","Lv_6w","Lv_6s","Lv_7", #震度別回数
"n_over4","n_over1"#震度4以上、震度1以上(全て)のカウント列
)
#Dateクラスの日付列と、解析使用のタイムポイント列を追加
dat<-dat%>%mutate(date=as.Date(as.POSIXct(seq(as.Date("2011-03-11"), as.Date("2016-03-11"), by="days"))))%>%
#順序を並び替え
dplyr::select(date,-date_jap,Lv_1,Lv_2,Lv_3,Lv_4,Lv_5w,Lv_5s,Lv_6w,Lv_6s,Lv_7,n_over4,n_over1)
#xts型に変換(dygraph用)
dat_xts<-as.xts(read.zoo(dat))2. データの可視化と確認
formattableでデータの確認
続いてデータの確認を行います。表と図で見てみます。まずは表からです。
震度1から7までのそれぞれの発生回数と、震度4以上(n_over4)および全ての震度の発生回数(n_over1)が日ごとに得られています。なお、2011年の3月11日だけ、14時以降のみのデータとなっています。
formattable(dat[1:10,])| date | Lv_1 | Lv_2 | Lv_3 | Lv_4 | Lv_5w | Lv_5s | Lv_6w | Lv_6s | Lv_7 | n_over4 | n_over1 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 2011-03-11 | 127 | 145 | 93 | 41 | 9 | 2 | 0 | 1 | 1 | 54 | 419 |
| 2011-03-12 | 259 | 158 | 51 | 16 | 1 | 0 | 0 | 0 | 0 | 17 | 485 |
| 2011-03-13 | 154 | 70 | 23 | 2 | 1 | 0 | 0 | 0 | 0 | 3 | 250 |
| 2011-03-14 | 107 | 67 | 17 | 3 | 1 | 0 | 0 | 0 | 0 | 4 | 195 |
| 2011-03-15 | 96 | 37 | 19 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 153 |
| 2011-03-16 | 94 | 31 | 17 | 3 | 1 | 0 | 0 | 0 | 0 | 4 | 146 |
| 2011-03-17 | 84 | 38 | 12 | 3 | 0 | 0 | 0 | 0 | 0 | 3 | 137 |
| 2011-03-18 | 70 | 30 | 9 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 110 |
| 2011-03-19 | 96 | 29 | 3 | 3 | 0 | 1 | 0 | 0 | 0 | 4 | 132 |
| 2011-03-20 | 80 | 36 | 7 | 2 | 0 | 0 | 0 | 0 | 0 | 2 | 125 |
dygraphsでの可視化(震度別)
時系列データを図示するときに有効なパッケージ、dygraphsで確認します。
#震度別推移
sumCols<-c(ncol(dat_xts),ncol(dat_xts)-1)
dygraph(dat_xts[,-sumCols],
main = "東日本大震災発生後5年間の地震回数の推移(震度別)")%>%
dySeries("Lv_1",label = "震度1")%>%
dySeries("Lv_2",label = "震度2")%>%
dySeries("Lv_3",label = "震度3")%>%
dySeries("Lv_4",label = "震度4")%>%
dySeries("Lv_5w",label = "震度5弱")%>%
dySeries("Lv_5s",label = "震度5強")%>%
dySeries("Lv_6w",label = "震度6弱")%>%
dySeries("Lv_6s",label = "震度6強")%>%
dySeries("Lv_7",label = "震度7")%>%
dyEvent(x = "2012-03-11", "一年後 ", labelLoc = "top")%>%
dyEvent(x = "2013-03-11", "二年後 ", labelLoc = "top")%>%
dyEvent(x = "2014-03-11", "三年後 ", labelLoc = "top")%>%
dyEvent(x = "2015-03-11", "四年後 ", labelLoc = "top")%>%
dyEvent(x = "2016-03-11", "五年後 ", labelLoc = "top")%>%
dyOptions(fillGraph = TRUE, fillAlpha = 0.3) %>%
dyRangeSelector(height = 40)%>%
dyLegend(width = 600)
dygraphsでの可視化(震度4以上と全震度)
今回の分析で対象とするのは、全震度の地震発生回数と震度4以上(本稿で“強い地震”として扱う)の日別発生回数です。
可視化してみます。日毎の地震発生回数は、5年経過して当時より少なくなっていることが確認されます。
#全震度と震度4以上
dygraph(dat_xts[,sumCols],main = "東日本大震災発生後5年間の地震回数の推移(全震度と震度4以上)")%>%
dySeries("n_over1",label = "全震度合計")%>%
dySeries("n_over4",label = "震度4以上合計")%>%
dyEvent(x = "2012-03-11", "一年後 ", labelLoc = "top")%>%
dyEvent(x = "2013-03-11", "二年後 ", labelLoc = "top")%>%
dyEvent(x = "2014-03-11", "三年後 ", labelLoc = "top")%>%
dyEvent(x = "2015-03-11", "四年後 ", labelLoc = "top")%>%
dyEvent(x = "2016-03-11", "五年後 ", labelLoc = "top")%>%
dyOptions(stepPlot = TRUE,drawPoints = FALSE, pointSize = 2,fillGraph = TRUE, fillAlpha = 0.3,stackedGraph = TRUE) %>%
dyRangeSelector(height = 40)%>%
dyLegend(width = 500)%>%
dyHighlight(highlightSeriesOpts = list(strokeWidth = 2))3. メカニズムの想像
本稿では次の2つの観点から現象を考えてみたいと思います。モデルなので、こういう見方が正解というのはありません。しかし、本稿の考え方よりも妥当なモデルというのは存在するはずです(データ取得初期の地震発生回数はとても多い。ポアソン分布を仮定するのは過分散だと思います、一回差分だけというのもシンプルすぎます。周期性などがあるはずです、どんな周期かはわかりませんが)。色々とデータ生成メカニズムを想像し、検証していただけたら幸いです。
【１】全震度の発生回数をポアソン分布の状態空間モデルで表現(上の図の緑線)
【２】震度４以上の発生回数(上図の青線)は、全震度の発生回数(所与する,上図緑線)と「強い地震の生起確率」をパラメータとする二項分布の状態空間モデルで表現
本稿でメインに推定したいのは、「強い地震の生起確率の推移」となります。この部分だけは、ローデータからそのトレンドが見えていません。
4. Stanモデル(ポアソン分布と二項分布の状態空間モデル)
全震度の発生回数推移を表現するポアソン分布の状態空間モデルと、「強い地震の発生確率」のトレンドを捉える二項分布の状態空間モデル。
これらを同時に扱ってコードにしています。理由は、データとしては手元に得られていない未来の予測(forecast)を31日分行うためです。
generated quantitiesブロックでforecastを行っています。オフセット項わざを使っています。
//PoissonBinomialDynamic.stan
data {
int t;
int n_all[t];
int n_over4[t];
}
parameters {
vector<upper=0>[t+31] mu_p;
vector<lower=0> [t+31] lambda;
real <lower=0> s_p;
real <lower=0> s_lambda;
}
transformed parameters {
vector<lower=0, upper=1>[t+31] p;
p = inv_logit(mu_p);
}
model {
//variance parameters
s_p ~ normal(0,1);
s_lambda ~ normal(0,1);
mu_p[1] ~ normal(-2,s_p); //p=0.1192029
lambda[1] ~ normal(log(n_all[1]),s_lambda);
//observation
for(i in 1:t){
n_over4[i] ~ binomial(n_all[i],p[i]);
n_all[i] ~ poisson_log(lambda[i]);
}
//state(including forecast term)
for(j in 2 : (t+31)){
mu_p[j] ~ normal(mu_p[j-1],s_p);
lambda[j] ~ normal(lambda[j-1],s_lambda);
}
}
generated quantities {
vector<lower=0>[t+31] pred_over4;
vector<lower=0>[t+31] pred_all;
for (k in 1 : (t+31)){
pred_over4[k] = poisson_rng(exp(lambda[k])*p[k]);
pred_all[k] = poisson_log_rng(lambda[k]);
}
}
【dataブロック】
tが時点数、n_all[t]が全震度合計回数、n_over4[t]が震度４以上の地震発生回数です。
【parametersブロック】
mu_pはt+31の長さをもつベクトルで定義した、強い地震の発生確率の、、、なんて言ったら良いんでしょう、“もと”です。正規分布を仮定して、後に逆ロジット変換で0-1の区間に畳みます。+31というのは、31日分の未来の発生確率を考えています。本稿では、未来31日分を欠損値として扱って、推定対象とします。
lambdaは、全震度地震の発生回数のポアソン平均です。
s_pは二項分布モデルのシステムノイズ、s_lambdaはポアソン分布モデルのシステムノイズです。
【transformed parametersブロック】
p = inv_logit(mu_p);
で、強い地震の発生確率を定義しています。inv_logitは逆ロジット変換です。引数に自動でマイナス1が乗算されるのでお気をつけて。
【modelブロック】
観測方程式と状態方程式を分けてfor文で記述しています。一回差分のモデルです。
推定の実行
では、推定の実行です。
事前に収束は確認しております。
datastan = list(n_all = dat$n_over1, n_over4 = dat$n_over4, t = nrow(dat))
#from .stan file
fit<-stan(file="PoissonBinomialDynamic.stan",
data=datastan,
iter=12000,warmup = 2000,chain=4,thin = 25,seed = 1234)MCMC(*´Д`)ﾊｧﾊｧ
5. 解析結果
まずは、サンプルの取り出しを行います。
コードが冗長で大変申し訳ありません。ggmcmcのggs関数を使えば、キレイに取り出せると知りました。 株式会社ホクソエムのブログにある、Kentaro Matsuuraさんのご発表内で解説されています。
p_long<-rstan::extract(fit)$p%>%
as.data.frame()%>%
gather()
lambda_long<-rstan::extract(fit)$lambda%>%
as.data.frame()%>%
gather()
pred_over4_long<-rstan::extract(fit)$pred_over4%>%
as.data.frame()%>%
gather()
pred_all_long<-rstan::extract(fit)$pred_all%>%
as.data.frame()%>%
gather()
df_res<-cbind(p_long[,2],lambda_long[,2],pred_all_long[,2],pred_over4_long[,2])%>%
as.data.frame()
names(df_res)<-c("p","lambda","pred_all","pred_over4")
res_summary<-df_res%>%
mutate(TimePoint=rep(1:1859,each=1600))%>%
group_by(TimePoint)%>%
summarize(p_EAP=mean(p),
p_up=quantile(p,0.975),
p_lo=quantile(p,0.025),
p_up50=quantile(p,0.75),
p_lo50=quantile(p,0.25),
lambda_EAP=exp(mean(lambda)),
lambda_up=exp(quantile(lambda,0.975)),
lambda_lo=exp(quantile(lambda,0.025)),
lambda_up50=exp(quantile(lambda,0.75)),
lambda_lo50=exp(quantile(lambda,0.25)),
pred_over4_EAP=mean(pred_over4),
pred_over4_up=quantile(pred_over4,0.975),
pred_over4_lo=quantile(pred_over4,0.025),
pred_over4_up50=quantile(pred_over4,0.75),
pred_over4_lo50=quantile(pred_over4,0.25),
pred_all_EAP=mean(pred_all),
pred_all_up=quantile(pred_all,0.975),
pred_all_lo=quantile(pred_all,0.025),
pred_all_up50=quantile(pred_all,0.75),
pred_all_lo50=quantile(pred_all,0.25)
)%>%
ungroup()%>%
dplyr::mutate(date=as.POSIXct(seq(as.Date("2011-03-11"), as.Date("2016-04-11"), by="days")),
n_over4 = c(dat$n_over4,rep(NA,nrow(.)-nrow(dat))),
n_over1 = c(dat$n_over1,rep(NA,nrow(.)-nrow(dat)))
)%>%
dplyr::select(date,everything())%>%
as.data.frame()結果の可視化1
まずは、ポアソン分布の状態空間モデルのトレンドを可視化します。ポアソン平均lambdaの95%信用区間を帯で表現します。
#描画したいデータの取り出し
res_all<-c("lambda_EAP","lambda_lo","lambda_up","n_over1")
#全震度と震度4以上
dygraph(res_xts[,res_all],main = "東日本大震災発生後5年間で発生した地震が震度4以上である確率(日別)")%>%
dySeries("n_over1",label = "全震度地震発生回数(ローデータ)")%>%
dySeries(c("lambda_lo","lambda_EAP","lambda_up"),label = "全震度地震発生回数の推定平均")%>%
dyEvent(x = "2012-03-11", "一年後 ", labelLoc = "top")%>%
dyEvent(x = "2013-03-11", "二年後 ", labelLoc = "top")%>%
dyEvent(x = "2014-03-11", "三年後 ", labelLoc = "top")%>%
dyEvent(x = "2015-03-11", "四年後 ", labelLoc = "top")%>%
dyEvent(x = "2016-03-11", "五年後 ", labelLoc = "top")%>%
dyOptions(stepPlot = FALSE,drawPoints = FALSE, pointSize = 2,fillGraph = FALSE, fillAlpha = 0.3,stackedGraph = FALSE,
colors = RColorBrewer::brewer.pal(5, "Set1")) %>%
dyRangeSelector(height = 40)%>%
dyLegend(width = 600)%>%
dyHighlight(highlightSeriesOpts = list(strokeWidth = 2))結果の可視化2
続いて、全震度の地震発生回数をパラメータとして処よした時の強い地震の発生確率トレンドを可視化します。
#xts型に変換(dygraphやhighcharter用)
res_xts<-as.xts(read.zoo(res_summary))
#描画したいデータの取り出し
res_prob<-c("p_up","p_lo","p_EAP")
#全震度と震度4以上
dygraph(res_xts[,res_prob],main = "東日本大震災発生後5年間で発生した地震が震度4以上である確率(日別)")%>%
#dySeries(c("p_lo50","p_EAP","p_up50"))%>%
dySeries(c("p_lo","p_EAP","p_up"),label = "震度4以上の確率EAP(95%CI)")%>%
dyEvent(x = "2012-03-11", "一年後 ", labelLoc = "top")%>%
dyEvent(x = "2013-03-11", "二年後 ", labelLoc = "top")%>%
dyEvent(x = "2014-03-11", "三年後 ", labelLoc = "top")%>%
dyEvent(x = "2015-03-11", "四年後 ", labelLoc = "top")%>%
dyEvent(x = "2016-03-11", "五年後 ", labelLoc = "top")%>%
dyOptions(stepPlot = FALSE,drawPoints = FALSE, pointSize = 2,fillGraph = FALSE, fillAlpha = 0.3,stackedGraph = FALSE,
colors = RColorBrewer::brewer.pal(3, "Set2")) %>%
dyRangeSelector(height = 40)%>%
dyLegend(width = 600)おわりに
地震そのものの回数は減っています。
しかし、地震発生回数に占める強い地震(震度4以上と本稿では定義)の発生確率は、震災後あまり変わっていないことがわかりました。
強い地震の発生には、不等間隔な周期があるようにみえなくもありません。
もう一歩踏み込んだ解析が必要だと考えています。