100年に一度の洪水を予測する: 極値統計学入門

はじめに

極値統計学とは、あるデータが取る最大値や最小値の分布を考える理論です。ところで対照的に、統計学の多くの理論では、データの平均やばらつきに関するものが多いです。特に、大抵のデータの平均値や合計値といった統計量は、中心極限定理によって正規分布に従ってくれるため非常に便利ですよね。このように、あるデータの統計量が別の分布に従うことを分布収束と言いますが、これと同じことが極値統計学でも起こります。すなわち、ある確率分布に従うデータの最大値や最小値が、別の確率分布に従うという現象を扱う理論です。

このアニメーションの上側では「指数分布に従う乱数を50個生成する」を繰り返しています。下側では、その50個の最大値を追加してヒストグラムを描いています。指数分布の最大値は、ガンベル分布に分布収束します。

最大値や最小値の分布を考えることは有益です。単に今まで得られたデータの最大値だけを見ても、今後それを上回る値が観測される可能性がどの程度あるのかも、どの程度上回るリスクがあるのかも分かりません。平均的な挙動ではなく、最大値を超えないように設定しなければならないタスク（津波を防ぐ防壁の高さの設定など）では重宝しますし、オリンピック記録の最大値の分布を調べることで、人間の限界を推定できれば面白そうではありませんか？

本記事ではそんな極値統計学の基礎理論を学んだ上で、いくつかのオープンデータに対して適用します。

redteaの思う極値統計学の目的

個人的な解釈ですが^[1]、極値統計学は将来のリスク評価です。今後起こり得る最悪の事象が起きるとすればどんな規模で起きるのか、何年以内を目安に起きるのかを知りたいモチベーションです。

基本的に過去に事例がないことを予測するのは難しいですし、将来のことは誰にも分かりません。そんな不確実性の中であっても、リスク評価くらいはしたいです。過去データでは最大30mの津波が来たという記録があったとき、30mの壁ではちょっと不安ですよね。その倍の60mあれば枕を高くして眠れますか?こういったことを、最大値の分布を調べることで少しでも精緻にリスク評価したいというモチベーションがあります。応用的には、ネットワーク遅延の最大値や、並列システムの稼働時間、機器の故障までの時間など、応用範囲は実に広いです。

ところで、わざわざ最大値の分布を調べなくとも、元のデータの分布から、0.9999999分位点とかを見ればよくないかと思う方もいらっしゃるかもしれません。しかし、基本的に分布の尾部のごくわずかなデータは外れ値です。その他多くのデータ点によって推定された分布の末尾にはほとんど情報がなく、推定誤差が大きくなります。

というわけで、元のデータ分布の端から端までを調べ尽くすのではなく、最大値の分布を調べて、その0.95分位点を見るアプローチを取ります。

統計の復習

分布関数の計算に慣れている方はこの章を飛ばしてください。極値統計学における計算は複雑になりがちなので、計算までしっかり理解したい方向けに復習から入ります。

確率密度関数

ここでは連続型関数を扱います。確率密度関数$f(x)$とは、xという値の出やすさを表します。例えば一様分布でしたら$f(x)$の値はxに依存せずに一定ですし、正規分布でしたら平均を中心に左右対称のベル型曲線を描きます。

正規分布と一様分布の確率密度関数

若干の注意点としては、確率密度関数は定義域全体で積分すると1になります^[2]が、$f(x)>1$な点も存在し得るということです。連続型関数の場合、特定のxの値を取る確率は0で、ある区間$[x, x+h]$の間に取る確率しか扱えません。そういうわけで、瞬時的に$f(x)>1$であっても、$\int_x^{x+h} f(t) dt <= 1$であるので、矛盾は起きません。

（累積）分布関数

分布関数の定義は、「確率変数$X$がある値以下（$x$）の値となる確率（$X \leq x$）」を表す関数です。よく大文字で表されるので、本記事でも$F(x)$と表記します。

$$F(x) = P[X \leq x]$$

分布関数を描くと以下のようになります。$\lim_{x \rightarrow \infin} F(x) = 1$と$\lim_{x \rightarrow 0} F(x) = 0$が成り立ちます^[3]（定義から明らかですね）。

正規分布と一様分布の分布関数

分布関数の計算

分布関数は結構簡単に計算できます。例えば、$F_X(x) = P[X \leq x]$が分かっているとき、$Y=X^2$の分布関数$F_Y(x)$は以下のように考えられます。

$$\begin{aligned} F_Y(y) &= P[Y \leq y] \\ &= P[X^2 \leq y] \\ &= P[-\sqrt{y} \leq X \leq \sqrt{y}] \\ &= F_X(\sqrt{y}) - F_X(-\sqrt{y}) \\ \end{aligned}$$

確率密度関数が知りたければ、この分布関数を$y$で微分すれば良いです。

今回は変数変換を例にしましたが、同様に和や差の分布も考えられます。平均や総和に関してはデータ数が一定あれば中心極限定理で正規分布近づいていきます。

極値統計学の考え方

ここからが本題で、いよいよ極値分布を求めていきます。元の確率変数$X$に対して、その最小値の分布や最大値の分布を考えます。最小値を考えるときは、元のデータに$-1$を掛けた上で最大値の分布を求めればよいので、本記事では最大値の分布のみを考えます。

最大値の分布を定義通りに計算する

$n$個のiid^[4]なデータ$X$が得られているとき、その最大値$M=\max{\{X_1, \ldots X_n\}}$の分布は、

$$\begin{aligned} F_M(m) &= P[M \leq m] \\ &= P[\max{\{X_1, \ldots X_n\}} \leq m] \\ &= P[X_1 \leq m, \cdots, X_n \leq m] \\ &= P[X_1 \leq m] \cdots P[X_n \leq m] \\ &= \{F_X(m)\}^n \\ \end{aligned}$$

ここで、

• 2行目から3行目にかけては、$X$の最大値が$m$以下である必要があるので、$X_1, \ldots X_n$の全てが$m$以下であれば良いという風に書き換えています

