逆変換法 (3) - 人工知能に関する断創録

Pythonによるモンテカルロ法入門（2014/6/20）

逆変換法 (2)のつづき。今回は、ガンマ分布とベータ分布に従う乱数を一様分布から生成してみます。

ガンマ分布

ガンマ分布のpdfとcdfは下の式。今までの分布と比べると少し複雑だな。

f:id:aidiary:20140625214018p:plain
f:id:aidiary:20140625214025p:plain
http://en.wikipedia.org/wiki/Gamma_distribution

ここで、 $\Gamma$ はガンマ関数、 $\gamma$ は不完全ガンマ関数で下の式で計算できます。

f:id:aidiary:20140625214310p:plain
f:id:aidiary:20140625214315p:plain

どちらもscipy.specialモジュール（特殊関数）で関数が定義されているのでそれを使おう。ガンマ関数は、gamma()、不完全ガンマ関数は、gammainc()で計算できます。ただ、gammaincは上の式と少し違って分母に $\Gamma(s)$ がある定義になっているので注意！
$\displaystyle gammainc(s, x) = \frac{1}{\Gamma(s)} \int_0^x t^{s-1} e^{-t} dt$
これに気付かなくて少しはまった・・・

では、逆変換法を使って一様乱数からガンマ乱数を生成してみます。前回と同じようにガンマ分布の累積分布関数の左辺をuと置いてxについて解きます。このとき、xが不完全ガンマ関数の中にあり、しかも積分の上限というとんでもない位置にあるのでどうやって取り出そう？と詰まったのですが、手抜きしてgammaincinv()を使いました。結局、xについて解くと

$x = \theta \cdot gammaincinv(k, u)$

となります。さっそく、実装。

# 逆変換法で一様分布からガンマ分布を得る
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
import scipy.special

nbins = 50

# ガンマ分布のパラメータ
k = 2.0
theta = 2.0

# 逆変換法で一様分布からガンマ分布を得る
np.random.seed()
N = 100000
U = np.random.uniform(0.0, 1.0, N)

# ガンマ分布の累積分布関数の逆関数を用いて変換
print scipy.special.gamma(k)
X1 = theta * scipy.special.gammaincinv(k, U)

# 変換したガンマ分布と理想的なpdfを描画
plt.figure(1)
rv = scipy.stats.gamma(a=k, scale=theta)
plt.hist(X1, nbins, normed=True)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((rv.ppf(0.01), rv.ppf(0.99)))

# numpyのガンマ分布に従う乱数生成関数を利用した場合
plt.figure(2)
X2 = rv.rvs(N)
plt.hist(X2, nbins, normed=True)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((rv.ppf(0.01), rv.ppf(0.99)))

plt.show()

ガンマ分布のパラメータを変えていくつか試してみます。まずは、k=2, θ=2の場合。
f:id:aidiary:20140625220144p:plain
k=1, θ=2の場合、
f:id:aidiary:20140625220206p:plain
k=9, θ=0.5の場合、
f:id:aidiary:20140625220220p:plain
生成した乱数の分布は理論値と一致していることがわかります。Wikipediaの例とも同じ結果になることが確認できました。

ベータ分布

次はベータ分布。ベータ分布のpdfとcdfは

f:id:aidiary:20140625222713p:plain
f:id:aidiary:20140625222719p:plain
http://en.wikipedia.org/wiki/Beta_distribution

ここで、 $B(\alpha,\beta)$ はベータ関数、 $B(x;\alpha,\beta)$ は不完全ベータ関数で下の式で計算できます。

f:id:aidiary:20140625222825p:plain
f:id:aidiary:20140625222828p:plain

どちらもscipy.specialモジュールで関数が定義されています。ベータ関数は、beta()、不完全ベータ関数は、betainc()で計算できます。不完全ガンマ関数と同じようにbetainc()は、分母のベータ関数を含む形で定義されています。つまり、上の累積分布関数の $I_x$ がscipyのbetainc()に当たります。

では、逆変換法を使って一様乱数からベータ乱数を生成してみます。累積分布関数をxについて解くとこちらも手抜きを使って不完全ベータ関数の逆関数betaincinv()を使うと

$x = betaincinv(\alpha, \beta, u)$

で計算できます。ではさっそく実装。

# 逆変換法で一様分布からベータ分布を得る
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
import scipy.special

nbins = 50

# ベータ分布のパラメータ
alpha = 5.0
beta = 1.0

# 逆変換法で一様分布からガンマ分布を得る
np.random.seed()
N = 100000
U = np.random.uniform(0.0, 1.0, N)

# ベータ分布の累積分布関数の逆関数を用いて変換
X1 = scipy.special.betaincinv(alpha, beta, U)

# 変換したベータ分布と理想的なpdfを描画
plt.figure(1)
rv = scipy.stats.beta(a=alpha, b=beta)
plt.hist(X1, nbins, normed=True)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((rv.ppf(0.01), rv.ppf(0.99)))

# numpyのベータ分布に従う乱数生成関数を利用した場合
plt.figure(2)
X2 = rv.rvs(N)
plt.hist(X2, nbins, normed=True)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((rv.ppf(0.01), rv.ppf(0.99)))

plt.show()

ベータ分布のパラメータを変えていくつか試してみます。まずは、α=0.5、β=0/5の場合、
f:id:aidiary:20140625225047p:plain
α=5、β=1の場合、
f:id:aidiary:20140625225120p:plain
α=1、β=3の場合、
f:id:aidiary:20140625225132p:plain
α=2、β=2の場合、
f:id:aidiary:20140625225140p:plain
α=2、β=5の場合、
f:id:aidiary:20140625225154p:plain
ベータ分布はパラメータを変えるといろんな形になるんですね。特殊関数についてはもっと深く勉強しておく必要がありそう。

逆変換法は十分試したので次回は2.2節の一般変換法に進みたいと思います。