統計: ポアソン分布を使って今後の大地震が起こる確率を求めてみる

ポアソン分布というのは、ごくまれに起こるような事象の確率分布をいう。この説明を聞いても何のこっちゃという感じだけど、これを使うと滅多に起こらないようなことがある時間内にどれくらい起こりそうなのかが分かる。もちろん、別に未来を予知しているわけではなく、あくまで過去の生起確率からはそのように求まるというのに過ぎない。

定義

ポアソン分布は次の数式で求められる。

$P(X = k) = \frac{\lambda^{k} e^{- \lambda}}{k!}$

まず $\lambda$ が、その滅多に起こらないような事象の単位時間あたりの生起確率になっている。そして $k$ が単位時間あたりにその事象が何回起こるかを表す変数になっている。

日本で今後一年間に震度 7 の地震が起こる確率を求めてみよう

それでは、手始めにポアソン分布を使って今後一年間に震度 7 の地震が日本の何処かで起こる確率を求めてみよう。ちなみに、実際に政府が今後 X 年間にとある地域でマグニチュード Y 以上の地震が起こる確率を計算するときはポアソン分布を使っているらしい。

まず、ポアソン分布を使うにはその事象が単位時間あたり何回起こっているかを調べる必要がある。今回は Wikipedia にある震度 7 の項目をベースに計算してみよう。ここには、日本で過去に発生した震度 7 の地震の記録がある。

震度7 - Wikipedia

上記を見ると震度 7 の地震は 1995 年から 2016 年の 21 年間で 5 回発生していることが分かった。生起確率は 1 年あたり $\frac{5}{21}$ ということになる。細かくやるなら震度 7 の定義ができたタイミングから今日までで計算すべきなんだろうけど、ひとまずざっくりということで。

まずは、今後一年間に震度 7 の地震が起こらない確率を求める

単位時間あたりの生起確率が分かったところで、まずは今後一年間に震度 7 の地震が起こらない確率から求める。この計算をする理由は後述する。

事象が一度も起こらない場合なので、回数には $k = 0$ を代入する。そして生起確率は $\lambda = \frac{5}{21}$ になる。実際にポアソン分布に代入して式を整理してみよう。

$P(0) = \frac{\lambda^{k} e^{-\lambda}}{k!} = \frac{\frac{5}{21}^{0} e^{-\frac{5}{21}}}{0!} = e^{-\frac{5}{21}}$

すると $P(0) = e^{-\frac{5}{21}}$ になることが分かった。

実際に値を計算してみる

実際に、上記の値を Python で計算してみよう。

まずは Python インタプリタを起動する。

$ python --version
Python 3.5.2
$ python

そして数学系の math パッケージをインポートしたら、先ほどの数式になるように式を入力する。

>>> import math
>>> math.e ** -(5/21)
0.7881276277453111

すると 0.78812... という値が得られた。これはつまり、今後一年間で震度 7 の地震が起きない確率は約 78% ということだ。

何で起きない確率を計算したのか？

ポアソン分布は単位時間あたりに事象が $k$ 回起きる確率を求める分布だった。つまり、単位時間あたりに少なくとも 1 回以上起きる確率を愚直に計算しようとすると、こうなってしまう。

$P(X \gt 0) = \displaystyle \sum_{k = 1}^{\infty} \frac{\lambda^{k} e^{-\lambda}}{k!}$

上記は単位時間あたりに事象が $k$ 回起こる確率を、1 から無限まで足し合わせていく計算になっている。

試しに Python を使って、実際に $k$ 回起こる確率を 0 から 4 までの間で見てみよう。

>>> f = lambda n: ((5 / 21) ** n) * (math.e ** -(5 / 21)) / math.factorial(n)
>>> f(0)
0.7881276277453111
>>> f(1)
0.187649435177455
>>> f(2)
0.022339218473506547
>>> f(3)
0.0017729538471036941
>>> f(4)
0.00010553296708950561