• 3行目から4行目にかけては、データが独立なので、確率の積に書き換えています

意外と早く分布関数が求まってしまいました。さらに$m$で微分すると

$$\frac{d}{dm}F_M(m) = n\{F_X(m)\}^{n-1}f_X(m)$$

が得られるので、こちらも元の確率密度関数や分布関数が分かれば書き下せます。一旦これで極値分布の主な考え方は終わりです。ただ、実はこれだけではダメで、もう一工夫が必要です。それは次節で説明します。

さらに後で詳しく解説しますが、このようにして最大値の分布が求められる場合^[5]、求めた分布は、元の確率分布に対応して、3種類いずれかの分布に、分布収束することが知られています。平均値の分布は特殊な場合を除いて正規分布という1つの分布に分布収束しますが、極値分布の場合は収束先が1つではない点に注意が必要です。

分布の退化

退化の問題

前節で、最大値の分布関数が$\{F_X(m)\}^n$と求まりましたが、これは不完全です。$n$回試行して得られた最大値を$M_n$としたとき、その分布関数の極限は以下のようになります。

$$\lim_{n \rightarrow \infin} P[M_n \leq x] = \begin{cases} 1 & (x \ge c) \\ 0 & (x < c) \end{cases}$$

上限よりも少しでも小さい値$x$を選ぶと、無限回試行すれば$M_n$は必ずその$x$を超えてしまいます。最大値が上限$c$を超えることはないので、$x \ge c$のときは分布関数の確率は1になります。「サイコロの出目の最大値の分布」を考えると直観的に理解できます。この時、上限は$c=6$です。サイコロを$n=1$回振ったときの最大値が$x=3$以下である確率はそこそこあります。そして、サイコロを$n=6$回振ったときの最大値が$x=3$以下である確率は$(3/6)^6 = 1/64$です。このように、サイコロを振る回数$n$が大きくなればなるほど、最大値が$x=3$以下になる確率は0に近づきます（$n \rightarrow \infin$となると、最大値が$x=3$以下になる確率は0でしょう）。

このように、確率変数$X$がある値$x_0$をとる確率が1になってしまうとき、$X$と対応する分布は退化していると言います。こうなってしまうと確率の計算ができないので役に立ちません。そこで、極値統計量に位置尺度変換^[6]を施すことで、退化しない分布に収束するように調整してやります^[7][8]。

退化しないための変換

少々天下り的^[9]ですが、共通の分布関数$F$を持つiid確率変数$X$について、正数列$a_n>0$と実数列$b_n$と、退化していない分布$G(x)$が存在する時、$G(x)$を極値分布と呼びます。

$$\begin{aligned} P \left[\frac{M_n - b_n}{a_n} < x \right] &= P[M_n \leq a_n x + b_n] \\ &= F^{n}(a_n x + b_n) \\ &\rightarrow G(x) \end{aligned}$$

いきなり$a_n,b_n$で最大値を変換すれば退化しないと言われてもピンとこないかもしれません。私の解釈ですが、$n$が大きくなるにつれて、取り得る最大値は大きくなったり、安定性が変わったりしそうですよね。それを$n$に関係する定数で引いたり（$-b_n$）、スケールを調整（$\frac{1}{a_n}$）したりすることで、1点に退化しないようにしているイメージです。例えば0から1の間を取る連続^[10]一様分布ですと、$n$が増えるにつれて最大値は1に限りなく近づいていきます。その1に近づく速度と同じくらい尺度を引き伸ばせば、1点に収束せずに分布に収束してくれるようなイメージです。なので、↓の方に一様分布の例も書いていますが、$a_n = 1/n$にして、最大値の分布を引き延ばしています。

収束先$G(x)$は3通りしかありません。そして、元の分布$F$に応じて定数$a_n, b_n$は決まってくれます（後段で説明がありますが、実務上$a_n, b_n$を求める必要はないのでご安心ください）。定数$a_n, b_n$の決め方は一意ではありませんが、どう選んでも収束先$G(x)$の形は同じになります。ちょっと例を見てみましょう。

指数分布

期待値が1の指数分布^[11]を考えます。確率密度関数は$f(x)=e^{-x} (x \geq 0)$です（$x < 0$ のとき$f(x)=0$）。以下、$x < 0$のときは省略します。分布関数は

$$F(x) = 1 - e^{-x}$$

で与えられるので、例えば$a_n=1, b_n=\log n$とすると、

$$\begin{aligned} P \left[ \frac{M_n - b_n}{a_n} \leq x \right] &= F^n(a_n x + b_n) \\ &= F^n(x+ \log n) \\ &= \left(1- e^{-(x+\log n)} \right)^n \\ &= \left( 1 + \frac{-e^{-x}}{n} \right)^n \\ &\rightarrow e^{-e^{-x}} \end{aligned}$$

また、同様の指数分布の設定で、例えば$a_n=2, b_n=\log 3n$とすると、

$$F^n(2x + \log 3n) \rightarrow e^{-\frac{1}{3}e^{-2x}}$$

となり、別の分布になりますが、いずれもガンベル分布と言われる分布に収束します。このように、定数列$a_n, b_n$の決め方は一意ではありませんが、パラメータの異なる同じ極値分布$G_1, G_2$について、$G_1(ax + b) = G_2(x)$が成り立つので（証明略）、この違いは問題になりません。適当に収束する定数列$a_n, b_n$が1つ見つかれば良いです。

以下、いくつか有名な確率分布の極値分布の例を置いておきます。

正規分布

正規分布

標準正規分布 $N(0,1)$ を考えます。確率密度関数は
$f(x)=\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$、分布関数を $F$ とします。
このとき上側の裾確率について

$$1-F(x) \sim \frac{1}{x\sqrt{2\pi}}e^{-x^2/2}\qquad (x\to\infty)$$

が成り立つので、$b_n$ を

$$n\{1-F(b_n)\}=1$$

を満たすように取ると、近似的に

$$b_n = \sqrt{2\log n} -\frac{\log\log n+\log(4\pi)}{2\sqrt{2\log n}} \quad(\text{より高次は省略})$$

