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

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

逆変換法（2014/6/22）のつづき。今回は、指数分布ではなく、ロジスティック分布とコーシー分布に従う乱数を一様乱数から変換して生成してみます。Rによるモンテカルロ法入門の練習問題2.2のPython実装です。

ロジスティック分布

ロジスティック分布の確率密度関数（pdf）と累積分布関数（cdf）は下の式になります。

f:id:aidiary:20140623213002p:plain
f:id:aidiary:20140623213006p:plain
引用: http://en.wikipedia.org/wiki/Logistic_distribution

ロジスティック分布は、 $\mu$ と $s$ の2つのパラメータを取ることがわかります。scipy.stats.logisticではlocとscaleという名前になっています。ロジスティック分布のcdfがニューラルネットでおなじみのロジスティック・シグモイド関数なのか。初めて知った。

では、逆変換法を使って一様乱数からロジスティック分布に従う乱数を生成してみます。前回と同じようにロジスティック分布の累積分布関数の左辺をuと置いてxについて解くと

$\displaystyle x = \mu + s \log \frac{u}{1-u}$

となります。この逆関数を使って一様乱数を変換します。

# 逆変換法で一様分布からロジスティック分布を得る
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

nbins = 50

# ロジスティック分布のパラメータ
mu = 0
s = 1

# 逆変換法で一様分布からロジスティック分布を得る
np.random.seed()
N = 100000
U = np.random.uniform(0.0, 1.0, N)

# ロジスティック分布の累積分布関数の逆関数を用いて変換
X1 = mu + s * np.log(U / (1 - U))

# 変換したロジスティック分布と理想的なpdfを描画
plt.figure(1)
rv = scipy.stats.logistic(loc=mu, scale=s)
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()

実行すると

f:id:aidiary:20140623213917p:plain

となります。真のロジスティック分布の確率密度関数とほぼ一致していることも確認できました。

コーシー分布

次はコーシー分布です。正規分布と似てますが、より裾野が広い分布。コーシー分布のpdfとcdfは

f:id:aidiary:20140623222349p:plain
f:id:aidiary:20140623222355p:plain
http://en.wikipedia.org/wiki/Cauchy_distribution

逆変換法を使って一様乱数からコーシー分布に従う乱数を生成してみます。コーシー分布の累積分布関数の左辺をuと置いてxについて解くと

$\displaystyle x = x_0 + \gamma \tan (\pi (u - \frac{1}{2}))$

となりました。scipy.stats.cauchyでは、 $x_0$ がloc、 $\gamma$ がscaleに当たります。では、実装！

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

nbins = 50

# コーシー分布のパラメータ
x0 = 0
gamma = 1

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

# コーシー分布の累積分布関数の逆関数を用いて変換
# コーシー分布は裾野が広く極端な値も出やすいためplotするときは
# それらのデータを切り捨てている
X1 = x0 + gamma * np.tan(np.pi * (U - 0.5))
X1 = X1[(X1>-10) & (X1<10)]

# 変換したコーシー分布と理想的なPDFを描画
plt.figure(1)
rv = scipy.stats.cauchy(loc=x0, scale=gamma)
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((-10, 10))

# numpyのコーシー分布に従う乱数生成関数を利用した場合
plt.figure(2)
X2 = rv.rvs(N)
X2 = X2[(X2>-10) & (X2<10)]
plt.hist(X2, nbins, normed=True)
x = np.linspace(-10, 10, 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((-10, 10))

plt.show()

コーシー分布を描画するときは要注意。確率密度関数の裾野が広いので極端に小さい、大きい値もわりと多く生成されます。実際にX1の値でmin()とmax()を取ってみるとそのような値も生成されていることがわかります。そのままヒストグラムを描画すると真っ平らになってしまうため描画するときは極端に大きな乱数は削除しています。

実行結果は、下のようになり理論値とも一致しています。

f:id:aidiary:20140623223008p:plain

参考

http://yagays.github.io/blog/2012/10/21/imcmr-2-1/