ほぼ毎年のように聞く「100年に一度の大雨」とは何なのか？を極値統計学で考える


# はじめに
極値統計学については、Qiitaには[極値統計学に軽く触れてみる](https://qiita.com/cotton-gluon/items/bd1a12e22cf4c2d8b805)や[正規分布の最大値、最小値の分布](https://qiita.com/raidofuwa/items/5bbb2c58a4ac457992bc)といった素晴らしい記事がすでにアップされており、
[Zennにも記事](https://zenn.dev/imaimai/scraps/259efba2e8ccf4e8787c)がアップされているので、こちらも参照いただいたほうが良いです。

# 「100年に一度」と極値統計学
よくマスコミや政府が使う「**100年に一度の大雨**」という表現、最近この表現が使われすぎていてやや陳腐化しているきらいがあります。
あまりにも頻繁に聞くものですから、「100年に一度なんて言って大したことないじゃないか、不安を煽りたいだけないんじゃないか？」と思ってしまいがち。
でも良く考えてみると、「この雨は100年に一度の大雨なんです！」というのは何を根拠に言っているのでしょう？
**実は、このことは極値統計学によって出しています！**
 他には、[腐食材と石油タンク底板の残存寿命](https://ms-laboratory.jp/pdf/kyoku/kt2.htm)の推定や、[市場収益のモデル化](http://infoshako.sk.tsukuba.ac.jp/~databank/pdf/902.pdf)、[陸上競技においてどこまで人類は記録を伸ばせるか?](https://www.tandfonline.com/doi/abs/10.1198/016214508000000698)という抽象的な問に対しても極値統計学は一定の解をもたらしています
(これによると100 m走は9秒29まで短縮可能らしい)。

なぜこのようなことができるのかというと、**極値統計学とは最大値や最小値が得られる確率を求めるための手法だからです。**
この記事では、極値統計学の簡単な理論と、アメダスのデータを用いた「100年に一度の降水量」の推定を行おうと思います。

# 極値分布

「最大値や最小値が得られる確率」と言われても具体的にはどういうものなのかわからないので、少し実験してみましょう。

手元にサイコロがあるとします。これらの出る目は同様に確からしく(つまり離散一様分布)、サイコロを3回降ったときの最大値を記録するという操作を2000回行ったとします。
ここで、**どの目が最も最大値として記録されていそうだろうか?** という問いを立てます。
平均であれば3.5になりますが、最大値が6になる確率は非常に高そうです。逆に、おそらく最大値が1となる確率は非常に低いでしょう。

実際にjuliaで書いてみると、

```setup.jl
using Plots
using Distributions
using StatsBase
using StatsPlots
```

```maximum_Uniform.jl

sample_max = zeros(Int, 6)
anim = Animation()

for i in 1:2000
    max = maximum(sample(1:6, 3, replace=true))
    sample_max[max] += 1
    plt = bar(sample_max ./ i, xlims=(0.5, 6.5), ylims=(0, 0.5), label="", title = "Uniform dist: $i")
    frame(anim, plt)

end

@time gif(anim, "maximum_Uniform.gif", fps = 20)

```
![maximum_Uniform.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/269013/7747d4c2-90e8-aece-5f91-d40f2e63c400.gif)

最大値が6になる確率は実際に高いようです、およそ0.4ぐらい。
しかも、なにかの分布に収束していそうな雰囲気があります。
これは偶然でしょうか？他の分布から得られた場合ではどうだろう、ということで正規分布から3回データをとり、その最大値を取るとしましょう。

```maximum_Normal.jl
anim = Animation()
sample_max = []

for i in 1:2000
    push!(sample_max, maximum(rand(Normal(0, 1), 3)))
    plt = histogram(sample_max, normalize=:probability,bin = -2:0.25:5,title = "Normal dist: $i")
    frame(anim, plt)
end

@time gif(anim, "maximum_Normal.gif", fps = 20)

```
![maximum_Normal.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/269013/67ce3c4e-eb1e-68d3-528a-0903f29771e3.gif)

正規分布のような雰囲気がするが、やはり収束していそう。
指数分布やコーシー分布でも同様に描いてみます。
julia言語なら、randの中身を少し変えるだけで良いので便利。


```maximum_Exponential.jl
anim = Animation()
sample_max = []

for i in 1:2000
    push!(sample_max, maximum(rand(Exponential(), 3)))
    plt = histogram(sample_max, normalize=:probability,bin = 0:10, title = "Exponential dist: $i")
    frame(anim, plt)
end

@time gif(anim, "maximum_Exponential.gif", fps = 20)

```

![maximum_Exponential.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/269013/b277c38f-75c6-d8c0-71a9-1d7abd419249.gif)

```maximum_Cauchy.jl
anim = Animation()
sample_max = []

for i in 1:2000
    push!(sample_max, maximum(rand(Cauchy(10), 3)))
    plt = histogram(sample_max, normalize=:probability, title = "Cauchy dist: $i")
    frame(anim, plt)
end

@time gif(anim, "maximum_Cauchy.gif", fps = 20)

```

![maximum_Cauchy.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/269013/e3b8ba20-5784-98c8-61a1-8b955c08cd0a.gif)


また正規分布由来の分布とは異なる収束の仕方をしていますね。色々趣があって面白い！

このように、$\\{X_i\\}_{i=1}^{n}, n = 1,2,3, ...$で考えられる確率変数列の極値

(今回の場合は最大値、すなわち $ Z_{n} = max(X_1, X_2, ... , X_n)$)を考えるとその分布はどの様になるのか？ということを考える学問が極値統計学です。
このように分布が収束する場合、この分布を極値分布と呼びます。
端的に言えば「**ある確率分布にしたがう独立な$N$個のサンプルがあったときに、そのサンプルの極値 (最大値や最小値など) の分布**」のことを指しているということです。

# アメダスのデータを使った「100年に一度の大雨」推定

極値分布は、理論的には**Gumbel分布、Fréchet分布、そして負のWeibull分布の3種類のどれか**になることがわかっています。
f(・)を確率密度関数、F(・)を累積分布関数とすると以下のような式で記述されます(この辺りは統計検定とか受けている人向けですかね？)。

+ Gumbel分布: $f(x; \alpha, \beta) = \frac{1}{\beta}\exp(-\exp(\frac{x-\alpha}{\beta})+\frac{x-\alpha}{\beta}), F(x;α,β)=\exp(−\exp(-\frac{x−\alpha}{\beta}))$
(もとの分布が指数分布や正規分布でこの分布に収束)

+ Fréchet 分布: $f(x; \alpha) = \alpha x^{-\alpha-1}\exp(-x^{-\alpha}), F(x; \alpha) = \exp(-x^{-\alpha})$
(もとの分布がコーシー分布だとこの分布に収束)

+ 負のWeibull 分布: $f(x; \alpha, \beta) =\frac{\alpha} {\beta^\alpha}x^{\alpha-1}\exp(-(\frac{x}{\beta})^\alpha), F(x;\alpha, \beta) = 1-\exp(-(\frac{x}{\beta})^\alpha) $
(もとの分布が一様分布だとこの分布に収束)

少し統計に詳しい人ですと、中心極限定理をご存知かと思われます。
雑に言えば、「ある分布から得られた平均値の分布は正規分布に従う」という定理ですが、コーシー分布はその中心極限定理に従わないことで有名です。
そのコーシー分布であっても、「ある分布から得られた最大値の分布はちゃんと収束している、と言うのは個人的に非常に興味深いところです。

更にこれらを統合した関数は**一般化極値分布（Generalized Extreme Value distribution; GEV）** と呼ばれます。
しかし、今回はGEVまでには立ち入らず、上の分布を利用して東京で得られた[アメダスのデータ](https://www.data.jma.go.jp/risk/obsdl/index.php)を用いて「100年に1度の大雨」とは何かを考えてみましょう。
まず、アメダスのデータ(1872/01/01 ~ 2021-12-31)を[気象庁のHP](https://www.data.jma.go.jp/risk/obsdl/index.php)からダウンロードし、東京都での日々の降水量をplotします。(ここからはR言語で描いています、ご了承下さい)
2022-12-28追記: [2022年末、Rでcsvを読む](https://zenn.dev/eitsupi/articles/r-read-csv-2022
)を参考に一部書き直しました。

```data_plot.R
library(tidyverse)
library(magrittr)
library(fs)
library(here)
library(lubridate)

df <- 
read_csv(
  dir_ls(
    path = here(),
    glob = "*.csv"
  ) , 
  locale = locale(encoding = "cp932"),
  skip = 5, 
  col_names = FALSE
)%>% 
  rename(
    "date" = "X1", "降水量_mm" = "X2"
  ) %>% 
  na.omit() %>%   
  select(date, 降水量_mm) %>% 
  mutate(date = as_date(.$date)) %>% 
  arrange(desc(降水量_mm))

df %>% 
ggplot(aes(date, 降水量_mm)) + geom_point()
```

![降水量.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/269013/f825b1e4-2fc0-41c5-bd53-5f589a4db3d2.png)


1日に300 mm以上降ったのは[狩野川台風](https://typhoon.yahoo.co.jp/weather/calendar/108/)という台風の際に記録されたそうです。

さて、今回知りたいのは「100年に1度の大雨」なので、**年間の最大雨量**を集計します。

```max_rain.R
df_max <- 
    df %>% 
    mutate(rain_year = year(date)) %>%
    group_by(rain_year) %>% 
    summarise(max = max(降水量_mm))

> df_max
# A tibble: 122 × 2
   rain_year   max
       <dbl> <dbl>
 1      1901  80.5
 2      1902  66.3
 3      1903 165.
 4      1904 126.
 5      1905 107. 
 6      1906 172.
 7      1907 112.
 8      1908  89.2
 9      1909  74.1
10      1910 147.
# … with 112 more rows
# ℹ Use `print(n = ...)` to see more rows
```

さて上述した通り、最大値を推定するモデルとして、理論的にはGumbel分布、Fréchet分布、負のWeibull分布の3種類のどれかになることがわかっています。
ここで、降雨量は指数分布に従うと仮定すれば、最大値を推定するモデルとしてGumbel分布が有効になるはずです。
(なぜ指数分布にしたかは割愛します、すいません。)

なので、最大値をヒストグラムで描画してみましょう。

```histo.R
df_max %>%
ggplot(aes(max)) +
geom_histogram(
    aes(y = ..density..),
    position = "identity",
    lwd = 0.2
    )+
    xlab("年間最大降水量") +
    ylab("probability density")
```
![降水量_histo.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/269013/75e8ef82-e15b-61e0-7ea0-31f8b9c44998.png)

確かに、Gumbel分布に近い形を取っているように見えます。なので、パラメータを最尤法で推定しましょう。
Rのパッケージに{evd}と呼ばれるものがあり、Gumbel分布のパラメータ推定とggplotによるグラフの描写がサポートされています。

```estimate.R
library(evd)
fgumbelx(df_max$max %>% na.omit(), std.err = FALSE, method = "Nelder-Mead")

Call: fgumbelx(x = df_max$max %>% na.omit(), std.err = FALSE, method = "Nelder-Mead") 
Deviance: 1186.366

Estimates
  loc1  scale1    loc2  scale2
 73.65   51.27   73.65   25.24

Optimization Information
  Convergence: successful
  Function Evaluations: 299
```

このloc1, scale1と呼ばれるものが推定されたパラメーターです。早速カーブフィッティングしてみましょう。

```curve_fitting.R
df_max %>%
ggplot(aes(max)) +
geom_histogram(
    aes(y = ..density..),
    position = "identity",
    lwd = 0.2
    )+
    xlab("年間最大降水量") +
    ylab("probability density") +
    stat_function(fun = dgumbel, args = list(loc = 73.90, scale = 51.27))

```

![降水量_histo_fit.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/269013/56c1a38b-7d47-3816-7844-14c0db7cd05b.png)

いい感じですかね？では、100年に一度の大雨がどの程度の雨量になるかを推定しましょう。
100年に1度ということなので、上側1%の降水量を求めれば良いはずです。
```qgumbel()```という関数で導出可能ですので、出してみましょう。

```qgumbel.R
qgumbel(0.01, loc = 73.65, scale = 51.27, lower.tail = FALSE)
[1] 309.4997
```
ということで、東京における「100年に1度の大雨」とは、**およそ309 mmの雨が1日に降ること**と推定できます！
実は気象庁の[異常気象リスクマップ](https://www.data.jma.go.jp/cpdinfo/riskmap/qt_table.html)に各地点での100年確率降水量が出ています。
そこのデータを見ると、**289 mm**となっています。計算に使っている期間が違う(気象庁は1901～2006年で集計)ので少しずれてしまっていますが、おおよそ近い値になっています。

# 最後に

ということで、今回は「100年に1度の大雨」について極値統計学を用いて推定しました。
これ以外にも様々な分野で扱われていますので、是非興味があったら他の記事や教科書で学んでみるのも良いかもしれません。