となります。また $a_n=1/b_n$ とおくと、任意の $x\in\mathbb{R}$ に対し

$$\begin{aligned} P\left[\frac{M_n-b_n}{a_n}\le x\right] &=F^n(a_n x+b_n)\\ &=\left\{1-(1-F(a_nx+b_n))\right\}^n\\ &\sim \left\{1-\frac{1}{n}e^{-x}\right\}^n\\ &\rightarrow e^{-e^{-x}} \end{aligned}$$

となり、正規分布の最大値もガンベル分布に収束します。

パレート分布

パレート分布

裾の重い例として、形状母数 $\alpha>0$ のパレート分布を考えます。
分布関数は

$$F(x)=1-x^{-\alpha}\qquad (x\ge 1)$$

で与えられます（$x<1$ のとき $F(x)=0$）。
ここで $a_n=n^{1/\alpha},\ b_n=0$ とおくと

$$\begin{aligned} P\left[\frac{M_n-b_n}{a_n}\le x\right] &=F^n(a_n x)\\ &=\left(1-(a_nx)^{-\alpha}\right)^n\\ &=\left(1-\frac{x^{-\alpha}}{n}\right)^n\\ &\rightarrow e^{-x^{-\alpha}} \qquad (x>0) \end{aligned}$$

となります。極限分布 $e^{-x^{-\alpha}}$ はフレシェ分布と呼ばれ、
指数分布や正規分布とは別の（裾の重さを反映した）型に収束します。
なお $x\le 0$ のときは確率 $0$ です。

一様分布

一様分布

有界台の例として、一様分布 $\mathrm{Unif}(0,1)$ を考えます。
分布関数は

$$F(x)=x\qquad (0\le x\le 1)$$

で与えられます（$x<0$ のとき $F(x)=0,\ x>1$ のとき $F(x)=1$）。
ここで上側端点が $1$ なので、$a_n=1/n,\ b_n=1$ とおくと

$$\begin{aligned} P\left[\frac{M_n-b_n}{a_n}\le x\right] &=F^n\!\left(1+\frac{x}{n}\right)\\ &=\left(1+\frac{x}{n}\right)^n\\ &\rightarrow e^{x} \qquad (x\le 0) \end{aligned}$$

となります（$x>0$ では $1+x/n>1$ となり $F=1$ なので極限は $1$）。
したがって極限分布は

$$G(x)= \begin{cases} e^{x}, & x\le 0,\\ 1, & x>0, \end{cases}$$

で、これはワイブル（有界端点型）の極値分布です。

極値分布の一般化

実は、極値分布の分布関数は3種類しかありません。元の分布が決まれば、分布収束先も1つに決まります。さらに実務上においては、元の分布を気にすることなく、1つの一般的な分布のパラメータを推定することで、どの形の極値分布か判定できます。本節では、まずは3種類の極値分布を見て、それぞれの特徴を理解します。その上で、一般化された極値分布が3種類の極値分布とどのように対応するかを確認します。

極値分布の3つの型

極値分布の分布関数は、以下3つの型のいずれかになります。それぞれ順番にタイプ1、タイプ2、タイプ3の極値分布と呼ばれています。

• ガンベル分布（Gumbel distribution）
• フレシェ分布（Fréchet distribution）
• ワイブル分布（Weibull distribution）

それぞれの確率密度関数の概形は以下の図の通りです。

3種類の極値分布の確率密度関数の概形。ワイブル分布（青線）では2以上の値を取る確率が0になっています。これは大きい値を取るリスクが小さいときに現れます。ガンベル分布（黒点線）は減衰が速いので、大きな値を取る確率は小さいです。最後にフレシェ分布（赤線）ですが、右側の裾が厚く、大きな値を取る確率がそこそこあります。なお、凡例にある$\xi$の意味は次章「一般極値分布」で登場します。

フレシェ > グンベル > ワイブル の順に大きな値を取る確率が高く、リスクが高いと考えられます。例えば水位や気温の最大値がどの分布の形になっているかを知ることは、リスク評価の上で非常に重要です。

3種類の極値分布の分布関数以下の通りです。確率密度関数は複雑になるうえ、本記事では使わないので割愛します。気になる方はwikipediaを閲覧するか、自力で微分してください。

ガンベル分布

$$F(x) = \exp\left( - \exp\left( - \frac{x - \mu}{\sigma} \right) \right)$$

フレシェ分布

$$F(x) = \exp\left( - \left( \frac{x - \mu}{\sigma} \right)^{-\alpha} \right)$$

ワイブル分布

$$F(x) = 1 - \exp\left( - \left( \frac{x - \mu}{\sigma} \right)^{\alpha} \right)$$

一般極値分布

先ほど紹介した3種の極値分布は1つの式、一般極値分布（Generalized Extreme Value Distribution）で表せます。この1つの一般極値分布のパラメータをデータから推定して、その形状パラメータから、3種のどの分布に対応するかを考えられます。確率密度関数と分布関数は以下の通りです。

$$f(x) = \frac{1}{\sigma} \left[ 1 + \xi \left( \frac{x - \mu}{\sigma} \right) \right]^{(-1/\xi) - 1} \exp\left\{ - \left[ 1 + \xi \left( \frac{x - \mu}{\sigma} \right) \right]^{-1/\xi} \right\}$$

$$F(x) = \exp\left\{ - \left[ 1 + \xi \left( \frac{x - \mu}{\sigma} \right) \right]^{-1/\xi} \right\}$$

一般極値分布の形状パラメータ$\xi$と3つの極値分布の対応は以下の通りです。

形状パラメータ $\xi$	対応する極値分布の型	分布関数の形	裾（テール）/サポートの特徴
$\xi = 0$	ガンベル分布（Type I）	$F(x)=\exp\left\{-\exp\left(-\frac{x-\mu}{\sigma}\right)\right\}$	右裾は指数的に減衰。上限なし。
$\xi > 0$	フレシェ分布（Type II）	$F(x)=\exp\left\{-\left(1+\xi\frac{x-\mu}{\sigma}\right)^{-1/\xi}\right\}$	右裾が厚い（heavy tail）。上限なし。
$\xi < 0$	ワイブル分布（Type III）	$F(x)=\exp\left\{-\left(1+\xi\frac{x-\mu}{\sigma}\right)^{-1/\xi}\right\}$	上限あり（有限端点）。$x < \mu-\sigma/\xi$ の範囲で定義。

