手帳と試行

学んだことをアウトプットしていきます。 日々、ノートあるのみ。

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

基本的なMCMC法であるメトロポリス・ヘイスティングス法を導入する。

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

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

準備

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

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

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

アルゴリズム

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

Metropolis-Hastings Algorithm
1.X{}2.for t{1,2,,n}:3.yQ(yx)4.rU(r0,1)5.if A(yx)>r:6.xy7.XX{x}8.returnX\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}

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

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

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

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

と表現できる。

注意事項

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

Q(yx)=NN(y0,I)exp ⁣(12y22)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(yx)A(y | x) をできるだけ 11 に近づけたいという要望もある。なぜなら、この確率が小さすぎると、同じサンプルを頻繁に出力してしまうようになるため、乱数として取り扱えるようになるのに時間がかかりすぎるからである。

make A(yx)=P(y)Q(xy)P(x)Q(yx) as large as possible!\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法というわけである。