0から始めた極値理論 - プログラムdeタマゴ

　どうも、お久しぶりです。
　さて、今回記事にする極値理論をご存じでしょうか。主に土木工学等で発展した統計理論で、21世紀になってから良解説書が相次いで出版されたことで、世界的に様々な分野で（日本を除く）ブレイク中の学問だそうです。

概要

　観測しうる全てのデータにマッチする確率モデルを生成し、それを元にした予測を行えば、万象は確率モデルの手のひらの中となる。全ては予測可能である。

　なんてことは全くなく、全てのデータに対してマッチする確率モデルなんて作れませんし、仮に作ったとしたら間違いなく過学習が起こって予測精度は悲惨なことになるだろう。ということは、利用するデータは、減らす、統合する、選択するなど何らかの人工的な操作を行ってからモデル化するよりほかない。
　極値理論はその中でも、ある区間の中での最大値（もしくは、ある区間内で閾値以上となる値）という操作を適応したデータに対して、その出現確率を求めるという話です。
　元々は土木工学の分野から発達した統計理論で、「土手が決壊しないですむ十分な高さはどのくらいか？」という問題を解く為に出発したそうです。土手が決壊するかどうかは、普段の水位を考慮してつくったって、意味がありません。重要なのは大雨等で増水した時に、水位が取り得る最大値だけです。そこで、ある区間の中での最大値を用いる統計論、というのが出てきたわけですね。

一般極値分布関数（GEV分布関数）

　さて、ある期間(長さn）の観測データ{X1,X2,...,Xn}の最大値Znについて考えてみます。

$Z_n = \max\{X_1,X_2,...,X_n\}$

　このZnが閾値xよりも小さくなる可能性について考えてみよう、というのがここでの主題です。
$P(Z_n < x) = P(\max\{X_1,X_2,...,X_n\} < x)=P(X_1 < x,X_2 < x,...,X_n < x)$
　上の式は、要するに、最大値がx以下ってことは、全ての観測値もx以下だと言っているだけです。簡単ですね。

　さて、さらに観測値{X1,X2,...}は互いに独立であるとします。なんでそんな仮定をするかと言われると、数学的にデータを扱う場合、データが独立じゃないと、ほとんどの場合計算が全く出来ないからです。そうすると、式は以下のようになります。
$P(Z_n < x )=\Pi_{i=1}^{n}P(X_i < x)$

　さらに、この観測値がある分布関数Fに従っているとすると、以下の式になります。（分布関数＝確率密度関数を積分したもの。F(x)はデータがx以下の値になる確率）
$P(Z_n < x)=\Pi_{i=1}^{n}F(x)=F^n(x)$
　ここまで単純化されました。

　さて、nが中途半端な値というのは非常に数学的に扱いにくいです。扱いやすくする為に、n→∞と極限を取りたくなるのが人の子というものです。しかし、安易に極限をとると、F(x) < 1ならばP(Zn<x) = 0、F(x)=1ならば、P(Zn<x)=1となってしまい、非常に極端な確率（退化していると言います）が得られるだけです。この極端な確率が言うところは、無限のサンプルをとったとき、サンプルの最大値は取り得る値の最大値になるということです。なぜなら、無限にサンプルをとった場合、無限のサンプルはとりうる値全てを埋め尽くしてしまうからです。

　で、んなことは、わざわざ数式にしなくったって直感的に分かります。そんな情報はいりません。そこで、何とかサンプルの最大値Znがnの極限を取ったときに、取り得る値の最大値に収束してしまわないように、スケーリングした $Y_n$ を用意します。

$Y_n = \frac{Z_n-b_n}{a_n}$

　これが収束しなくなるような、都合の良い数列 $a_n,b_n$ が存在する理由は全く知りません。何か知らんけど、存在するそうですよ。いいんだよ、こまけぇこたぁ
　というわけで、Znの代わりにこのYnを使った $P(Y_n < x)$ は、n→∞の極限でも退化していない確率関数Gが得られるはずです。次からは $P(Y_n < x)$ についてみていきます。

　まずは、先度と同様に $P(Y_n < x)$ をFで表現してみましょう。
$P(Y_n < x) = P(\frac{Z_n-b_n}{a_n} < x)=P(Z_n < a_nx+b_n)=F^n(a_nx+b_n)$
　ZnをYnで置き換えましたが、途中式を見ると「 $x$ を $a_nx+b_n$ で置き換えた」とも取ることが出来ますね。

　さて、先も書きましたが、n→∞の極限で得られる関数Gは、何らかの退化していない関数になります。

$G(x) = \lim_{n\to\infty}F^n(a_nx+b_n)$

　で、このGって何なの？ってことですが、過去の偉大な3名の学者様により、3種類の分布関数にしかならないことが示されております。理由は知りません。

$G(x)=\exp\{-\exp(-x)\}$ .... Gumbel型
$G(x)=\exp(-x^{-\alpha})$ 　　(x≧0,α>0)...Frechet型
$G(x)=\exp\{-(-x)^\alpha\}$ 　　(x≦0,α>0)...Weibull型

　この3つを一つの式で表現すると以下のようになります。(γ=0の時にGumbel型になるのは、高校生の頃に自然対数eの定義で散々やりましたね。)

$G_\gamma(x) = \exp\{-(1+\gamma x)^{-\frac{1}{\gamma}}\}$ 　　（1+γx>0）

　いやぁ、ようやくここで求めたかった $P(Z_n < x)$ の結論が出ました。つまり、Znを適当にスケーリングすると

$\lim_{n\to\infty}P(\frac{Z_n-b_n}{a_n} < x) = G_\gamma(x)$

となる。これを求めたかったのだよ。

　ついでに、Gは分布関数ですので、微分することで確率密度関数gも得ることが出来ます。
$g_\gamma(x) = (1+\gamma x)^{-\frac{1}{\gamma}-1}\exp\{-(1+\gamma x)^{-\frac{1}{\gamma}}\}$ 　　　(γ≠0)
$g_\gamma(x) = \exp\{-x-\exp(-x)\}$ 　　　(γ=0)

　さらに、xをそのまま用いるのではなく、変数zとパラメタμ、σを用いてx=(z-μ)/σと表現し直します。σ>0です。ついでに、後の区別の為にγをζと置き直します。（γ＝ζ）
$G_\zeta(z,\mu,\sigma) = \exp\{-[1+\zeta\frac{z-\mu}{\sigma}]^{-\frac{1}{\zeta}}\}$
　こいつを一般極値分布関数といいます。一般極値分布関数の確率密度関数は
$g_\zeta(z,\mu,\sigma)=\frac{dG_\zeta(z,\mu,\sigma)}{dz}=\frac{dx}{dz}\frac{dG_\gamma(x)}{dx}=\frac{1}{\sigma}g_\gamma(x)=\frac{1}{\sigma}g_\gamma(\frac{z-\mu}{\sigma})$
となります。

　次回、極値理論を用いて推測を行ってみましょう