形状パラメータ$\xi$の値に応じた一般極値分布の確率密度関数は以下のような概形になります。

形状パラメータ$\xi$に応じた一般極値分布の確率密度関数の概形。loc=0, scale=1で指定（$\mu=0, \sigma=1$）。

以下、少々数式が続きますが、分布関数に着目して、一般極値分布が3つの極値分布と一致することを確認していきます。

$\xi$の値に応じて3つの極値分布と一致し、↑の表のような対応関係になることを信じてくださる方は次の計算は読み飛ばしても差し支えございません。

ガンベル分布

$\xi=0$のとき、一般極値分布はガンベル分布になります。分布関数

$$F(x) = \exp\left\{ - \left[ 1 + \xi \left( \frac{x - \mu}{\sigma} \right) \right]^{-1/\xi} \right\}$$

において、$z=\frac{x-\mu}{\sigma}$とおくと

$$F(x)=\exp\left\{-(1+\xi z)^{-1/\xi}\right\}$$

です。ここで$\xi\to 0$の極限を取ります。

まず、

$$(1+\xi z)^{-1/\xi} = \exp\left\{-\frac{1}{\xi}\log(1+\xi z)\right\}$$

と書けます。$\xi\to 0$で$\log(1+\xi z)=\xi z+o(\xi)$より、

$$-\frac{1}{\xi}\log(1+\xi z) = -\frac{1}{\xi}(\xi z+o(\xi)) \to -z$$

したがって

$$(1+\xi z)^{-1/\xi}\to \exp(-z)$$

となり、

$$\begin{aligned} F(x) &= \exp\left\{-(1+\xi z)^{-1/\xi}\right\} \\ &\to \exp\left\{-\exp(-z)\right\} \\ &= \exp\left\{ - \exp\left( -\frac{x-\mu}{\sigma} \right) \right\} \end{aligned}$$

を得ます。よって$\xi=0$の一般極値分布はガンベル分布と一致します。

フレシェ分布

$\xi>0$のとき、フレシェ分布$\exp\left( - \left( \frac{x - \mu}{\sigma} \right)^{-\alpha} \right)$に対し、$\mu=-\frac{1}{\xi}, \sigma=\frac{1}{\xi}, \alpha=\frac{1}{\xi}$とおくと、

$$\begin{aligned} F(x) &= \exp\left( - \left( \frac{x - \mu}{\sigma} \right)^{-\alpha} \right) \\ &= \exp\left( - \left( \frac{x + \frac{1}{\xi}}{\frac{1}{\xi}} \right)^{-\frac{1}{\xi}} \right) \\ &= \exp\left\{ - \left( 1 + \xi x \right)^{-1/\xi} \right\}, x > -\frac{1}{\xi} \end{aligned}$$

となるので、これは一般極値分布の位置尺度変換（平均$\mu$を引いて、スケール$\sigma$で割る）したものと一致しています。

ワイブル分布

$\xi<0$のとき、（上限を持つ）ワイブル分布の分布関数

$$F(x) = \exp\left( - \left( \frac{\mu - x}{\sigma} \right)^{\alpha} \right), \quad x \le \mu$$

に対し、$\mu=-\frac{1}{\xi}, \sigma=-\frac{1}{\xi}, \alpha=-\frac{1}{\xi}$とおくと（$\xi<0$なので$\sigma>0,\alpha>0$）、

$$\begin{aligned} F(x) &= \exp\left( - \left( \frac{\mu - x}{\sigma} \right)^{\alpha} \right) \\ &= \exp\left( - \left( \frac{-\frac{1}{\xi} - x}{-\frac{1}{\xi}} \right)^{-\frac{1}{\xi}} \right) \\ &= \exp\left\{ - \left( 1 + \xi x \right)^{-1/\xi} \right\}, \quad x < -\frac{1}{\xi} \end{aligned}$$

となるので、これは一般極値分布の位置尺度変換（平均$\mu$を引いて、スケール$\sigma$で割る）したものと一致しています。

データ取得方法と極値分布の推定法

一般極値分布のパラメータを推定するためのデータがどのように取得されるかによって、一般極値分布の推定法が変わります。データ取得方法は主に以下の3つです。

• ブロック最大値
• 上位$r$個データ
• 閾値超過データ

ブロック最大値

データ取得方法

ブロック最大値とは、例えば100年分、毎日の降水量のデータが蓄積されているとします。パラメータの推定にはいくつかの最大値のデータが欲しいので、例えば1年ごとにブロックを区切って、年間の最大値を100個求めます。一般的な言い方をすると、全データを$K$個のブロックに区切り、それぞれの最大値$m_1,\cdots m_K$を得て、これらのデータで一般極値分布のパラメータを最尤推定します。

ブロック最大値のイメージ。適当な大きさ（年間とか1日とか、データに合わせて決めます）でブロックを区切り、その中の最大値を用いるシンプルな方法です。

場合によってはお手元のデータが既にブロック最大値になっている可能性もあります。例えば日々の最高気温データがあれば、それは日をブロックとしたブロック最大値データです。

なお、ブロック中のデータにいくつか欠損があっても基本的に気にしなくても良いです。最大値が欠損するような状況であれば話は別ですが、大半のデータを捨てることになるので、多少欠損があっても問題なしです。

最後に、$K$の選び方です。100年分のデータを1年ごとに区切って100個の最大値を作るのか、10年ごとに区切って10個の最大値を作るのか、という話です。これは完全なトレードオフになってしまうので、一般的な良い方法ありません。$K$を小さくすればブロック内でしっかりと最大値が取れるので、極値分布への収束がよくなります。一方、粗く分割するので、$K$個しかデータを取得できず、推定量のばらつきが大きくなってしまいます。$K$を大きくすると逆のことが起きます。

