メトロポリス・ヘイスティングス法

メトロポリス・ヘイスティングス法

M-H法を用いて、$p(x) \propto \exp(-x^4)$ からのサンプリングを行った例。$p(x)$ の分配関数を計算することなくサンプリングを行うことができる。

メトロポリス・ヘイスティングス法 (the Metropolis-Hastings algorithm; M-H法) はMCMC法の一つ。直接サンプリングすることが困難な確率分布からサンプルを生成するために使用される。「提案」という操作のアルゴリズムを工夫することで、メトロポリス法、ギブスサンプリング、ハミルトニアン・モンテカルロ法などに帰着する。M-H法はこれらの一般的な形であると解釈すれば良いだろう。

準備

まず「提案分布」と呼ばれる条件付き確率分布 $Q(y | x)$ を準備する。続いて「受理確率」と呼ばれる関数 $A(y | x)$ を準備する。この関数は次式で表される。

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

正確には、「確率」と呼ぶためには、これを $\min$ でラップして $P_\mathrm{accept} (y|x) = \min( A(y | x), 1 )$ とでもするべきだが、アルゴリズム的には何も変化しないので省略した。

アルゴリズム

確率分布 $P(x)$ に従う乱数列 $X$ を生成する。

$$\begin{darray}{rl} 1. & X \gets \{\} \\ 2. & \text{for } t \in \{1, 2, \dots, T\} : \\ 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}$$

なお、以上において $\mathcal U(0, 1)$ は $0$ 以上 $1$ 未満の連続一様分布を表す。また4. 以降については、次のような内容 (メトロポリステスト) を意味している。

• 確率 $\min( A(y | x), 1 )$ で提案を受理 (accept)、$x \gets y$。
• さもなくば棄却 (reject)、yを捨てる。

したがって、上記をもう少し端折って書くと

1. 提案: $y \sim Q(y | x)$
2. メトロポリステスト: $\min( A(y | x), 1 )$ の確率で提案を受理

と表現できる。

何がすごいの？

アルゴリズムをよく見ると、確率分布を計算する箇所が、$A(y | x) = \dfrac{P(y) Q(x | y)}{P(x) Q(y | x)}$ のみになっていることがわかる。

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

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

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

注意事項

このアルゴリズムは、提案分布 $Q(y | x)$ から乱数をサンプルする操作が要求される。したがって、可能な限りサンプルが容易な分布であることが望ましい。例えば次のような標準正規分布などはシンプルで使いやすい。

$$Q(\bm y | \bm x) = \mathcal N_N(\bm y | \bm 0, \bm I) \propto \exp \! \left( -\frac{1}{2} \| \bm y \|_2^2 \right)$$

ただし、単にシンプルなだけでは不十分で、受理確率 $A(y | x)$ をできるだけ $1$ に近づけたいという要望もある。なぜなら、この確率が小さすぎると、同じサンプルを頻繁に出力してしまうようになるため、乱数として取り扱えるようになるのに時間がかかりすぎるからである。

$$\begin{aligned} \text{make } A(y | x) = \frac{P(y) Q(x | y)}{P(x) Q(y | x)} \text{ as large as possible!} \end{aligned}$$

提案分布をシンプルにしつつ、受理確率も高くするようなものが、優れたM-H法というわけである。

実装例

スカラーをサンプリングするM-H法の Python での実装例。

from typing import Callable

class MetropolisHastings:
    def __init__(self,
                 target_distribution: Callable[[float], float],
                 proposal_distribution: Callable[[float, float], float],
                 proposal_sampler: Callable[[], float]
                 ):
        """
        Args:
            target_distribution (Callable): 目標とする確率分布
            proposal_distribution (Callable): 提案分布
            proposal_sampler (Callable): 提案分布からのサンプリングを行う関数
        """

        def calc_acceptance_ratio(proposed_x, current_x):
            """
            Args:
                proposed_x (float)
                current_x (float)

            Returns:
                float: 受理確率
            """
            acceptance_ratio = \
                (target_distribution(proposed_x) * proposal_distribution(current_x, proposed_x)) \
                / (target_distribution(current_x) * proposal_distribution(proposed_x, current_x))
            return acceptance_ratio

        self.calc_acceptance_ratio = calc_acceptance_ratio
        self.proposal_sampler = proposal_sampler
        self.samples = []

    def sample(self, target_samples, initial_x):
        """
        Args:
            target_samples (int): 目標サンプル数
            initial_x (float): サンプリングの初期値
        """
        self.samples.append(initial_x)
        for _ in range(target_samples - 1):
            current_x = self.samples[-1]
            proposed_x = self.proposal_sampler()

            acceptance_ratio = self.calc_acceptance_ratio(proposed_x, current_x)

            if np.random.rand() < acceptance_ratio:
                self.samples.append(proposed_x)
            else:
                self.samples.append(current_x)

提案分布 proposal_distribution() と、提案分布からのサンプリングを行う関数 proposal_sampler() は辻褄があっていないといけない。もしも一方が正規分布を指定するなら、もう一方も全く同じ正規分布を指定していなくてはいけない。もしも一方が連続一様分布なら、もう一方も連続一様分布を指定する。

具体的に、たとえば正規分布を指定する場合を見てみよう:

import numpy as np
import scipy.stats as stats

def target_distribution(x):
    return np.exp(-x**4)

def proposal_distribution(proposed_x, current_x):
    return stats.norm.pdf(proposed_x, loc=0, scale=1)

def proposal_sampler():
    return stats.norm.rvs(loc=0, scale=1)

np.random.seed(0)

initial_x = 0
iterations = 10000

mcmc_sampler = MetropolisHastings(target_distribution, proposal_distribution, proposal_sampler)
mcmc_sampler.sample(iterations, initial_x)

ここでは、提案分布として SciPy の stats.norm.pdf() (正規分布の確率密度関数) を指定し、提案分布からのサンプリングには stats.norm.rvs() (正規分布からのサンプリングを行う関数) を指定した。これをヒストグラムに可視化するには、たとえば次のようにする。

import matplotlib.pyplot as plt

true_x = np.linspace(-2, 2, 1000)
true_y = target_distribution(true_x)

fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(6,3))

ax1.hist(mcmc_sampler.samples, bins=30, density=True, color='w', linewidth=2, edgecolor='k')
ax1.set_xlabel('Sampled values')
ax1.set_ylabel('Frequency')
ax1.set_xlim(-2, 2)
ax1.set_ylim(bottom=0)
ax1.set_title('Sample')

ax2.plot(true_x, true_y, color='k')
ax2.set_xlabel('$x$')
ax2.set_xlim(-2, 2)
ax2.set_ylim(bottom=0)
ax2.set_title('Target distribution')

plt.show()

ここまでを実行すれば、冒頭で示した図のような結果が得られる。