例によって、本日もメモ代わりです。表題の通り、正規分布の最大値、最小値の「分布」について。こんなの、ちゃんと勉強した人にとっては当たり前なのかもしれませんが、素人にしてみると、「あー、そんなことくらい…えー、どうなるんだろ？」という感じの問題ですね。
＜追記＞
例によって、今回も、私以外の人がちゃんとした記事を書いていました。数式が苦手ではない方は、下記の記事を参照することをお勧めします。極値の分布に関する一般論について解説がなされています。
どうやら、計算の過程で「正規化列」なるものが出てくるそうで、そいつがわければ収束先の分布のパラメータ―もおのずとわかるようです。結局、極値というのは、元の分布の裾P(x>あるかず）という確率を積算して計算されるもののようです。だから、極値分布のパラメーターは元の確率分布の累積密度関数から計算できるのだと思います。
＜追記終わり＞
ですが、世の偉い数学者はそのくらいのこと結構昔からすでに知っており、極値理論という分野まであるそうです。彼らの成果によると正規分布の最大値、最小値の分布はガンベル分布に従う、ということが分かっています。たぶん、正規分布の畳み込みてきなことを無限回やったら収束するのがそれ、みたいな計算なのだと思いますが、数学音痴の私には理解不能でしょう。でんでんむしはともかく、絶対にょろにょろやdxと出会いそうな予感がします。ただし、正規分布の最大値、最小値がガンベル分布に収束する、ということはいたるところに書かれているし、証明っぽいものも探せばあると思います。
が、しかし、です。我々一般市民が知りたいのは、「ある」正規分布（平均と分散が既知）が与えられたとき、「その」正規分布の最大値と最小値の分布がどうなるのか？です。どうせ、計算なんてできないので、答えの方が知りたいということで…（怠け者）。というわけで調べてみたところ、下記のブログに答えがありました。
英語が読める人はこの記事を読んだらいいと、思いまーす。ただし、元の記事にも証明とかはなくて、
It was difficult to find an online reference that shows how to derive the Gumbel parameters for a normal sample of size n. I finally found a formula in the book Extreme Value Distributions (Kotz and Nadarajah, 2000, p. 9).
などと書かれているので、これって意外とマニアックなことなのかもしれません。ちなみに、ブログの著者はRick Wicklin, PhDという方になっています。PhDだぜ。ってことは、いわゆる「せんもんか」です。その、「せんもんか」でも本を見ないと分からんくらい面倒な計算なのだろうか？？なんか、素人の我々には、正規分布の確率密度の式が分かってるんだから、数学が得意な人なら、ノートの端っこでチャチャっと計算してしまいそうなイメージなんですが。ただ、本の9ページって書いてるから初歩的なことなのかも。著者の方も、一応大事を取って（？）文献に当たったのかもしれません。
ま、難しい理論は置いておいて、とりあえずコードです。コードを書けば、分かった気になりますからね。
import numpy as np
from scipy.stats import norm
from scipy.stats import gumbel_r
def norm_max_distribution(n):
mu_n = norm.ppf(1-1/n)
sigma_n = norm.ppf(1-1/n*np.exp(-1)) - mu_n
return gumbel_r(loc=mu_n, scale=sigma_n)
見ての通り、サンプル数nがなければ計算できません。考えてみれば当たり前のことですね。N(1,0)から5個引いたとき、5という数字が出てくる確率は低いですが、1000回引けば結構な確率（標準偏差が5なので、それでも一回か二回くらい？）で5が出てくるイメージです。平均0、分散（又は標準偏差）が1の正規分布の場合、その最大値（又は最小値）が従うガンベル分布のパラメータは、loc=F(1-1/n)、scale=F(1-1/n*exp(-1)) - locとなるそうです。Fというのは正規分布の累積密度関数です。要するに、正規分布のパーセンタイルからガンベル分布のパラメータが求められるのだそうです。ちなみに、パラメータをlocとかscaleと呼んでいるのは、分かりやすいのとscipyのパラメータがそういう名前になっているからです。通常、locはμやα、scaleはβなどと呼ばれることが多いとのこと。
さて、理屈上の正解は分かりましたが、私の理解が正しいとは限らないので、一応実験してみます。
# N(0,1)に対する最大値の分布を得る
g=norm_max_distribution(1000)
# N(0,1)の最大値の分布を実際に作ってみる。
vals=[]
for i in range(1000):
vals.append(np.max(np.random.normal(0., 1., 1000)))
import matplotlib.pyplot as plt
# 両方をヒストグラムにプロットする。
plt.hist(vals, bins=256)
plt.hist(g.rvs(1000), bins=256)
plt.show()
以下、プロットした結果。
おー！ちゃんとなっていそうですね！
ちなみに、言うまでもないことですが、正規分布は対称なので最大値の分布が分かれば最小値の分布もわかります。マイナスつければよいだけです。
で、大元の正規分布のパラメータが任意の場合ですが…えーと、そのときは、平均が0でない場合はその分を足せばよく、分散が1でない場合は、標準偏差に直して掛け算すればよい…と思う（自信なさげ）。というわけで、N(0,1)の最大値、最小値の分布が分かれば、一般的な正規分布の最大値、最小値の分布もわかる…はずです（自信なさげ）。やってみましょう。
全然だめじゃねーか。てか、よーく考えてみればすぐにわかりますね。だって、平均0、標準偏差1の正規分布があったとして、標準偏差が二倍になったら、とっても少ない確率で出現する極値（たとえば5とか）の大きさも二倍になると思う？正規分布ってそういう形してますか？と訊かれたら、違うと思うもの。あー、なるほどね。だから、ガンベル分布のパラメータは正規分布の累積確率密度で与えられるわけか。うすぼんやりなわたしにも、なんとなく答えが見えてきたぜ。こーゆーことを「やってみないと」分からないのが、数学音痴のダメなところです、というわけで、若人の皆さんはちゃんと数学勉強しようね。算数じゃなくて、数学ね。
では、この知見（？）に従ってコードを治します。
def norm_max_distribution(n,m,s):
mu_n = norm(m,s).ppf(1-1/n)
sigma_n = norm(m, s).ppf(1-1/n*np.exp(-1)) - mu_n
return gumbel_r(loc=mu_n, scale=sigma_n)
上記の関数を使って、平均10で、標準偏差13の正規分布について、最大値の分布を計算してみます。
あってそうですね。よかったよかった。最小値の方も、似たようなもので下記のように計算できます（たぶん）。
def norm_min_distribution(n,m,s):
mu_n = norm(m,s).ppf(1/n)
sigma_n = norm(m, s).ppf(1-1/n*np.exp(-1)) - norm(m,s).ppf(1-1/n)
return gumbel_l(loc=mu_n, scale=sigma_n)
パラメータlocに関しては、最小値なので正規分布のテールの反対側からとってくる。ガンベル分布の広がり方は分散に依るので、最大値の分布が欲しいか、最小値の分布が欲しいかには依らないはず、そして分布は丁度ひっくり返ったものになるはず、という、うすぼんやりの理解で書いたものです。まぁ、とりあえず、最大／最小の分布が見てみたかっただけなので、こんなもんと言うことで…。こっちの関数は、私がでっち上げたもんなので、ちょっと信用がおけないのですが、多分あってると思う（実験しました！いや、分かってるけど！計算できないから実験するしかないのだ！）。
ただ、ガンベル分布って名前は難しそうですが、確率密度関数の形は簡単でgumbel_rがexp(-(x + exp(-x)))、gumbel_lがexp(x - exp(x)))だそうです。いや、簡単そうに見えるだけで簡単ではないか。他のscipyの確率分布と同じで、パラメータの意味はy = (x - loc) / scaleだそうなので、多分、考え方はあってると思う…やっぱり、こういうとき、数学ができるとちょちょっと計算して答えを出すということができるので、勉強は真面目にしておくべきですね。
「正規分布」なんて超基本なのに、じゃあ、その最大値の分布は？と訊かれると意外とわからなかったので、計算方法くらいは知っておきたいと思ってやってみた次第です。
＜追記＞
極値統計について、詳しく説明されたありがたいスライドがありました。そこに、最大値の分布から最小値の分布を求める方法があったので抜粋します。
参照元：https://ibisml.org/ibis2009/pdf-invited/takahashi1.pdf
というわけで、極値の分布をG(x)とすると、1-G(-x)で求まるということです。ちなみに、このスライドには、いくつかの確率分布に対する極値の分布の導出も載っています。ふわっとした理解ですが、計算の手順としては：
- 確率変数x(n)の数列を考えて、x(n) = a(n)*x(n) + b(n)という確率変数の数列（？）を作る。
- このときP(max(X(i) <= x(i)))なるものを考えると、これが最大値X(i)が出現する確率になる。
- X(i)がなんとかより大きくなる確率というのは累積確率密度関数F(x)を考えればよい。
- 後は、そいつを掛け算すればよいのでF(x(n))^nとおき、n→∞の極限をとる。
- その収束先Gが最大値の分布。
というロジックのようです。たぶん難しいのは、ちゃんと収束するようなa(n)やb(n)を見つけるところで、やっぱり素人の我々は、ノートの切れ端で計算するより、答えを探した方が早そうです（なげやり）。