はじめに

ネットで、ギブスサンプリングについて調べていた時に、「そういえばサンプリングってそんなにしっかり勉強したことないな」と思って、何かいい本ないか探してたら、たまたま手元にPRMLがあった（pdfは持ち運びやすくていいですね）のでやって見ることにしました。
どうせなら、pythonでプログラム書いたり、本の内容を自分なりに説明したらより理解できるし、どこかで必要としている人の役に立つかな？って思ったので投稿してみました。
周りに理論をしっかりやってる人がいなく、独学ですので理解が不十分なところが多々あると思います。
ですので、ご指摘等していただけると嬉しいです！

続きの章・節に関しましても挫折しないところまではやろうと思っています（笑）。

11.1.1 標準的な分布

zを区間(0,1)の一様分布従う変数とする。
ここで、zを変換して、ある分布p(y)に従う変数yを得たいとする。
zをyに変換する関数を$f(\cdot)$とすれば、
$$y = f(z)$$
と表せる。

では、変数yがp(y)に従うようにするためには、$f(\cdot)$をどのように選べば良いだろうか。
ここで、$p(y)$は、

p(y) = p(z) \left |\frac{dz}{dy} \right|

に従う。（今、zは一様分布に従うのでp(z)=1）
そこで、上記の式を積分すると、

\begin{align}
\int p(y)dy &= \int \left |\frac{dz}{dy} \right|dy \\
&= \int dz \\
&= z
\end{align}

よって、$p(y)$の不定積分を$h(y)$とおき、不定積分の表記の仕方をちょっと変えると、

\begin{align}
z = h(y) = \int_{-\infty}^{y} p(y)dy
\end{align}

と表せる。
ここで、逆関数をとると、

y = h^{-1}(z)

となる。
よって、変数yがp(y)に従うようにするためには、zを変換する関数$f(\cdot)$を以下のように定めれば良い。

f(z) = h^{-1}(z)

このままだと良く分からない（自分は良く分からなかった）ので、例題を見ていきます。

例1: 指数分布

例えば、p(y)を指数分布としよう。

p(y) = \lambda exp(-\lambda y)

つまり、「区間(0,1)の一様分布に従う変数zを指数分布p(y)に従うような変数yに変換したい」ということである。

zをyに変換する関数$f(z)$は以下のように表せた。

y = f(z) = h^{-1}(z)

ここで、p(y)の定義域は、$(0 \leq y \leq \infty)$ であるので、

\begin{align}
z = h(y) &= \int_{-\infty}^{y} p(y)dy \\
&= \int_{0}^{y} \lambda exp(-\lambda y)dy \\
&= 1 - exp(-\lambda y)
\end{align}

よって、zをyに変換する以下の式が得られる。

y = h^{-1}(z) = -\frac{1}{\lambda}log(1-z)

この変換から得られたyは、指数分布に従うらしい。

というわけで、実装して確認。

コード

z = np.random.uniform(0, 1, 1000)

# 指数分布のλ
lam = 1.5

# 指数分布の式
exp_dist = lambda x: lam * np.exp(-lam * x) 

y = -np.log(1 - z) / lam

# yのヒストグラム
plt.hist(y, bins=50, range=(0, 4), color='aqua', label='サンプリング')

# 指数分布のグラフ
x = np.linspace(0, 4, 1000)
plt.plot(x, exp_dist(x)*80, 'darkblue', label='指数分布')

plt.xlabel('y')
plt.legend(loc='best', prop={'size': 12})

ちゃんと指数分布に従っていますね！

例2: 標準コーシー分布

変換法が適用できるもう１つの例として、コーシー分布を挙げる。

コーシー分布は以下の式で表される。

p(x;x_0, \gamma) = \frac{1}{\pi} \frac{\gamma}{(x-x_0)^2 + \gamma}

もうちょっと細かいことは、ググるなり本見るなりして調べてください！

そして、$x_0 = 0, \gamma = 1$であるコーシー分布を標準コーシー分布と言って、以下の式で表す。

p(x;0, 1) = \frac{1}{\pi} \frac{1}{x^2 + 1}

この例では、変数zを標準コーシー分布に従う変数yに変換するための関数を求めたいと思います。

まず、yについての不定積分$h(y)$を求めたいと思います。
標準コーシー分布$p(y)$の定義域は、$-\infty \leq y \leq \infty$なので、

\begin{align}
z = h(y) = \int_{-\infty}^{y}p(y)dy &= \int_{-\infty}^{y} \frac{1}{\pi} \frac{1}{y^2 + 1} dy \\
&=  \frac{1}{\pi} \int_{\frac{\pi}{4}}^{tan^{-1}y} 
\frac{1}{1 + tan^2\theta} (1 + tan^2\theta) d\theta \\
&=  \frac{1}{\pi} \int_{\frac{\pi}{4}}^{tan^{-1}y} d\theta\\
&= \frac{1}{\pi}tan^{-1}y + \frac{\pi}{4}
\end{align}

となる。
(上記の式では、$y = tan\theta$ と置換して積分しました。)

よって、zからyへの変換関数は以下のようになる。

y = h^{-1}(z) = tan\left(z - \frac{1}{4}\right)\pi

コード

# 標準コーシー分布の式
cauchy_dist = lambda x: 1 / (np.pi * (1 + x**2)) 

z = np.random.uniform(0, 1, 1000)
y = np.tan(np.pi * (z - 0.25))

# yのヒストグラム
plt.hist(y, bins=50, range=(-10, 10), color='aqua', label='サンプリング')


