導入

昨今では、ベイズ的アプローチもはやっておらず、いかにして学習データを集めてDeep Learning にぶっこむかが流行っていますが、たまにはこういう古典的な計算をするのもありでしょう。

MCMC についても説明するつもりでしたが、単純なシミュレーションについて説明するのもそこそこしんどかったので、いったん単純なシミュレーションについて説明します。
MCMC についても、近いうちに書きたいなと個人的には思っています。

本記事では、大体 PRML の 11 章を参考にしています。

モンテカルロ法

機械学習に携わっていると、モンテカルロシミュレーション、もしくは、モンテカルロ法という単語をよく耳にすると思います。「あー、MCMCでしょ」とこの記事を書く前の僕は漠然と思っていたのですが、具体的に説明しようとすると答えに詰まってしまいました。こういうときは、Wikipedia さんに聞けば大体なことがわかりますので、とりあえず Wikipedia の記事見てみたいと思います。

モンテカルロ法（モンテカルロほう、（英: Monte Carlo method、MC）とはシミュレーションや数値計算を乱数を用いて行う手法の総称。
カジノで有名な国家モナコ公国の4つの地区（カルティ）の1つであるモンテカルロから名付けられた。

モンテカルロはただの地名でした。カジノという言葉を聞くと、カジノで勝とうとしている人たちがシミュレーションを行おうとして生み出したのかなと

計算理論の分野において、モンテカルロ法とは誤答する確率の上界が与えられる乱択アルゴリズム（ランダム・アルゴリズム）と定義される[1]。

つまり、乱数を用いれば大抵のことはモンテカルロシミュレーションやモンテカルロ法と名乗ることができるわけです。実際、「ミラー-ラビン素数判定法」という確率的に素数かどうかを判断するアルゴリズムもモンテカルロ法であると、記事には書いてあります。ですが、今回は機械学習に主眼を置いており、素数判定法なんぞにうつつを抜かしている暇はありません。

そのため、考える対象を以下に限定したいと思います。

ある確率分布 P を決めたときに、P の分布に従っている乱数 X を生成するアルゴリズムを考える

上の問題設定に出てくる乱数 $X$ を生成することができれば、様々なシミュレーションや数値計算に応用することができます。

目次

変数変換法
- 指数分布
- 正規分布
- コーシー分布
棄却サンプリング法
- ガンマ分布

全体を通して、$[0, 1]$ の一様分布の乱数 $U$ を生成することは可能であるとします。この $U$ から有名な確率分布に従う乱数を

変数変換法

生成したい確率分布を $P$ とします。
一様分布の乱数 $U$ に対して、「$Y = F_{P}(U)$ が $P$ の分布に従う」ような関数 $F_{P}$ を探す手法のことを変数変換法と言います。
「$Y = F_{P}(U)$ が $P$ の分布に従う」と書きましたが、これは以下の数式と同じ意味です。

$$ du = p(y) dy $$

$u, y$ が同時に出てきているため、わかったようなわからないような気分にさせてくれます。ですが、$Y = F_{P}(U)$ という関係があるので、$F_{p}$ の逆関数 $F_{p}^{-1}$ を用いると、$U = F_{p}^{-1}(Y)$ と書けます。この式を用いると、$p(y) = du / dy = d(F_{p}^{-1}(y)) / dy$ と書けることがわかります。
実際、この式変形はかなりトリッキーだと思います。少なくとも初見では理解しづらいです。なぜ、逆関数が現れるのか不思議だと思いますが、以下の例を見ながら少しずつ慣れていきましょう。

例: 指数分布

指数分布の確率密度関数 $p$ は以下の式で表されます。

$$p(x)=\frac{1}{\lambda} e^{-x/\lambda}$$

$G(x) = \int_{0}^{x} p(x)dx$ とおきます。積分した関数を微分すると元の関数に戻るので、$dG/dx = p(x)$ が成り立ちます。$p(y) = d(F_{p}^{-1}(y)) / dy$ を満たす関数を見つけたかったので、なんとなく $G$ の逆関数が $F$ である気がします (実際そうです)。
$G$ を具体的に計算すると、
$$ u = G(x) = \bigg[-\exp \bigg(-\frac{x}{\lambda}\bigg) \bigg]^{x} = 1 - \exp \bigg(-\frac{x}{\lambda}\bigg) $$
となるので、逆関数 $F$ は、
$$ F(u) = - \lambda \log(1 - u) $$
となります。
この関数 $F$ を用いて、乱数を生成する Python スクリプトを実際に書いてみました。

import numpy as np
import matplotlib.pyplot as plt

Lambda = 0.1
N = 100000
u = np.random.rand(N)
y = - Lambda * np.log(1 - u)

plt.hist(y, bins=200)
plt.show()

ヒストグラムの結果を見る限り、指数分布の確率密度関数と似た形になっており、指数分布に従う乱数が生成できていることがわかります。

例: 正規分布

次に、標準正規分布の乱数を作成する有名なアルゴリズムである「Box-Muller アルゴリズム」を紹介します。
半径 1 の円板上の一様分布 $(z_1, z_2)$ から

天下り的に式を与えたのは、この式を導出することを諦めてしまったからです。。

この関数 $F$ を用いて、乱数を生成する Python スクリプトを実際に書いてみました。

import numpy as np
import matplotlib.pyplot as plt

N = 100000
z = np.random.rand(N * 2, 2) * 2 - 1
z = z[z[:, 0] ** 2 + z[:, 1] ** 2 < 1]

z = z[:N]
r = z[:, 0] ** 2 + z[:, 1] ** 2
y = z[:, 0] / (r ** 0.5) * ((-2 * np.log(r)) ** 0.5)