パラメータ推定法

ブロック最大値で得たデータに対して一般極値分布を推定するには最尤法を用いれば良いです。対数尤度関数は以下の通りです。

$$z_i=\frac{x_i-\mu}{\sigma},\qquad t_i=1+\xi z_i=1+\xi\left(\frac{x_i-\mu}{\sigma}\right)$$

$$L(\mu,\sigma,\xi)=\prod_{i=1}^K \left[ \frac{1}{\sigma}\, t_i^{-1/\xi-1}\, \exp\!\left(-t_i^{-1/\xi}\right) \right]$$

$$\ell(\mu,\sigma,\xi)=\log L(\mu,\sigma,\xi) = -K\log\sigma +\left(-\frac{1}{\xi}-1\right)\sum_{i=1}^K \log t_i -\sum_{i=1}^K t_i^{-1/\xi}$$

$$\text{（制約）}\quad t_i>0\ \ (\forall i),\ \sigma>0$$

ただ、この偏微分=0は解析的に解けないので、ニュートン法などを用います。

上位$r$個データ

データ取得方法

先ほどのブロック最大値では、各ブロックごとの最大値を用いていましたが、「上位$r$個データ」では、最大値だけでなく、上位$r$件使います。最大値だけではブロックを区切る位置によって推定結果が変わってしまったり、十分なデータ件数を取得できなかったりする際に効果的です。

上位$r$個データのイメージ。適当な大きさ（年間とか1日とか、データに合わせて決めます）でブロックを区切り、その中の上位r件（この図は$r=3$）を用いる方法です。オリンピックの銅メダル、銀メダル、金メダルのスコアを記録するようなイメージです。

パラメータ推定法

上位$r$件ずつ取るので、同時密度関数が変わります。各ブロックに$n$個のデータがあり、それらの確率変数列$X_1, \ldots X_n$の順序統計量を

$$X_{(n)} \geq X_{(n-1)} \ldots \geq X_{(1)}$$

とします。

最大値の分布は既にみてきたように、

$$\frac{X_{(n)} - b_n}{a_n} \sim G_{\xi}$$

と、極値分布に従います。同じ定数列$a_n, b_n$を用いて、上位$r$個の順序統計量を標準化した$z_1, \ldots z_r$の同時分布関数の極限を求めると、

$$G_{\xi}(z_1, \ldots z_r) = \lim_{n \rightarrow \infin} P \left( \frac{X_{(n)} - b_n}{a_n} \leq z_1, \ldots , \frac{X_{(n-r+1)} - b_n}{a_n} \leq z_r \right)$$

と表せます。この極限分布の同時密度関数の計算は大変なので省略しますが、結果は以下の通りです。

$$g_{\xi}(z_1, \ldots z_r) = \begin{cases} \dfrac{1}{\sigma^r} \exp\!\left[ -\left\{ 1+\xi\left(\dfrac{z_r-\mu}{\sigma}\right) \right\}^{-\frac{1}{\xi}} \right] \displaystyle\prod_{i=1}^{r} \left\{ 1+\xi\left(\dfrac{z_i-\mu}{\sigma}\right) \right\}^{-\frac{1}{\xi}-1}, & \xi \neq 0, \\[1.5em] \dfrac{1}{\sigma^r} \exp\!\left[ -\exp\!\left\{ -\left(\dfrac{z_r-\mu}{\sigma}\right) \right\} \right] \displaystyle\prod_{i=1}^{r} \exp\!\left\{ -\left(\dfrac{z_i-\mu}{\sigma}\right) \right\}, & \xi = 0 \end{cases}$$

これも対数尤度=0は解けないので、ニュートン法によって数値計算的に解くことになります。

閾値超過データ

データ取得方法

閾値超過法（Peak Over Threshold）では、ある閾値を超えたデータのみを用いて極値分布を推定します。ある閾値を超えた「異常データ」のみが蓄積されているケースもあり得ますし、風速や地震などが一定値以上となる場合、大きな問題になると考えれば、このようなデータの取得方法が自然とも思えます^[12]。

閾値超過データのイメージ。適切な閾値を設定して、その値を超えたデータを用いる方法です。

パラメータ推定法

閾値を$u$とします。すなわち、$u$以上のデータ$X$のみが記録されている（そのデータのみ使う）とします。データは閾値$u$を超えているという条件下でのみ観測されるので、分布関数は条件付き確率の式で以下のように書けます。

$$\begin{aligned} F_u(y) &= P(X - u \leq y | X > u) \\ &= \frac{F(y+u) - F(u)}{1-F(u)} \\ &= 1 - \frac{1-F(y+u)}{1-F(u)}, \quad y \geq 0 \end{aligned}$$

なお、$F(x)$は$X$が従う分布の分布関数です。ここで、iidな確率変数$X_1, \ldots, X_n$のうち、閾値$u$よりも大きいデータが$X_{(1)}, \ldots, X_{(m)}$の$m$個であるとき、閾値との差を

$$\begin{aligned} Y_1 &= X_{(1)}-u \\ Y_2 &= X_{(2)}-u \\ \cdots \\ Y_m &= X_{(m)}-u \\ \end{aligned}$$

とおくと、これらは（$X$が独立なので）独立であり、$u$が十分に大きければ漸近的に一般化パレート分布に従うことが知られています。一般パレート分布の分布関数$F(x)$と確率密度関数$f(x)$は以下の通りです。

$$F(z)=1-\left(1+\xi z\right)^{-1/\xi} \\ f(x)=\frac{1}{\sigma}\left(1+\xi\frac{x-\mu}{\sigma}\right)^{-\left(\frac{1}{\xi}+1\right)}$$

今、$Y = X - u$で、$z = \frac{X - u}{\sigma}$と考えれば、一般化パレート分布の定義と符合します。ここで重要なことは、一般極値分布のパラメータ$\xi$と、この式の一般化パレート分布のパラメータ$\xi$は同じです。これによって、十分に大きな$u$を選べば、一般化パレート分布の最尤推定よって極値分布を求めることができます。ただ、例のごとく解けないので数値計算的に解くことになります。