もちろん、実際にはこんな計算をしていく必要はない。全ての事象が起こる確率は足せば 1 になるはず。

$\displaystyle \sum_{k = 0}^{\infty} \frac{\lambda^{k} e^{-\lambda}}{k!} = 1$

だとすると、事象が少なくとも一回以上起こる確率は 1 から起こらない確率を引けば良いことが分かる。

$P(X \gt 0) = 1 - P(0)$

このために、まずは起こらない確率を計算したというわけ。

今後一年間に震度 7 の地震が起こる確率を求める

ということで、今度は起こる確率を求めてみる。

先ほど記述した通り、少なくとも一回以上起こる確率は 1 から起こらない確率を引けば良い。

$P(X \gt 0) = 1 - P(0) = 1 - \frac{\lambda^{0} e^{-\lambda}}{0!} = 1 - e^{-\lambda} = 1 - e^{- \frac{5}{21}}$

上記を元に、今度は起こる確率を計算する。

>>> 1 - math.e ** -(5/21)
0.2118723722546889

これで、今後一年間に震度 7 の地震が日本で起こる確率は約 21% と分かった。

うんうん…で？

まあ、といっても上記に関してはポアソン分布を使うまでもないでしょうという感じ。だって、過去の単位時間あたりに起こる生起確率は分かっているんだから。

$\frac{5}{21} = 0.23809523809523808$

単純に一年あたり何回起こっているかを計算した結果と大して変わりがない。

日本で今後 10 年間に震度 7 の地震が起こる確率を求めてみよう

ポアソン分布を使うなら単位時間をもっと伸ばしたり縮めたりしないと面白みがない。次は試しに 10 年間で起こる確率を求めてみよう。

既に分かっている通り、起こる確率は 1 から起こらない確率を引いて計算する。ただし、今度は期間が 10 倍になっているので生起確率も 10 倍にする。

$P(X \gt 0) = 1 - P(0) = 1 - e^{-\lambda} = 1 - e^{-\frac{10 \times 5}{21}}$

それでは、実際に Python で計算してみよう。

>>> 1 - math.e ** -(10 * 5/21)
0.90753752393708

ということで、今後 10 年間に震度 7 の地震が日本で起こる確率は約 90% と分かった。

今後 30 年間なら？

これを 30 年間まで伸ばすと、なんと確率は約 99.92% まで上がる。

>>> 1 - math.e ** -(30 * 5/21)
0.99920950967688

地震が起こる世界線と起こらない世界線があるとしたら、なんと起こらない世界線に到達できる可能性は $\frac{1}{1265}$ だ。

>>> 1 / math.e ** -(30 * 5/21)
1265.0376238043307

あくまで日本の何処か、とはいえ残りの人生の中で震度 7 の地震はほぼ間違いなく起こることが分かる。日々の備えが必要だな…。

明日、震度 7 の地震が起こる確率は？

そして、ポアソン分布の面白いところは期間を縮めてもちゃんと計算できるところ。明日の確率を求めたいなら、単位時間を一日にする。

$P(X \gt 0) = 1 - e^{-\frac{5}{21 \times 365}}$

Python で計算してみよう。

>>> 1 - math.e ** -(5/(21*365))
0.0006521030091632962

上記から、明日震度 7 の地震が日本の何処かで起こる確率は約 0.06% と分かった。

今後一週間なら？

生起確率の単位時間を一日にして 7 倍する。

$P(X \gt 0) = 1 - e^{-\frac{7 \times 5}{21 \times 365}}$

>>> 1 - math.e ** -((7*5)/(21*365))
0.004555800758262785

一週間だと確率は約 0.4% になることが分かった。

まとめ

ポアソン分布を使うと滅多に起こらない事象が今後どれくらいの確率で起こりそうなのかを計算できる。

スマートPythonプログラミング: Pythonのより良い書き方を学ぶ

作者: もみじあめ
発売日: 2016/03/12
メディア: Kindle版
この商品を含むブログを見る

CUBE SUGAR CONTAINER

技術系のこと書きます。