一般変換法 - 人工知能に関する断創録

Pythonによるモンテカルロ法入門（2014/6/20）といいながらいまだモンテカルロ法が出てきていないという・・・

今日は一般変換法を実験。Rによるモンテカルロ法入門2.2節では、指数分布に特定の変換を施すことで生成できる標準的な分布としてカイ二乗分布、ガンマ分布、ベータ分布が紹介されています。以下の数式で $X_j$ はすべてパラメータ $\lambda=1$ の指数分布から生成された乱数です。

指数乱数をカイ二乗乱数に変換

$\displaystyle Y = 2 \sum_{j=1}^\nu X_j \sim \chi_{2 \nu}^2$

まだ理屈はよくわからないけれどこの式でできるみたい。たとえば、自由度6のカイ二乗分布 $\chi_6^2$ に従う乱数が欲しい場合は、 $\nu=3$ なので3つの指数分布から生成された乱数を足し合わせることで1つのカイ二乗分布に従う乱数ができます。さっそく検証してみよう。

指数乱数はscipy.stats.expon.rvs()で直接生成できますが、せっかくなので復習をかねて一様乱数から逆変換法で生成（6/22）してみました。つまり、一様乱数→指数乱数→カイ二乗乱数の順に変換しています。

# カイ二乗分布に従う乱数を指数乱数から生成
import numpy as np
from scipy.stats import uniform, chi2
import matplotlib.pyplot as plt

nbins = 50

# カイ二乗分布のパラメータ（自由度: degrees of freedom）
# この方法では偶数のみ
df = 6
if df % 2 != 0:
    print u"ERROR: 自由度は偶数のみ"
    exit(0)

nu = int(df / 2)

# 逆変換法で一様乱数を指数乱数に変換
# 1つのカイ二乗乱数を得るためにはnu個の指数乱数が必要なので
# 行列化して和を取りやすくする
np.random.seed()
N = 100000
rv = uniform(loc=0.0, scale=1.0)
U = rv.rvs(nu * N).reshape((nu, -1))
X = - np.log(1 - U)

# 指数乱数からカイ二乗乱数を得る
# 各列を足しあわせて1つのカイ二乗乱数を得る
Y = 2 * np.sum(X, axis=0)

# 変換したカイ二乗分布の乱数とpdfを描画
plt.figure(1)
plt.hist(Y, nbins, normed=True)
rv = chi2(df=df)
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()

f:id:aidiary:20140626213453p:plain

変換した乱数の分布（ヒストグラム）は理論値と一致していることが確認できます。これはすごい！

指数乱数をガンマ乱数に変換

$\displaystyle Y = \theta \sum_{j=1}^k X_j \sim G(k, \theta)$

ガンマ分布の第1パラメータk個分の指数乱数を足し合わせるとガンマ乱数が1つ生成できるみたい。これも理屈がわからないけどさっそく検証してみよう。

# ガンマ分布に従う乱数を指数乱数から生成
import numpy as np
from scipy.stats import uniform, gamma
import matplotlib.pyplot as plt

nbins = 50

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

# 逆変換法で一様乱数を指数乱数に変換
# 1つのガンマ乱数を得るためにはk個の指数乱数が必要なので
# 行列化して和を取りやすくする
np.random.seed()
N = 100000
rv = uniform(loc=0.0, scale=1.0)
U = rv.rvs(k * N).reshape((k, -1))
X = - np.log(1 - U)

# 指数乱数からガンマ乱数を得る
# 各列を足しあわせて1つのガンマ乱数を得る
Y = theta * np.sum(X, axis=0)

# 変換したガンマ分布の乱数とpdfを描画
plt.figure(1)
plt.hist(Y, nbins, normed=True)
rv = gamma(a=k, scale=theta)
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()

f:id:aidiary:20140626213502p:plain

指数乱数をベータ乱数に変換

$\displaystyle Y = \frac{\sum_{j=1}^\alpha X_j}{\sum_{j=1}^{\alpha+\beta} X_j} \sim Be(\alpha, \beta)$

これは少し複雑。さっそくプログラムを書いて検証。

# ベータ分布に従う乱数を指数乱数から生成
import numpy as np
from scipy.stats import uniform, beta
import matplotlib.pyplot as plt

nbins = 50

# ベータ分布のパラメータ
a = 5.0
b = 1.0

# 逆変換法で一様乱数を指数乱数に変換
# 1つのガンマ乱数を得るためにはk個の指数乱数が必要なので
# 行列化して和を取りやすくする
np.random.seed()
N = 100000
rv = uniform(loc=0.0, scale=1.0)
U = rv.rvs((a+b) * N).reshape((a+b, -1))
X = - np.log(1 - U)

# 指数乱数からベータ乱数を得る
Y = np.sum(X[0:a,], axis=0) / np.sum(X, axis=0)

# 変換したベータ分布の乱数とpdfを描画
plt.figure(1)
plt.hist(Y, nbins, normed=True)
rv = beta(a=a, b=b)
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()

f:id:aidiary:20140626213508p:plain