サンプリングアルゴリズムについていくつか紹介します. その中で理解を深めるために Python で簡単な実装と可視化を行います.
第一回目は 逆関数法 についてです. 逆関数法では [0,1] 区間の一様分布(=$U(0,1)$)から得られた乱数を変換することで任意の確率分布に従う乱数を得る方法です.
逆関数法は次のような問題に答えを与えてくれます
問題
確率密度関数 $f(x)$ とその累積分布関数の逆関数 $F(x)^{-1}$ の関数形がわかっているが、 python のライブラリに $f(x)$ に従う乱数を出力するものが実装されていないとき、 どのように $f(x)$ に従う乱数を生成するのが良いでしょう?
サイコロを作ろう(モンテカルロ法)
一様乱数を用いて出目の出現確率が等価なサイコロを作るにはどうすればよいでしょう？答えは
- $0\leq u <\frac{1}{6}$ のとき1
- $\frac{1}{6}\leq u <\frac{2}{6}$ のとき2
- ・・・
- $\frac{5}{6}\leq u <\frac{6}{6}$ のとき6
のような振り分けを行えば作れます.
実際にPythonで実装してみましょう.
# 必要なライブラリのインポート
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
パラメータの設定
# パラメータの設定
np.random.seed(0)
N = 100000
# 確率分布関数を rv に格納
# [0.0, 1.0] の一様分布からの確率変数
rv = uniform(loc = 0.0, scale = 1.0)
指定した区間に入ると1を返す関数を定義します. これがサイコロを作るうえで重要な関数です.
def section(u,a,b):
if (a<=u)&(u<b):
return 1
else:
return 0
乱数を生成してどの目が出たのかを 0 or 1 で管理
dice = [[section(u,num/6,(num+1)/6) for u in rv.rvs(size = N)]
for num in range(0,6)]
# 出目の出現頻度を確認
np.array([np.sum(dice[num]) for num in range(0,6) ] )/ N
# np.array([0.166608, 0.167029, 0.166872, 0.166501, 0.1672 , 0.166226])
だいたい6等分の確率になっており、サイコロが実装できました.
ここで重要なのは [0,1] 区間を適切に区切ることで任意の確率分布を作成することができることです. これは, 区間の幅が大きければ大きいほど実現確率が高くなり, 小さいものは実現確率が小さくなるように設定することです. これが確率の本質的な性質で, ある事象の発生確率は全体のうちに占める面積の割合で決まることに由来しています.
例えば $0\leq u <\frac{1}{2}$ のとき1が出るようにして残りの区間を2~6の出目が出るように均等に割り振ることで1の目が出やすいインチキサイコロも作成できます.
このように一様乱数の値を適切な区間に区切ることで目的の確率分布を作成できる事がわかる, これがモンテカルロ法です.
逆関数法
それでは逆関数法について解説します.
まず、確率密度関数の復習です. まず
- $f(x)$は確率密度関数
- $F(x)$は累積分布関数
とします. 確率変数が連続値を取る場合, $x$ が出現する確率はある微小な値 $\delta $ に対して
p(x) = F(x+\delta) - F(x) = \int_{x}^{x+\delta} dx \,\, f(x) \sim f(x)\delta
で与えられます. この式を見れば分かる通り, 確率の本質は区間 [x,x+\delta] における面積(=積分)で表されることです.
累積分布関数の値域は$[0,1]$であることに注意すると, 逆関数 $F^{-1}(u)$ に対して $[0,1]$区間の一様分布の乱数 $u \sim U(0,1)$を引数とすることができます.
モンテカルロ法の重要な点は欲しい分布関数を得るためには, それに合致するように適切な区間を区切って一様分布の値を変換することでした. ここでは分布の関数形がわかっている状態において, 出力される値に対して$x\leq F^{-1}(u) < x+\delta$ のとき $x$ を観測したとみなします. これが逆関数法でありうまく確率の区分を区切ることに相当します. つまり, 先程のサイコロを実装した場合と同じことを行っています. 例えば $u \in [0,1/6)$ のとき 1 の目が出たと解釈するのと同じです.
逆関数法の手順をまとめると:
- 大量の一様乱数 $u$ を生成
- $F^{-1}(u)$を計算
- 値が区間 $[x, x+\delta ]$ に入っている場合$x$が発生したとして扱う
- 区間の観測された件数を総数で割ることで $p(x)$ が得られる
です.
(注意)
累積分布関数$F(x)$は単調増加関数なので必ず逆関数が存在します.
例1：指数分布
具体例を通じて確認しましょう.
パラメータ $\lambda$ の指数分布
f(x|\lambda) = \lambda e^{-\lambda x}
について, その累積分布関数は
F(x| \lambda) = 1 - e^{-\lambda x}
です.
その逆関数を求めるにはまず,
F(x) = u
として $x$ について解き直すと
x = F^{-1}(u) = -\frac{1}{\lambda}\log(1-u)\,,
として得られます.
$u \sim U(0,1)$ を $F^{-1}(u)$に代入した値が指数分布からのサンプリングの結果として解釈できます.
Python で実装
先程の内容を Python で実装して指数分布に従う乱数を逆関数法で生成してみます.
まず必要なライブラリのインポート:
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats
# シードの固定
np.random.seed(0)
パラメータの設定:
nbins = 50
# 指数分布のパラメータ(scale = 1/lambda)
scale = 1.0
N = 100000
一様分布の生成してプロット:
U = scipy.stats.uniform(loc = 0.0, scale = 1.0).rvs(size = N)
# 生成した一様乱数を描画
plt.figure()
plt.hist(U, nbins, density = True)
plt.show()
次のような画像が生成されるはずです.
[0,1]区間のどの値も同じ確率で発生している事がわかります.
次に, 先程生成した一様分布に従う乱数 $u$ を累積分布関数の逆関数 $F^{-1}(u) = -\frac{1}{\lambda}\log(1-u)$ に代入して変換します.
# 指数分布の累積分布関数の逆関数を用いて変換
X1 = - scale * np.log(1-U)
変換された値をヒストグラムにして描画:
# 変換した指数分布の乱数と理想的な確率密度関数(PDF)を描画
plt.figure()
# 変換した分布関数のプロット
plt.hist(X1, int(nbins * 10 /5), density = True)
# 指数分布の確率密度関数を定義
rv = scipy.stats.expon(scale = scale)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
# 真の確率密度関数(pdf)をプロット
plt.plot(x, rv.pdf(x), 'r-', lw = 2)
# 上端と下端を設定
plt.xlim(rv.ppf(0.01),rv.ppf(0.99))
plt.show()
この画像を見れば分かるように逆関数法によって生成された乱数が真の確率密度関数に従っている事がわかります.
逆関数法の理解
ここで先程のコードを使いながら逆関数法がなぜ有効な乱数精製方法であるかを深堀りしてみましょう.
$u = 0.0001$ に対して$\delta = 0.0001 $ の幅をもった二点$u, u + \delta$を累積分布関数の逆関数に入れると
- scale * np.log(1-np.array([0.0001,0.0002,]))
# np.array([0.00010001, 0.00020002])
となり、出力結果の幅はおおよそ $0.0001$ 程度になります.
一方で幅は $\delta = 0.0001 $ で同じだが, $u$の値を0.9990に変えてみましょう.
- scale * np.log(1-np.array([0.9990,0.9991,]))
# np.array([6.90775528, 7.01311579])
となり、出力結果の幅はおおよそ $0.1$ 程度になる.
つまり, 幅 $\delta$ が一定でも一様分布に従う乱数の値 $u$ によって $F^{-1}(u)$ の出力結果の幅$[F^{-1}(u),F^{-1}(u)+\delta]$が異なっており, この差がヒストグラムにあらわれてきているのです. つまり, 指数分布の場合
- $u$の値が0に近い→$F_{X}^{-1}(u)$ の値の幅が小さくなる
- $u$の値が1に近い→$F_{X}^{-1}(u)$ の値の幅が大きくなる
ということです.
「$u$ の値が0に近いほど$F_{X}^{-1}(u)$ の値の幅が小さくなる」
ということは
「$F_{X}^{-1}(u)$ の区間を固定するとその区間に入る $u$ の件数が多くなる」
ということであり、まさに指数分布の性質を現しています.
下のグラフは $F^{-1}(u)$ をプロットしたものです.
青色と緑色の矢印で示された区間を比較してみると, 縦軸の $F_{X}^{-1}(u)$ の区間を固定していますが, $F_{X}^{-1}(u)$ の値が大きくなる (つまり, $f(x)$の値が大きくなる) ほど 横軸の $u$ の区間が小さくなっています. これはまさに指数分布の性質を反映しています.
横軸の $u$ の値が小さいものほど縦軸の世界ではほとんど同じ値に丸め込まれておりこれはつまり同じ値が実現していると解釈できます.
これが逆関数法の本質であり、サイコロの場合と同様に,実現確率が高い値の区間ほど発生件数が高いと判定されるように, 一様分布を適切な区間で区切ることがてきました.
例2：ロジスティック分布
もう一つ同様に具体例でチェックしてみましょう.
ロジスティクス分布の確率密度関数、累積分布関数、累積分布関数の逆関数は
\begin{align}
f(x|\mu,s) &= \frac{1}{s(1+e^{-\frac{x-\mu}{\mu}})}e^{-\frac{x-\mu}{\mu}} \\
F(x|\mu,s) &= \frac{1}{1+e^{-\frac{x-\mu}{\mu}}}\\
F^{-1}(u|\mu,s) &= \mu + s \ln(\frac{u}{1-u})
\end{align}
です.
先ほどと同様に逆関数法を実装してみましょう
パラメータの設定:
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))
得られた結果を描画
plt.figure()
# 逆関数法で得られた確率密度関数をプロット
plt.hist(X1, nbins, density = True)
# 真の確率密度関数をプロット
rv = scipy.stats.logistic(loc = mu, scale = s)
# 真の確率密度関数
x = np.linspace(rv.ppf(0.01),rv.ppf(0.99),1000)
plt.plot(x,rv.pdf(x), 'r-', lw = 2)
plt.xlim(rv.ppf(0.01),rv.ppf(0.99))
true_pdf(rv)
plt.show()
これが逆関数法で生成されたヒストグラムです.
一方で numpy で出力されるのロジスティック分布に従う乱数生成関数を利用してみましょう.
plt.figure()
X2 = rv.rvs(N)
plt.hist(X2,nbins, density = True)
# 真の確率密度関数をプロット
true_pdf(rv)
plt.show()
ほとんど同じヒストグラムが生成されているはずです。