読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

大森・宇津公式を強度関数とした非定常ポアソン過程の重ね合わせによる熊本地震の余震の分析

スクレイピング

使用するツールは R です。

データは 気象庁|地震情報 から取得しました。

rvest パッケージを使うとかんたんです。

また XML パッケージの readHTMLTable 関数を使う方法もあります。

R から HTML の表を読み込む - 廿TT

速報は一定時間たつと流れていってしまうため、数日おきにデータをダウンロードしてきて dplyr の full_join でマージしています。

データを読み込んだら Nippon パッケージを使って全角英数字を半角になおしたりして整形しました。

最終的には震央地名が熊本県の行のみを抜き出しています。

library(rvest)
html1 <- read_html("http://www.jma.go.jp/jp/quake/quake_local_index.html")

dat_row <- html_table(html1,header = TRUE)[[4]]

library(Nippon)
Mag <-sapply(tab_row[,4],zen2han)
Mag <-as.numeric(substr(Mag,2,4))

times <-sapply(tab_row[,2],zen2han)
day <- paste0("2016/04/",substr(times,1,2))
hour <-paste(substr(times,4,5),substr(times,7,8),sep=":")
times <-as.POSIXct(paste(day,hour))
flag <-grepl("熊本県",tab_row[,3])
kumamoto <-data.frame(time=times[flag],magnitude=Mag[flag]) %>%
  arrange(time)%>%
  distinct(time)%>%
  mutate(time_int = as.integer(time-min(time)))

加工したデータは from http://www.jma.go.jp/jp/quake/quake_local_index.html · GitHub にこっそりおいて置きます。

より本格的に地震の分析をしたい方は 気象庁|震源リスト などからデータを取得したほうが正確だと思います。

記述統計

最初にいくつかグラフを書きました。

横軸を時間、縦軸をマグニチュードとしたマーク付き点過程チックなグラフ、

f:id:abrahamcow:20160429005928p:plain

横軸を時間、縦軸を地震の累積発生回数とした計数過程チックなグラフ、

f:id:abrahamcow:20160429010005p:plain

マグニチュードヒストグラム

f:id:abrahamcow:20160429010911p:plain

地震の発生間隔のヒストグラム

f:id:abrahamcow:20160429011003p:plain

地震の発生日時のヒストグラム

f:id:abrahamcow:20160429011030p:plain

です。

library(ggplot2)
theme_set(theme_bw(20,"HiraMaruProN-W4"))
ggplot(kumamoto)+
  geom_segment(aes(x=time,y=magnitude,xend=time,yend=0))+
  scale_x_datetime(date_labels= "%m/%d")

ggplot(kumamoto)+
  geom_point(aes(x=time,y=1:length(time),size=magnitude),alpha=0.6)+
  scale_x_datetime(date_labels= "%m/%d")+
  ylab("cumulative counts")+
  geom_rug(aes(x=time))

ggplot(kumamoto) +
  geom_histogram(aes(x=magnitude),fill="white",color="black")
ggplot() +
  geom_histogram(aes(x=diff(kumamoto$time_int)/3600),fill="white",color="black")+
  xlab("発生間隔(時間)")

ggplot(kumamoto) +
  geom_histogram(aes(x=time),fill="white",color="black")+
  scale_x_datetime(date_labels = "%d日%H時")+
  xlab("発生日時")

累積発生回数のグラフを見ると4月16日ごろに地震の多発がいったん収まってきてから、また盛り上がっていることがわかります。

このダイナミクスを再現できるようなモデルを考えます。

大森・宇津公式を強度関数とした非定常ポアソン過程

余震に関する経験法則として、大森・宇津公式(修正大森公式)というものが知られています。

これは本震からの経過時間 t における単位時間あたりの余震回数が次のようになるというものです。

\displaystyle \frac{K}{(t+c)^p}

各パラメータは

  • K:余震の多さ。
  • c:本震直後の余震の少なさ。0.1ぐらいの値をとる場合が多い。
  • p:時間経過に伴う減衰度。1ぐらいの値をとる場合が多い。

と説明されます。

余震 - Wikipedia

余震活動を非定常ポアソン過程と考え、その強度関数(intensity function)に大森・宇津公式を採用して、パラメータを推定します。

時刻 t における強度が  \lambda(t) の非定常ポアソン過程を、区間 (0, T] で観測した場合、尤度関数は以下のようになります。

\displaystyle L(\beta,m)= \prod_{i=1}^{n} \lambda ( t_{i}) \exp\left(-\Lambda(T)\right)

今回は、

\displaystyle \lambda(t)= \frac{K}{(t+c)^p}

です。そのため

\displaystyle \Lambda(t) = \int^{t}_{0} \lambda(u) \, du= K\{(c+x)^{1-p}\}/(1-p)

となります。

各パラメータの推定値は K = 0.14、c = -2.28、p = 0.50 となりました。

f:id:abrahamcow:20160429013019p:plain

黒い線が実際の地震の累積発生回数、赤い線が大森・宇津公式による推定値です。

なにを持って前震、本震、余震とするかむずかしいところではありますが、4月14日を本震と見ると、あまり当てはまりはよくありません。

dat <- kumamoto$time_int[-1]

LLOmori <- function(par){ #尤度関数
  K <- par[1]
  c <- par[2]
  p <- par[3]
  sum(log(K)-p*log(dat+c))-K*((c+max(dat))^(1-p))/(1-p)
}
fit1 <-optim(c(1,0,0.1),LLOmori,control = list(fnscale=-1,maxit=10000))
Lambda <- function(x,K,c,p){
  K*((c+x)^(1-p))/(1-p)
}
lambda <- function(x,K,c,p){
  K/((x-c)^p)
}