なお、閾値$u$の決め方ですが、これもブロック最大値の$K$の決め方同様トレードオフの関係です。$u$が大きい方がパレート分布へ分布収束しやすいですが、データ件数が減ってバリアンスが高くなり、$u$が小さいと逆のことが起きます。

再現レベル

「100年に一度の出来事」とはどんなレベルの事象なのか、を考えるのが再現レベルです。ある確率$p$に対して、一般極値分布の上側確率点は、

$$1-p = \exp \left[ - \left( 1+\xi \frac{z_p - \mu}{\sigma} \right)^{-1/\xi} \right]$$

を満たす点$z_p$です。これは$z_p$について解けます。

$$z_p = \begin{cases} \mu + \frac{\sigma}{\xi} \{(-\log(1-p))^{-\xi} -1\} & (\xi \neq 0) \\ \mu - \sigma \log(-\log(1-p)) & (\xi \rightarrow 0) \end{cases}$$

この確率点$z_p$を、再現期間$\frac{1}{p}$の再現レベルと呼びます。年間最大値のデータを用いて推定した極値分布に対して、$p=1/100$を設定すれば、再現期間100年の再現レベル、すなわち100年に一度観測されるような値が分かります。

時系列データにおける極限分布

本記事では既にブロック最大値のような方法で、時系列データに対する極値分布を求めました。実際、時系列データに対してもブロック最大値のような考え方を適用すれば良いのですが、iidであるという仮定を置いていました。時系列データの多くはiidを満たさないものが多いです。例えば最高気温は周期性がありますし、前日との自己相関も強そうです。また、オリンピックの記録はトレーニング方法や道具の改良によって年々向上しています。

ミキシング性

時系列データに対する分析は通常、定常性を仮定する場合が多いですが、極値統計学ではミキシング性を仮定した理論が作られています。ミキシングの定義は難しいので、本記事では雰囲気で書きます。時系列$\{X_i\}$に対して、$\mathcal{M_1} = \{X_1, X_2, \ldots X_p\}$と$\mathcal{M_2} = \{X_{p+l}, X_{p+l+1}, \ldots\}$という風に、$l$期離れた2つの系列を考えます。そして、$l$が十分大きいとき（2つの系列が十分離れているとき）、$\mathcal{M_1},\mathcal{M_2}$が独立であることをミキシング性がある言います。

時系列データへの適用

実は多少自己相関があってもブロック最大値の考え方は、ブロックがある程度大きければミキシング性が成り立つと考えられるので、そのまま適用できます^[13]。また、トレンドや周期性は事前にその成分を取り除いた系列を時系列モデルから作ったり、差分系列を取ることで定常過程にすれば極値統計学の理論をそのまま適用できます。もちろん、差分系列に対する極値分布は、「差分の最大値の分布」になるので、大きな変化がどれくらい起こるかを示した分布であることに注意です。実務上、トレンドなどがあるデータに対して差分の最大値を知ることは有用です。

Pythonによる極値統計学の適用

R, Python共に極値分布を扱うライブラリは整備されています。今回はPythonのScipyを使って、種々の極値分布を推定します。

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.genextreme.html

PythonのScipyでは、形状パラメータを$c=-\xi$で定義しています。したがって、本記事で紹介した一般極値分布のパラメータ$\xi$と正負の対応関係が逆なので要注意です。

なお、Scipyは極値統計学の一部しか扱えず、ブロック最大値くらいしかGEV分布を推定できません。上位$r$個データや、閾値超過データに対する分布推定をするにはpyextremesという専用ライブラリがお勧めです。pyextremesも、パラメータ$\xi$と形状パラメータcは$c=-\xi$の関係になっているので要注意です。pyextremesを使うのも簡単で、生成AIがある程度書いてくれるはずなので、本記事ではpyextremesの実装は割愛します。

https://georgebv.github.io/pyextremes/

データの取得

気象庁の「過去の気象データ・ダウンロード」サイトからデータをダウンロードできます。

https://www.data.jma.go.jp/risk/obsdl/

以下の図のように、各地域の最大風速、最高気温、最大降水量などのデータを取得できます。

気象庁より。データ取得時の設定。

ここでダウンロードしたデータを適当な場所において（data/data.csv）、次項のPythonスクリプトを実行してみてください。

なお、地球の環境変化などによって、最高気温の時系列データが定常性があるとは言えないかもしれませんが、本記事では定常性があると思い込んで分析します。

一般極値分布の推定

推定はライブラリにやらせましょう。コア部分は以下の通りです。maximaにはarray-likeなデータならなんでもOKです。最大値$\{m_1, \cdots m_K\}$を入れておけばOKです。このshapeが負ならフレシェ分布、正ならワイブル分布、0ならガンベル分布です（繰り返しになりますが、本記事の$\xi$とscipyのshapeは正負が逆です）。

from scipy.stats import genextreme
shape, loc, scale = genextreme.fit(maxima)

パラメータ推定した結果を以下に掲載します。当てはまりの良さを確認することは重要です。ヒストグラムと理論分布を重ねるのでも良いし、Q-Qプロットを見るのも良いでしょう。いずれの場合でも、明確な基準は存在しません。実務では第一種過誤と第二種過誤に注意しながら、数値やグラフだけでなく対象となるデータの背景も加味して判断する必要があります。

最大風速の極値分布です。パラメータshapeが0に近く、ガンベル分布っぽいです。低い確率で、今まで観測したことのない25(m/s)を超える風速が観測される可能性があります。

降水量の極値分布です。こちらもshapeパラメータが0に近いですが、負なので、フレシェ分布よりです。グラフを見ると裾が厚いようには見えませんが、低い確率で大きな値を取り得ます。

最高気温の極値分布です。shapeパラメータが正なので、大きな値を取りにくいです。最高気温が突拍子もない値が出ないだろうという直感と符合します。