# 標準コーシー分布のグラフ
x = np.linspace(-10, 10, 1000)
plt.plot(x, cauchy_dist(x)*360, 'darkblue', label='標準コーシー分布')

plt.xlabel('y')
plt.legend(loc='best', prop={'size': 10})

ちゃんと標準コーシー分布に従ってるのがわかりますね。

Box-Muller法

一様分布に従う変数$z_1, z_2$から、標準ガウス分布(平均０、分散１のガウス分布）に従う変数$y_1, y_2$へ変換する方法

この時のyは、以下で定義される。

y_1 = \sqrt{-2log z_1}cos{2\pi z_2} \\
y_2 = \sqrt{-2log z_1}sin{2\pi z_2}

では、このように定義された$y_1, y_2$が本当に標準ガウス分布に従うか調べて見る。

まず、$y_1, y_2$の同時分布は以下のように定義される。

p(y_1, y_2) = p(z_1, z_2) \left| \frac{\partial(z_1, z_2)}{\partial (y_1, y_2)} \right| \\

ここで、$ \frac{\partial(z_1, z_2)}{\partial (y_1, y_2)} $は、ヤコビアン（ヤコビ行列の行列式）で、以下のように定義される。

\frac{\partial(z_1, z_2)}{\partial (y_1, y_2)} = \left|
\begin{array}{cc}
\frac{\partial z_1}{\partial y_1} & \frac{\partial z_2}{\partial y_1} \\
\frac{\partial z_1}{\partial y_2} & \frac{\partial z_2}{\partial y_2} \\
\end{array}
\right|

従って、ヤコビ行列の各編微分を求めて行列式を計算すればよい。

\begin{align}
y_1^2 + y_2^2 &= -2log{z_1} \\
\therefore z_1 &= exp\left( -\frac{y_1^2 + y_2^2}{2} \right) \\
\\
\frac{y_2}{y_1} &= tan\left(2\pi z_2\right) \\
\therefore z_2 &= \frac{1}{2\pi} tan^{-1}\left (\frac{y_2}{y_1} \right)
\end{align}

上記に注意して、

\begin{align}

\frac{\partial z_1}{\partial y_1} &= -y_1 exp\left( -\frac{y_1^2 + y_2^2}{2} \right) \\
\frac{\partial z_1}{\partial y_2} &= -y_2 exp\left( -\frac{y_1^2 + y_2^2}{2} \right) \\
\frac{\partial z_2}{\partial y_1} &= -\frac{1}{2\pi}\frac{y_2}{y_1^2 + y_2^2} \\
\frac{\partial z_2}{\partial y_2} &= \frac{1}{2\pi}\frac{y_1}{y_1^2 + y_2^2} \\
\end{align}

よって、ヤコビアンの絶対値は、

\begin{align}
\left| \frac{\partial(z_1, z_2)}{\partial (y_1, y_2)} \right|
&= \left| \frac{\partial z_1}{\partial y_1} \frac{\partial z_2}{\partial y_2} - \frac{\partial z_1}{\partial y_2} \frac{\partial z_2}{\partial y_1} \right| \\
&= \left| -\frac{1}{2\pi} exp\left( -\frac{y_1^2 + y_2^2}{2} \right) \right| \\
&= \left[ \frac{1}{\sqrt{2\pi}}exp\left(-\frac{y_1^2}{2} \right) \right]\left[ \frac{1}{\sqrt{2\pi}}exp\left(-\frac{y_2^2}{2} \right) \right]
\end{align}

従って、$y_1, y_2$の同時分布は、$p(z_1, z_2) = p(z_1)p(z_2) = 1$より、

\begin{align}
p(y_1, y_2) &= p(z_1, z_2) \left| \frac{\partial(z_1, z_2)}{\partial (y_1, y_2)} \right| \\
&= \left[ \frac{1}{\sqrt{2\pi}}exp\left(-\frac{y_1^2}{2} \right) \right]\left[ \frac{1}{\sqrt{2\pi}}exp\left(-\frac{y_2^2}{2} \right) \right]
\end{align}

以上より、$y_1, y_2$の同時分布は標準ガウス分布となる。

コード

gauss_pdf = lambda x: np.exp(-x**2 / 2) / np.sqrt(2 * np.pi)

z = np.random.uniform(0, 1, (1000, 2))
# 変換
y1 = np.sqrt(-2 * np.log(z[:,0])) * np.cos(2 * np.pi * z[:,1])
y2 = np.sqrt(-2 * np.log(z[:,0])) * np.sin(2 * np.pi * z[:,1])


fig, (axL, axR) = plt.subplots(ncols=2, figsize=(10,4), sharex=True)

# y1のヒストグラム
axL.hist(y1, bins=50, range=(-4, 4), color='aqua', label='y1')
# y2のヒストグラム
axR.hist(y2, bins=50, range=(-4, 4), color='aqua', label='y2')


# 標準ガウス分布のグラフ
x = np.linspace(-4, 4, 1000)
axL.plot(x, gauss_pdf(x)*170, 'darkblue', label='標準ガウス分布')
axR.plot(x, gauss_pdf(x)*170, 'darkblue', label='標準ガウス分布')

axL.legend(loc='best', prop={'size': 9})
axR.legend(loc='best', prop={'size': 9})

どっちの変数も標準ガウス分布に従ってますね。

まとめ

このセクションでは、求めたい分布の不定積分を計算して、その逆関数を求めることで確率変数を変換しました。
ただ、不定積分や逆関数を求めることが可能な分布は限られた数しかないです。
そこで、別のアプローチとして、次のセクションでは棄却サンプリングと重点サンプリングを扱います。