ggplot(kumamoto)+
  geom_step(aes(x=time,y=1:length(time)))+
  geom_line(aes(x=time,y=Lambda(time_int,fit1$par[1],fit1$par[2],fit1$par[3])),colour="red")+
  scale_x_datetime(date_labels= "%m/%d")+
  ylab("cumulative counts")+
  geom_rug(aes(x=time))

round(fit1$par,2)

非定常ポアソン過程の重ね合わせ

4月16日ごろに地震の多発がいったん収まってきてから、また盛り上がってくるダイナミクスをモデル化する方法の一つに、ポアソン過程をふたつ重ね合わせることが考えられます。

強度関数  \lambda _1(t)ポアソン過程と強度関数  \lambda _2(t)ポアソン過程の和は、強度関数  \lambda (t)=\lambda _1(t)+\lambda _2(t)ポアソン過程になります。

(こういうのをデュレットはスーパーポジションと呼んでいます。)

確率過程の基礎

確率過程の基礎

累積の強度関数は  \Lambda (t) = \int^{t}_{0} \lambda _1(u) \,du+\int^{t}_{0}\lambda _2(u) \,du となります。

変化点を τ として、\lambda _1(t)\lambda _2(t) に両方とも大森・宇津公式を採用します。

\displaystyle \lambda_1(t)= \frac{K_1}{(t+c_1)^{p_1}}

\displaystyle \lambda_2(t) =  \frac{K_2}{(t-\tau+c_2)^{p_2}}

対数尤度関数を R で表現すると、以下のようになります。

LLWOmori <- function(par){
  K1 <- par[1]
  c1 <- par[2]
  p1 <- par[3]
  K2 <- par[4]
  c2 <- par[5]
  p2 <- par[6]
  tmax <- max(dat)
  sum(ifelse(dat>tau,
             log(lambda(dat,K1,c1,p1)+lambda(dat-tau,K2,c2,p2)),
             log(lambda(dat,K1,c1,p1)))) -
    (Lambda(tmax,K1,c1,p1)+Lambda(tmax-tau,K2,c2,p2))
}

ただしこの対数尤度を最大化するようにパラメータを求めるのは難しいため、最小二乗法でパラメータを推定しました。

パラメータの推定は以下のよう手順で行いました。

  1. 最初の地震発生時点を τ として固定しそれ以外のパラメータを推定する
  2. 次の地震発生時点を τ として固定しそれ以外のパラメータを推定する
  3. 2 を繰り返し、中でも残差平方和が最小になったパラメータの組を最終的な推定値とする

いくつか収束してないやつもありますが、めんどうなのでしかたないものとしてそのままにしました。

f:id:abrahamcow:20160429031915p:plain

当てはまりはぼちぼちといったところですが、さきほどよりはましな推定値が得られた言えそうです。

AIC で比較すると、シンプルな大森・宇津公式のほうは 4835.01 大森・宇津公式をふたつ組み合わせたほうは 4784.858 でした。

AIC の小さいほうのモデルを採択すると、本震が二回起こっているとみなすのが自然ということになります。

ただ、ぼくの経験では AIC でモデルを選択すると、変化点ありのモデルが採択されやすい印象を持っており、AIC をこういうふうに使っていいのかは意見が欲しいところです。

パラメータの推定値は
 K_1=0.49, c_1= -2.91, p_1= 0.70
 K_2= 0.82, c_2= -0.94,p_2  =0.69
でした。

K で比較すると2番目の地震のほうが規模が大きいようです。

n <- length(dat)
Omorilist <- vector("list",n)
for(i in 1:n){
  resid <- resid0(dat[i])
  fit2 <- optim(c(fit1$par,fit1$par),resid)
  Omorilist[[i]] <- optim(fit2$par,resid,control = list(maxit=20000),method = "CG")
}
for(i in which(sapply(Omorilist, function(x)x$convergence)==1)){
  resid <- resid0(dat[i])
  fit2 <- Omorilist[[i]]
  Omorilist[[i]] <- optim(fit2$par,resid,control = list(maxit=20000),method = "CG")
}

table(sapply(Omorilist, function(x)x$convergence)) #収束判定
s <- which.min(sapply(Omorilist, function(x)x$value))
tau <- dat[s]
fit2 <- Omorilist[[s]]

ggplot(kumamoto)+
  geom_step(aes(x=time,y=1:length(time)))+
  geom_line(aes(x=time,y=Lambda2(time_int,fit2$par[1],fit2$par[2],fit2$par[3],fit2$par[4],fit2$par[5],fit2$par[6],dat[s])),colour="red")+
  scale_x_datetime(date_labels= "%m/%d")+
  ylab("cumulative counts")+
  geom_rug(aes(x=time))

#AIC
-2*fit1$value+2*3
-2*LLWOmori(fit2$par)+2*7

むすびに

以上が4月26日時点までのデータを使った結果ですが、4月27日以降さらに地震が頻発しており、変化点がもう一つあるとみなしたほうがいいような状況になっています。

現実の自然現象はむずかしい、一筋縄ではいかないものでした。

f:id:abrahamcow:20160429004845p:plain

f:id:abrahamcow:20160429033328p:plain

いいかげん余震おさまってほしい。