書記が数学やるだけ#785 乱数による確率分布の生成
マルコフ連鎖モンテカルロ法の理解を目標に，サンプリングについていくつか取り上げていく。
問題
説明
コンピュータは原則確定的な計算しか行うことはできず，完全な乱数を作ることはできない。その代わりとして，できるかぎり規則性をなくした擬似乱数が用いられる。その始祖は1946年頃ノイマンが提案した平方採中法で，以後さまざまな手法が開発されてきた。
上の手法で得られた一様乱数から，任意の確率分布を得る方法をいくつか挙げる。
逆関数法とは，累積分布関数の逆関数を用いて，標準一様分布に従う確率変数から所望の分布に従う確率変数を生成させる方法である。これは逆関数が陽に求められる分布に有効である。
変数変換法は，ある分布から別の分布を生成する方法で，特にBox-Muller変換が有名である：
棄却法は，目標分布にしたがう乱数を生成するために，代わりに提案分布の乱数を生成してある条件を満たす場合だけ生成した乱数を目標分布の乱数とみなす方法である。
解答
逆変換法の原理は至ってシンプルである。
逆変換法は離散確率分布にも適用できる。
Pythonで逆変換法を実装してみる。
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats
import scipy.optimize
from scipy.stats import uniform, norm, beta一様乱数は以下の通り。
#一様乱数の取得
np.random.seed(314)
N = 100000
U = stats.uniform(loc=0, scale=1).rvs(size=N)
nbins = 50
plt.hist(U, nbins, density=True)
一様乱数に対し，上に示した変換式を適用することで指数乱数が得られる。
#指数分布のパラメータ（scale = 1/lambda）
scale = 1.0
#指数分布の累積分布関数の逆関数を用いて変換
X1 = - scale * np.log(U)
#実データ：青 赤線：理論値
rv = stats.expon(scale=scale)
plt.hist(X1, nbins, density=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()
次にBox-Muller変換を実装する。これも一様乱数に変換式を適用するだけ。
#一様乱数
np.random.seed(314)
N = 100000
rv1 = uniform(loc=0, scale=1)
rv2 = uniform(loc=0, scale=1)
U1 = rv1.rvs(N)
U2 = rv2.rvs(N)
#Box-Muller変換
X1 = np.sqrt(-2 * np.log(U1)) * np.cos(2 * np.pi * U2)
X2 = np.sqrt(-2 * np.log(U1)) * np.sin(2 * np.pi * U2)
#理論値(赤線)
rv = norm(loc=0, scale=1)
#X1の分布
plt.subplot(2, 1, 1)
nbins = 50
plt.hist(X1, nbins, density=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)))
#X2の分布
plt.subplot(2, 1, 2)
plt.hist(X2, nbins, density=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()
最後に棄却法について，まずは理論的背景を確認しておく。
ベータ分布の例を以下に示す，ここで分布の最大値をRとしている。
実装において，まずRを決めるために最大値問題を解く。
#目標分布P
p = beta(a=2, b=6).pdf
#提案分布Q（一様乱数）
qv = uniform
q = qv.pdf
#分布の上限R
xopt = scipy.optimize.fmin(lambda x: -p(x), 0, disp=False)
R = p(xopt)[0]
print("R =", R)
#出力：R = 2.8131429777437695提案分布からの乱数Yと一様乱数Uを生成し，棄却アルゴリズムを適用する。
Nsim = 100000
#提案分布Qからの乱数Yを生成
Y = qv.rvs(size=Nsim)
#一様乱数UをNsim個生成
U = uniform.rvs(size=Nsim)
#棄却アルゴリズム
X = Y[U <= p(Y) / (M * q(Y))]
print("サンプル数: %d => %d" % (len(Y), len(X)))
print("実際の受理率 : %f" % (len(X) / float(len(Y))))
print("理論的な受理率: %f" % (1.0 / M))
#出力
#サンプル数: 100000 => 35547
#実際の受理率 : 0.355470
#理論的な受理率: 0.355474結果を以下に可視化する。
#目標分布（赤線）
x = np.linspace(0.0, 1.0, 1000)
y = p(x)
plt.plot(x, y, 'r-', lw=2)
#提案分布（緑線）
y = M * uniform.pdf(x)
plt.plot(x, y, 'g-', lw=2)
#実測値
plt.hist(X, bins=50, density=True)
plt.show()
本記事のもくじはこちら：
リニューアル
