Miyagi Earthquake & Inverse Gaussian Distribution

宮城県沖地震は定期的になんども大地震が起きている。

  • 1793/02/17
  • 1835/07/20
  • 1861/10/21
  • 1897/02/20
  • 1936/11/03
  • 1978/06/12
  • 2011/03/11

服部哲弥さんの「統計と確率の基礎」(2006年。学術図書出版社)という本を読んでいたところ、第4章「点推定の考え方」で点推定を用いて

  • 宮城県沖地震発生間隔
  • いまから何年後に大地震が起きるか?

を推定していた。

出版当時(2006年)は2011年の大地震が起こる前だったので、2011年3月11日の大地震を反映して、次の大地震を Rで推定してみる

MODEL
同書籍では次のようにモデルしている。
逆ガウス分布(Inverse Gaussian Distribution) \rho(t) を用いて、宮城県沖地震が発生してから次に発生するまでの時間間隔を表す確率変数を X と置くときそれが a 以上 b 以下となる確率は次の式で計算できる。
P[a \leq X \leq b ] = \int_a^b \rho(t) dt

ここで 逆ガウス分布(Inverse Gaussian Distribution) は次の密度関数を持つ連続分布
\rho(t) = \frac{1}{\sqrt{2 \pi v t^3 / m ^3}} e ^{-\frac{ (t-m) ^ 2}{2 v t / m}}, t > 0, m := 母集団の平均, v := 母集団の分散

R:点推定と分布へのあてはめ

赤が経験分布(観測された周期)。黒塗りが経験分布をもとにパラメータを点推定して逆ガウス分布に当てはめた結果。

library(ggplot2)
library(lubridate)

incident <- c("1793/02/17", "1835/07/20", "1861/10/21", "1897/02/20", "1936/11/03", "1978/06/12", "2011/03/11")
incident <- ymd(incident) # convert to date object
incident_interval <- (incident[-length(incident)] - incident[-1]) / dyears(1)
# 42.44384 26.27397 35.35890 39.72603 41.63288 32.76712
density <- function(t, m, v) {
   (1/ sqrt(2 * pi * v * t^3 / m^3)) * exp(- (t-m) ^ 2 / (2 * v * t / m))
}

m <- mean(incident_interval) # 36.36712
v <- var(incident_interval)  # 38.35701
year <- 1:80
dat <- data.frame(year=year, prob=density(year, m=m, v=v))
ggplot(dat, aes(x=year, y=prob)) + geom_area()+ geom_vline(xintercept=incident_interval, color='red') + opts(title='宮城県沖地震周期 + Inverse Gaussian Distribution')

lubridate
日付計算はlubridate を用いている。R 標準関数で日付計算をするのは、かなり面倒。次のチュートリアルは非常に有用

Garrett Grolemund, Hadley Wickham, “Dates and Times Made Easy with lubridate”, Vol. 40, Issue 3, Apr 2011, Journal of Statistical Software.

R:次の大地震は何年後?

横軸に時間をとって、次回の大地震発生までの発生確率を図にしたのが次。

dat$cumulative <- cumsum(dat$prob)
ggplot(dat, aes(x=year, y=cumulative)) + geom_area() + ylim(0, 1) + opts(title='When will the next Miyagi earthquake happen?') + xlab('year from 2011/03/11')

integrate
積分するなら
integrate(density, lower=foo, upper=bar)$value
のようにすればできる。

MEMO

Advertisements
Tagged with: , , , , , ,
Posted in algorithm, R

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Archives
%d bloggers like this: