Generator.gumbel – ガンベル分布から乱数配列を作成する方法

2021 5/06

2020年9月14日 2021年5月6日

Generator.gumbelは、ガンベル分布から乱数配列を生成するジェネレータメソッドです。

ガンベル分布（または「最小極値分布」または「最小極値分布タイプⅠ」）は、極値分布の分散の一つで、極値問題のモデル化のために使われます。ガンベル分布は、その初期には水文学の分野で、洪水の発生、最大風速、暴風雨をモデル化するために使われていました。なぜなら、ガンベル分布は、ガウス分布（正規分布）よりも裾の部分での分布が重くなっているからです。

ガンベル分布の前は、大洪水の発生頻度はガウス分布で見積もられていたのですが、実際はその見積もりよりも高い確率で発生したので、ガンベル分布が使われるようになったのです。

なお、極値分布には他にもワイブル分布やフレシェ分布（逆ワイブル分布）があります。

このページでは、ガンベル分布の乱数を生成するGenerator.gumbelについて解説します。

1. 書式

Generator.gumbel

書き方:

Generator.gumbel(loc=0.0, scale=1.0, size=None)

パラメーター:

loc:　float or array_like of floats, optional
分布様式の位置（平均値μ）。デフォルト値は0。

scale:　float or array_like of floats, optional
分布のスケール（β）。デフォルト値は1。正の数値のみ可。

size:　int or tuple of ints, optional
出力する配列のshapeを指定する。デフォルト値のNoneの場合で、かつ、locとscaleがスカラーの場合、1つの乱数の値を返す。それ以外の場合は、np.broadcast(loc, scale).sizeの乱数を生成する。

戻り値:

out:　ndarray or scalar
パラメータを設定したガンベル分布からの乱数を出力する。

2. サンプルコード

それでは、サンプルコードを見ていきましょう。

まずは、random.default_rng コンストラクタでジェネレータオブジェクトを作成します。『numpy.random.default_rng – 乱数生成のためのジェネレータオブジェクトの作成』に目を通しておいてください。

In [1]:

import numpy as  np
rng = np.random.default_rng()
rng

Out[1]:

Generator(PCG64) at 0x7FEB1035FE58

こうして作成したジェネレータオブジェクト rng に対して、Generator.gumbelを呼び出すことによって、ガンベル分布から乱数配列を取得することができます。

第一引数locには平均値μを、第二引数scaleにはβを、第三引数には配列のsizeを指定します。

In [2]:

mu, beta = 0, 0.1
rng.gumbel(mu, beta, 10)

Out[2]:

array([ 0.38308466,  0.13275393,  0.0252544 ,  0.10355636, -0.08596615,
        0.06370004,  0.14793568,  0.02068213,  0.35166948, -0.04179268])

取得した乱数のヒストグラムを、確率密度関数とともに描いてみましょう。

In [3]:

s = rng.gumbel(mu, beta, 10000)

import matplotlib.pyplot as plt
count, bins, ignored = plt.hist(s, 30, density=True)
plt.plot(bins, (1/beta)*np.exp(-(bins - mu)/beta)
         * np.exp( -np.exp( -(bins - mu) /beta) ),
         linewidth=2, color='r')
plt.show()

以下では、ガウス分布（黄）と一緒に描いています。ガンベル分布（赤）との極値の違いを比較してみましょう。

In [4]:

means = []
maxima = []
for i in range(0,1000) :
   a = rng.normal(mu, beta, 1000)
   means.append(a.mean())
   maxima.append(a.max())
count, bins, ignored = plt.hist(maxima, 30, density=True)
beta = np.std(maxima) * np.sqrt(6) / np.pi
mu = np.mean(maxima) - 0.57721*beta
plt.plot(bins, (1/beta)*np.exp(-(bins - mu)/beta)
         * np.exp(-np.exp(-(bins - mu)/beta)),
         linewidth=2, color='r')
plt.plot(bins, 1/(beta * np.sqrt(2 * np.pi))
         * np.exp(-(bins - mu)**2 / (2 * beta**2)),
         linewidth=2, color='y')
plt.show()