from extreme.bootstrap_ci import bootstrap_ciのパス設定さえしっかりすれば、そのまま動く（はず）のPythonスクリプト全体を以下に置いておきます。

from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import polars as pl
from scipy.stats import genextreme

from extreme.bootstrap_ci import bootstrap_ci

SYUIRO = "#F26649"
BLUE = "#1E90FF"
# 日本語フォントの設定
plt.rcParams["font.family"] = "MS Gothic"  # Windows環境用


def main():
    path = Path("data/data.csv")
    df = pl.read_csv(path, encoding="shift_jis")
    # 分析対象を以下から選ぶ
    # target = "日降水量の最大(mm)"
    target = "最大風速(m/s)"
    # target = "最高気温(℃)"

    # 年月列から日付を抽出し、年を取得
    df = (
        df.select(["年月", target])
        .with_columns(pl.col("年月").str.strip_chars().str.split("/").alias("ym_parts"))
        .with_columns(
            pl.col("ym_parts").list.get(0).cast(pl.Int32).alias("year"),
            pl.col("ym_parts").list.get(1).cast(pl.Int8).alias("month"),
        )
        .with_columns(pl.col(target).cast(pl.Float64, strict=False))
        .drop_nulls(subset=["year", "month", target])
        .drop("ym_parts")
    )

    # 年ごとの最高気温を集計
    block_maxima = (
        df.group_by("year").agg(pl.col(target).max().alias("annual_max")).sort("year")
    )

    maxima = block_maxima.get_column("annual_max").to_numpy()
    if maxima.size == 0:
        msg = "No annual maxima available after filtering; cannot fit GEV."
        raise ValueError(msg)

    shape, loc, scale = genextreme.fit(maxima)

    print(
        "GEV parameters (shape, location, scale): "
        f"({shape:.4f}, {loc:.2f}, {scale:.2f})"
    )

    # ブートストラップ法で信頼区間を計算
    ci_lower, ci_upper = bootstrap_ci(maxima)

    # 再現レベルを計算して表示
    return_periods = [10, 20, 50, 100, 200]
    print("\n再現レベル:")
    for period in return_periods:
        return_level = genextreme.ppf(1 - 1 / period, shape, loc=loc, scale=scale)
        print(f"{period}年再現レベル: {return_level:.2f}")

    plot_gev_fitting(maxima, shape, loc, scale, target)


def plot_gev_fitting(
    maxima: np.ndarray, shape: float, loc: float, scale: float, target: str
) -> None:
    # プロットの作成
    span = float(maxima.max() - maxima.min())
    padding = span * 0.1 if span > 0 else max(abs(float(maxima.max())) * 0.1, 1.0)
    x_min = float(maxima.min()) - padding
    x_max = float(maxima.max()) + padding
    x_values = np.linspace(x_min, x_max, 400)
    pdf_values = genextreme.pdf(x_values, shape, loc=loc, scale=scale)

    fig, ax = plt.subplots(figsize=(10, 6))
    ax.hist(
        maxima,
        density=True,
        alpha=0.5,
        label="年間最大値",
        color="lightgray",
        edgecolor="black",
    )
    ax.plot(
        x_values, pdf_values, color=SYUIRO, label="推定した一般極値分布の確率密度関数"
    )
    ax.set_xlabel(f"{target}の年間最大値")
    ax.set_ylabel("Density")
    ax.set_title(f"{target}の年間最大値のGEVフィッティング")
    ax.legend()

    # 推定したパラメータをテキストで図中に表示
    param_text = (
        f"Shape: {shape:.4f}\nLocation: {loc:.2f}\nScale: {scale:.2f}"
    )
    ax.text(
        0.95,
        0.95,
        param_text,
        transform=ax.transAxes,
        fontsize=10,
        verticalalignment="top",
        horizontalalignment="right",
        bbox=dict(boxstyle="round", facecolor="white", alpha=0.8),
    )
    fig.tight_layout()
    plt.show()

    # 保存用にファイル名を整形
    save_target = (
        target.replace("(", "")
        .replace(")", "")
        .replace("／", "_")
        .replace("/", "_")
        .replace(" ", "_")
    )
    fig.savefig(f"gev_fitting_illustration_{save_target}.png", dpi=300)


if __name__ == "__main__":
    main()

shape=0って現れないのでは？

ところで、shapeがちょうど0になることってほとんどあり得ないと思います。そうなると、ガンベル分布って現れないのでは🤔と思うかもしれませんが、これは問題ありません。最初に3種類の極値分布があるという話をしましたが、どれか1つになると考えるのではなく、「フレシェ分布よりのガンベル分布」のようなもので、分布はグラデーションのように考えるべきです。そういう意味ではパラメータは点推定ではなく、区間推定すべきです。極値分布のパラメータを区間推定する方法はいくつかありますが、実務的にはブートストラップCIを用いるのが無難と思います。このように区間推定して、形状パラメータshapeの区間が0をまたぐようであれば判断を保留にしたり、危険よりの形状でのリスク評価をしたりなど、柔軟な分析ができます。以下にサンプルコードを置いておきます。

import numpy as np
from scipy.stats import genextreme


def bootstrap_ci(maxima: np.ndarray, alpha=0.05) -> tuple[np.ndarray, np.ndarray]:
    # パラメータの信頼区間を計算（ブートストラップ法）
    n_bootstrap = 1000
    bootstrap_params = []
    rng = np.random.default_rng(seed=416)
    for _ in range(n_bootstrap):
        bootstrap_sample = rng.choice(maxima, size=len(maxima), replace=True)
        try:
            params = genextreme.fit(bootstrap_sample)
            bootstrap_params.append(params)
        except Exception:
            continue

    bootstrap_params = np.array(bootstrap_params)
    ci_lower = np.percentile(bootstrap_params, q=(alpha / 2) * 100, axis=0)
    ci_upper = np.percentile(bootstrap_params, q=(1 - alpha / 2) * 100, axis=0)
    level = (1 - alpha) * 100

    print(f"形状パラメータ: ({level:.0f}%CI: [{ci_lower[0]:.4f}, {ci_upper[0]:.4f}])")
    print(f"位置パラメータ: ({level:.0f}%CI: [{ci_lower[1]:.2f}, {ci_upper[1]:.2f}])")
    print(f"尺度パラメータ: ({level:.0f}%CI: [{ci_lower[2]:.2f}, {ci_upper[2]:.2f}])")
    return ci_lower, ci_upper