plt.hist(y, bins=200)
plt.show()

コーシー分布

最後に、コーシー分布に従う乱数を生成してみます。
コーシー分布は指数分布や正規分布と比較すると、統計学的にもほとんど利用されません。ですが、次節のガンマ分布の項でコーシー分布に従う乱数を利用するため、ここで紹介します。
コーシー分布の確率密度関数は以下の式で表されます。

p(x; c) = \frac{1}{1 + (x - c)^2}

コーシー分布の累積分布関数 $G$ は $\arctan$ という $\tan$ 関数の逆関数を用いれば、以下のように簡単に書くことができます。

u = G(x) = \int_{-\infty}^{x} p(x; c) = \frac{1}{\pi} \arctan (x - c) + \frac{1}{2}

よって、$G$ の逆関数 $F$ は以下のように書けます。

F(u) = \tan \bigg(\pi \bigg(u - \frac{1}{2} \bigg)\bigg) + c

つまり、$$ F(u) = \tan \bigg(\pi \bigg(u - \frac{1}{2} \bigg)\bigg) + c$$
を $[0, 1]$ の一様乱数を $U$ とすると、$F(U)$ がコーシー分布に従う乱数であることがわかります。

この関数 $F$ を用いて、実際に乱数を生成した結果を見てみましょう。

cauchy_F = lambda x: np.tan(x * (np.pi / 2 + np.arctan(c)) - np.arctan(c)) + c

N = 100000
u = np.random.rand(N)
y = cauchy_F(u)
y = y[y < 10]

plt.hist(y[y < 10], bins=200)
plt.show()

ヒストグラムの結果を見る限り、指数分布の確率密度関数と似た形になっており、指数分布に従う乱数が生成できていることがわかります。
ですが、スクリプトをよく見てみると、y = y[y < 10] という気になる一行が入っています。
この一行を外して実行すると、とんでもないことになります。

つまり、とんでもなく大きい外れ値が生成される可能性があることを示しています。この結果は、コーシー分布の平均値が発散するという事実を思い出してみると、なんとなく納得できるかと思います。

乱数値を 10 以下で制限すると平均値は 2 前後かと思うのですが、きわめて大きい値がまれに発生してしまうために、平均値を上側にずり上げてしまっています。収入の平均値と中央値が乖離している話を思い出しますね。

変数変換法のポイント

分布の特性をよく知っている必要があります。
大体の場合は累積分布関数が簡単な形に書ける必要があります。
それ以外の方法で求めるのは相当骨が折れるかと思います。
また、F が簡単に計算できる関数でなければ、結局意味がありません。

棄却サンプリング法

例: ガンマ分布

ガンマ分布の確率密度関数は以下の式で表されます。

p(x; a, b) = b^{a} z^{a-1} \exp(-bz) / \Gamma(a)

$ \Gamma $ 関数の値を求めるのは大変なので、$\Gamma $ の部分を無視した $\tilde{p}$ を利用します。

\tilde{p}(x; a) = z^{a-1} \exp(-z)

ここで、コーシー分布を超える k を求めましょう。手計算で求めるのが面倒だったので、今回はスクリプトに頼りました。

gamma_p_tilda = lambda x: (x ** (a - 1)) * np.exp(-x)
cauchy_p = lambda x: 1 / (1 + ((x - c) ** 2)) / (np.pi / 2 + np.arctan(c))

z = np.arange(0, 100, 0.01)
gammas = gamma_p_tilda(z)
cauchys = cauchy_p(z)

k = np.max(gammas / cauchys)
cauchys *= k

plt.plot(z, gammas, label='gamma')
plt.plot(z, cauchys, label='cauchy')

plt.legend()

plt.show()

上述のスクリプトに格納されている $k$ の値は $2.732761944808582$ でした。

N = 100000
u = np.random.rand(N * 3)
y = cauchy_F(u)

u0 = np.random.rand(N * 3)

u0 *= k * cauchy_p(y)

y = y[u0 < gamma_p_tilda(y)]
y = y[:N]

plt.hist(y, bins=200)
plt.show()

棄却サンプリングの適用条件

乱数を簡単に生成することができる提案分布 $Q$ を用意する必要があります。
また、$kq \geq p$ を満たす $k$ を算出しなければいけません。$k$ の値を不用意に大きくしすぎると、せっかく生成した乱数を棄却しなければいけないため、できるかぎり $k$ をぎりぎりの小さい値に取りたいわけです。
しかし、そうなると関数の具体的な形から導出するか、精密に値を計算する必要があり、そこそこの苦労を伴います。
また、$k$ の値を決める場合や棄却するかを計算するときに、$P$, $Q$ ともに確率密度関数の値を使うため、確率密度関数自体が簡単に計算できる必要があります。そのため、複数の変数が複雑に絡み合ったモデル等、確率密度関数の計算自体が難しい場合には、この手法を適用することはできません。

MCMC に向けて

ここまで、標準的な乱数生成手法を概観してきました。
一通り見てきたところ、大元の確率分布がある程度扱いやすい形をしている必要があることが見て取れるかと思います。

正規分布や指数分布、ガンマ分布等、統計学的に重要な乱数が一様分布から一通り生成できることはわかりました。これだけでも十分いろいろなことができるかと思いますが、世の中にはベイズ的アプローチを用いたより複雑なモデルがうじゃうじゃ存在しています。
今回の生成手法だけでは、複雑なベイズ的モデルに対抗するだけの武器を持ち合わせてはいません。

そこで、MCMC と呼ばれる、より複雑なモデルに対して、乱数を生成することができるアルゴリズムが出てきます。
次回は、MCMC を適用する主領域であるグラフィカルモデルと MCMC の具体的な手法について説明したいと思います。