メトロポリス法 - 手帳と試行

M-H法の一種であるメトロポリス法を導入する。

M-H法は次のようなアルゴリズムであった。

Metropolis-Hastings Algorithm

\begin{darray}{rl} 1. & X \gets \{\} \\ 2. & \text{for } t \in \{1, 2, \dots, n\} : \\ 3. & \qquad y \sim Q(y | x) \\ 4. & \qquad r \sim {\cal U}(r | 0, 1) \\ 5. & \qquad \text{if } A(y | x) \gt r : \\ 6. & \qquad \qquad x \gets y \\ 7. & \qquad X \gets X \cup \{x\} \\ 8. & \operatorname{return} X \end{darray}

ただし、 $A(y | x)$ は受理確率と呼ばれる関数で、次式で表される。

\begin{aligned} A(y | x) = \frac{P(y) Q(x | y)}{P(x) Q(y | x)} \end{aligned}

ここで、提案分布 $Q(y | x)$ を、関数の台 $\mathcal Y$ が

\begin{gathered} \mathcal Y = \{x + \lambda c | \lambda \in [1, -1]\} \end{gathered}

であるような連続一様分布、すなわち

\begin{aligned} Q(y | x) = \mathcal U(y | x-c, x+c) = \frac{1}{2c} \mathbb I\{x-c\le y \le x+c\} \end{aligned}

とおくと、

\begin{aligned} Q(x | y) ={}& \mathcal U(x | y-c, y+c) \\ ={}& \frac{1}{2c} \mathbb I\{y-c\le x \le y+c\} \\ ={}& \frac{1}{2c} \mathbb I\{x-c\le y \le x+c\} \\ ={}& Q(y | x) \end{aligned}

ゆえ、受理確率は次式に変化する。

\begin{aligned} A(y | x) = \frac{P(y) Q(x | y)}{P(x) Q(y | x)} = \frac{P(y)}{P(x)} \end{aligned}

これを用いれば、MH法は次のようなアルゴリズムに変化する。

\begin{darray}{rl} 1. & X \gets \{\} \\ 2. & \text{for } t \in \{1, 2, \dots, n\} : \\ 3. & \qquad y \sim {\cal U}(y | x-c, x+c) \\ 4. & \qquad r \sim {\cal U}(r | 0, 1) \\ 5. & \qquad \text{if } \frac{P(y)}{P(x)} \gt r : \\ 6. & \qquad \qquad x \gets y \\ 7. & \qquad X \gets X \cup \{x\} \\ 8. & \operatorname{return} X \end{darray}

このようなMH法をとくにメトロポリス法 (Metropolis algorithm) という。定数 $c$ はステップ幅と呼ばれる。

このアルゴリズムをちょっと書き換えれば、 $x$ がスカラーではなくベクトル $\bm x \in \R^n$ であるような場合のメトロポリス法も容易に作成できる。

Metropolis Algorithm

\begin{darray}{rl} 1. & X \gets \{\} \\ 2. & \text{for } t \in \{1, 2, \dots, n\} : \\ 3. & \qquad \text{for } i \in \{1, 2, \dots, N\} : \\ 4. & \qquad \qquad y_i \sim {\cal U}(y_i | x_i-c_i, x_i+c_i) \\ 5. & \qquad r \sim {\cal U}(r | 0, 1) \\ 6. & \qquad \text{if } \frac{P(\bm y)}{P(\bm x)} \gt r : \\ 7. & \qquad \qquad \bm x \gets \bm y \\ 8. & \qquad X \gets X \cup \{\bm x\} \\ 9. & \operatorname{return} X \end{darray}

具体的には、提案 $\bm y \in \R^n$ の各要素 $y_i$ を一様分布 $\mathcal U(y_i | x_i - c_i, x_i + c_i)$ からサンプルする。ただし $c_i$ は $\bm c \in \R^n$ の各要素である。

\begin{aligned} y_i \sim \mathcal U(y_i | x_i - c_i, x_i + c_i) \end{aligned}

このとき以下のように、ステップ幅の各要素が時間によって変化しなければ、 $\bm c$ の全要素が同一である必要はない。ただしアルゴリズムの途中で要素 $c_i$ が変化してはいけないので注意が必要である。

\begin{aligned} &\text{valid}: && \bm c = \begin{bmatrix} 1 \\ 1 \end{bmatrix} && \to && \bm c = \begin{bmatrix} 1 \\ 1 \end{bmatrix} && \to && \dots \\ \\ &\text{valid}: && \bm c = \begin{bmatrix} 1.5 \\ 0.5 \end{bmatrix} && \to && \bm c = \begin{bmatrix} 1.5 \\ 0.5 \end{bmatrix} && \to && \dots \\ \\ &\text{invalid}: && \bm c = \begin{bmatrix} 1.5 \\ 0.5 \end{bmatrix} && \to && \bm c = \begin{bmatrix} 0.5 \\ 1 \end{bmatrix} && \to && \dots \\ \\ \end{aligned}

もういちどよく見てみよう。

Metropolis Algorithm

\begin{darray}{rl} 1. & X \gets \{\} \\ 2. & \text{for } t \in \{1, 2, \dots, n\} : \\ 3. & \qquad \text{for } i \in \{1, 2, \dots, N\} : \\ 4. & \qquad \qquad y_i \sim {\cal U}(y_i | x_i-c_i, x_i+c_i) \\ 5. & \qquad r \sim {\cal U}(r | 0, 1) \\ 6. & \qquad \text{if } \frac{P(\bm y)}{P(\bm x)} \gt r : \\ 7. & \qquad \qquad \bm x \gets \bm y \\ 8. & \qquad X \gets X \cup \{\bm x\} \\ 9. & \operatorname{return} X \end{darray}

確率分布を計算する箇所が、 $\dfrac{P(\bm y)}{P(\bm x)}$ のみになっている。

これはなかなかすごいことで、この式により、

\begin{aligned} P(\bm x) = \frac{1}{Z} S(\bm x) \end{aligned}

などの形で書かれていて、しかも正規化定数 $Z = \int d\bm x S(\bm x)$ の計算が手計算では到底不可能であるような場合でも、 $P(\bm x)$ に従う乱数が生成できるのである。