実際、先のデータで適用すると以下のようになりました。この結果から、最高気温はワイブル分布型と言って問題なさそうで、突拍子もない最高気温が記録されるリスクはかなり低そうです。一方、降水量は裾があり、フレシェ分布よりのガンベル分布かもしれないので、今後より大きな値が観測される可能性があると考えられます。最大風速は0をしっかり跨いでいるので、ガンベル分布っぽいという判断で良いでしょう。

指標	パラメータ種別	95%信頼区間 (CI)
日降水量の最大 (mm)	形状パラメータ	[-3.5278, 0.0905]
最大風速 (m/s)	形状パラメータ	[-0.1463, 0.2371]
最高気温 (℃)	形状パラメータ	[0.0637, 0.2827]

モデル診断

モデルを推定したならば、それが正しく推定されているか、定量的や定性的に確認せねばなりません。scipyで推定したものでもこの図は普通に描けますが、pyextremesのplot_diagnostic()を使うと以下のような図が一発で出てきて便利です。定性的に見て、最大風速は程よく推定できているように見えます（同じような図を何枚も貼るのもアレなので最大風速だけにしました）。

pyextremesによる最大風速(m/s)のモデル診断図です。横軸に期間、縦軸に再現レベルを取った再現レベルプロット（左上）、経験分布と推定した確率密度関数を重ねたグラフ（右上）、Q-Qプロット（左下）、P-Pプロット（右下）。

結び

なにも極値統計学に限った話ではありませんが、極値統計で得られた分布の推定結果はあくまでも現象の近似であり、適切な近似になっているかどうかに常に注意を払わねばなりません。極値統計は特に、元データのiidへの依存が強い理論ですが、実務上のデータで綺麗なiidであることは稀でしょう。厳密に仮定が満たされていない状況下であっても、多角的に考えて、極値分布の理論を適用して良いかどうかを見定めることこそが、我々データサイエンティストの重要な仕事だと思っています。なお、当てはまりの良さを考えるのは個別具体の話になる上に、方法論はQ-Qプロットなど、極値統計学に限った話ではないので、本記事では割愛しました。

断っておきますが、統計学では100年に一度の災害は予測できません。極値統計学という武器を以てしても、所詮は手元にあるデータから最大限考えているだけに過ぎず、データが得られる前提が100年間変わらない理由はどこにもないからです^[14]。これは、手元のデータがiidであったり定常性、ミキシング性があったとしても、将来もそうとは限らないと言っています。

---

『統計学を拓いた異才たち』という書籍の第6章「100年に一度の洪水」に登場する一節を紹介します。

最も厄介な問題の一つは、紡いだばかりの新しい錦糸の強度についてだった。繊維を撚り合わせた一本の錦糸が切れるまでの張力は、同じ条件で撚られていも、そのつど大きく異なっていた。ティペットは慎重に実験を重ね、異なる張力を加えた錦糸の破断面を顕微鏡で調べた。そして、糸の破断は撚った糸の中で最も弱い繊維の強度に依存することを突き止めた。
（出典：『統計学を拓いた異才たち』 第6章「100年に一度の洪水」）

この研究は、錦糸を構成する各繊維がiidであることを用いた成果です。錦糸を構成する各繊維の強度の分布は、気象などと比べると確かに将来大きく変わりにくそうですし、分布の変化も計測によって確認できます。このように、極値統計学の適用範囲は想像より広く、身に着けることで分析の幅が大きく広がります。

最後に、本記事のタイトル「100年に一度の洪水」が回収されていないことにお気づきのあなたへ。『統計学を拓いた異才たち』の第6章のタイトルが「100年に一度の洪水」であるのにも関わらず、同書の中で「100年に一度の洪水」を掘り下げた話は出てきません。本記事はその肩すかしへのオマージュです。

参考文献

• デイヴィッド・サルツブルグ, 統計学を拓いた異才たち, 日本経済新聞出版, 2010.
• 久保川達也, 現代数理統計学の基礎, 共立出版株式会社, 2017.
• 西郷 達彦, Rによる極値統計学, オーム社, 2020.

脚注

1. 極値統計学を実務で使いまくっているわけでもなければ、勉強しまくったわけでもないので、保険のために枕詞を入れています。

2. いずれかの値をとるはずなので、全ての確率が1になるのは納得ですよね。

3. $0$や$\infin$は定義域の端っこくらいの意味で書いています。突っ込まないでください。

4. independent and identically distributed. 独立同分布に従うという意味です。

5. コーシー分布に従う確率変数の和が中心極限定理にならないように、例外はあります。

6. 平行移動とスケール調整くらいの意味です。

7. 位置尺度変換以外の変換をすると、あらゆる分布が作れてしまうので、位置尺度変換に限っています。

8. サイコロの例だと、元の分布に上限があるからこのようなことが起きると思われるかもしれませんが、上限のない分布であってもこのような位置尺度変換が必要です。

9. 問題を解く上で必要な数式や概念を、根拠を明示せずに導入するさま。コトバンクより

10. 先ほど挙げたサイコロの例では残念ながら離散なので退化してしまいます。例としては不適切だったかもしれませんが、まずは退化のイメージを掴んでいただくためサイコロの例を挙げていました。

11. 指数分布の定義には流派がありますが、期待値が1なら流派に関係なしでハッピーです。

12. この文言からなんとなく一般化パレート分布を思い浮かべた方は鋭いです。続きを楽しみにお読みください。

13. クラスターのあるような系列の場合、各クラスターから1点のみサンプリングするような処理が必要な場合があります。

14. 100年間計算機上でシミュレーションするような状況であれば「予測可能」かもしれませんね。