逆変換法

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

逆変換法（Inverse transform sampling）は、[0,1]区間の一様分布から得られた乱数（一様乱数）を変換することで任意の確率分布に従う乱数を得る手法とのこと。このとき必要となるのは変換先の確率分布の累積分布関数（cumulative distribution function: cdf）。累積分布関数（cdf）は確率密度関数（pdf）に比べると影が薄い気がしていましたがここでは大活躍です。逆変換法は、累積分布関数 $F_X$ の逆関数 $F_X^{-1}$ を使って一様乱数を変換します。下の例は、一様乱数を指数分布の乱数に変換した例になります。

なぜこれでよいのかわかりづらいのですが、@teramonagiさんの下記の資料が直感的に説明してくれています。

Rで学ぶ逆変換（逆関数）法 from tera monagi

指数分布の例

指数分布（exponential distribution）の確率密度関数と累積分布関数は

f:id:aidiary:20140622110156p:plain
f:id:aidiary:20140622110202p:plain
引用: http://en.wikipedia.org/wiki/Exponential_distribution

累積分布関数の逆関数を求めます。左辺をuと置いてxについて解くと

$\displaystyle x = - \frac{1}{\lambda} \log (1 - u)$

となりました。さっそく逆変換法で一様乱数から指数分布の乱数が生成できるか確かめてみます。

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

nbins = 50

# 指数分布のパラメータ（scale = 1/lambda）
scale = 1.0

# 逆変換法で一様乱数から指数分布の乱数を得る
np.random.seed()
N = 100000
U = scipy.stats.uniform(loc=0.0, scale=1.0).rvs(size=N)

# 指数分布の累積分布関数の逆関数を用いて変換
X1 = - scale * np.log(1 - U)

# 生成した元の一様乱数を描画
plt.figure(1)
plt.hist(U, nbins, normed=True)

# 変換した指数分布の乱数と理想的なPDFを描画
plt.figure(2)
rv = scipy.stats.expon(scale=scale)
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)))

plt.show()

実行すると下のようなグラフが得られました。変換した乱数の分布を見ると確かに指数分布に従うようです。赤色の線が本当の指数分布の確率密度関数です。

f:id:aidiary:20140622110814p:plain
f:id:aidiary:20140622110821p:plain

ちなみにですが、指数分布に従う乱数が必要なときはわざわざ一様乱数から自分で逆変換しなくてもよくてscipyにscipy.stats.exponという関数が用意されています。この実装に何を使っているのか気になる・・・逆変換法なのかな？

# scipyの乱数生成関数を利用した場合
plt.figure(3)
rv = scipy.stats.expon(scale=scale)
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:20140622111220p:plain

指数分布から一様分布に変換

逆変換法の逆（順変換？）ですが指数分布の乱数を累積分布関数で変換すると本当に一様乱数になるのかも実験で確かめてみました。累積分布関数はその逆関数とちがってscipy.stats.expon.cdfという関数がすでに用意されているのでそのまま使います。

# 指数分布においてU=F(X)を確認
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

# 指数分布に従うサンプルを生成
N = 100000
nbins = 50

figure(1)
X = scipy.stats.expon(scale=1).rvs(size=N)
plt.hist(X, nbins, normed=True)

# 指数分布の累積分布関数（cdf）で変換
figure(2)
U = scipy.stats.expon.cdf(X)
plt.hist(U, nbins, normed=True)

plt.show()