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

提案分布 Q はランダムウォークの粒子が次に移動する候補点を提案する。

数学や物理において、メトロポリス・ヘイスティングス法(もしくは M-H アルゴリズム)(Metropolis-Hastings algorithm)　は直接サンプリングするのが難しい確率分布から統計標本の配列を生成するのに用いられるマルコフ連鎖を構築するのに用いられる手法である。この配列はマルコフ連鎖モンテカルロ法において、目標分布の近似（ヒストグラム）として用いられたり、期待値のような積分計算を必要とするものに用いられる。

歴史

このアルゴリズムは1953年にボルツマン分布のための特殊形で発表したニコラス・メトロポリスと 1970年にそれを一般化した W.K. ヘイスティングスにちなんで名づけられた。ギブスサンプリング法はM-H アルゴリズムの特殊形であり、通常より高速で、簡易に利用することができるが、応用範囲は狭まる。

アルゴリズムの説明

M-H アルゴリズムは、確率分布$\displaystyle P(x)$ の確率密度関数に比例する関数が計算できるならば、 いかなる確率分布 $\displaystyle P(x)$ からでもサンプルを得ることができる。

ベイジアンの応用においては正規化係数を計算するのが非常に難しいため、得たい確率分布に一致する必要がなく、その密度関数に比例すればいい点はM-H アルゴリズムの大きな長所である。

M-H アルゴリズムはサンプル列を生成する。 サンプルがあればあるほど、そのサンプル分布は目標の分布P(x)を近似することになる。 サンプルは反復して生成され、次のサンプルが生成される確率は現在のサンプルにのみ依存する。（このサンプル列はマルコフ連鎖である。）

このサンプルの選択方法がアルゴリズムの重要な点である。候補となるサンプルはある確率で採択か棄却が決まる。 採択とは候補となるサンプルが次のステップで使用されることで、棄却とは次のステップでも現在のサンプルが使用されることである。この採択確率はP(x)に関する現在のサンプルと次のサンプルの尤度を比較して決定される。

以下では簡単な説明のために、M-H アルゴリズムの１例であるメトロポリスアルゴリズムを示す。

メトロポリスアルゴリズム

$f(x)$を目標分布$P(x)$に比例する関数とする。

1. 初期化: 初期値$x_0$と任意の確率密度関数$Q(x|y)$を決める。
   • 関数$Q(x|y)$は過去のサンプル$y$が与えられたときに、候補サンプル$x$を生成する関数である。メトロポリスアルゴリズムでは$Q$は対称でなければならない。つまり$Q(x|y) = Q(y|x)$を満たさなければならない。よくある$Q(x|y)$の選択として$y$を平均としたガウス分布がある。これは$y$の近くは次も選択しやすいことになる。これによりサンプル列はランダムウォークとなる。$Q$は提案分布やjumping distributionと呼ばれる。
2. $t$番目の更新:
   • 確率分布$Q(x^*|x_t)$から次のサンプル$x^*$ を生成する。
   • 採択率$\alpha = \frac{f ( x^* )}{f( x_t )}$を計算する。サンプルを採択するか棄却するかを決定するために使用される。$P$の確率密度分布は$f$に比例しているため、$\alpha=\frac{f(x^*)}{f(x_t)}=\frac{P(x^*)}{P(x_t)}$である。
   • $\alpha \geq 1$である場合、現在の候補$x_t$より次の候補の方がありそうだということを示している。そのため候補は採択される$x_{t+1}=x^*$。そうでなければ確率$\alpha$で候補を採択する。候補が棄却されれば$x_{t+1}=x_t$とする。

このアルゴリズムは移動したり、ある点にとどまったりを繰り返しランダムにサンプル空間を動きまわる。採択率は、提案分布と現在のサンプルにしたがって生成されるサンプルがどれだけ採択されうるかを表す。現在のサンプルよりも次のサンプルの確率が高ければ、必ず次のサンプルを採択する。しかし次のサンプルの確率が高くない場合、ある確率で移動せずにサンプルを棄却する。このため、高い確率密度の領域からは多くのサンプルを、また低い確率の領域は少ないサンプルを生成することになる。直感的にはこれがアルゴリズムの仕組みであり、目的の確率分布に従ったサンプルを生成する方法である。

分布から直接独立サンプルを生成するM-H アルゴリズムや他のMCMCアルゴリズムなどの適応的棄却サンプリング と比べいくつかの短所がある。

1. • サンプルが相関していること。

長い時間サンプルを生成しても、近接したサンプルは相関をもち、分布を正しく反映したものではない。これは独立したサンプルを得たい場合、多くのサンプルを捨てなければならないことを意味し、$n$個飛ばしでサンプルを集めたりする。この$n$の値は通常は近接サンプル間の自己相関を計算し決められる。サンプルされる領域を広げると自己相関は減少させることができるが、提案されたサンプルの棄却確率を増加させることになる。適切なジャンプサイズではない場合、遅い混合のマルコフ連鎖となる。これはつまり、サンプルが高い相関をもつことで目標とする分布を精確に推定するために非常に多くのサンプルが必要となることを意味する。このアルゴリズムで生成されたマルコフ連鎖が目標の分布に収束されることは保証されているが、 初期値が低い確率の領域であったり、初期値が悪いとうまくいかない場合がある。その対策としてburn-inの期間は通常は必要である。そのため、始めの数千、万サンプルは捨てられる。

多くの棄却サンプリング法は次元の呪いの影響を受ける。次元の呪いとは棄却確率が次元数の関数として指数関数的に増加することである。他のMCMC法と同様にM-Hアルゴリズムではこの次元の問題を同程度には持たないため、サンプル空間の次元が高い場合には唯一の可能な解法となる。そのためMCMC法は、現在では多くの分野で使用されている階層ベイズモデルや高次元な統計モデルからのサンプルを生成する方法として選ばれている。

次元が高い場合には、個々の次元ごとに異なった振る舞いをすることや、遅い混合を避けるためにすべての次元に関して適切なジャンプの大きさを決定することが問題となるため、提案分布を適切に選択することが自体が困難である。 そのような状況でしばしば取られる代替案としてはギブスサンプリングがある。ギブスサンプリングは、すべての次元から一度にサンプルするのではなく、個々の次元に関してサンプリングを行う。 これは多くの典型な階層モデルにあるように、少数の変数が他の変数を条件付けている場合には有効な方法である。 個々の変数は他の変数に関して条件づけされて１度にサンプリングされる。 他には適応的棄却サンプリング、一次元M-H ステップ、スライスサンプリングが考えられる。

提案密度 $\displaystyle Q(x', x^t)$ が $\displaystyle x'$ に一切依存しないことも可能であり、その場合はアルゴリズムは「独立連鎖メトロポリス・ヘイスティングス法」という(それ以外は「酔歩連鎖メトロポリス・ヘイスティングス法」である)。 ふさわしい提案分布を持った独立連鎖M-Hアルゴリズムは酔歩連鎖法よりも高い精度があるが、目標分布に関してアプリオリの知識を必要とする。

定式化

M-H アルゴリズムは目標確率分布$P(x)$に従ったサンプルの生成を行うことが目的である。 これを達成するために、漸近的に唯一の定常分布π(x)に収束するマルコフ連鎖を用いる。

マルコフ過程は、２つの状態間の遷移確率$P(x\rightarrow x')$によって一意に定義される。 ２つの条件が満たされるとき、マルコフ過程は定常分布$\pi(x)$をもつ。

1. 定常分布の存在：定常である分布$\pi(x)$が存在しなければならない。全ての状態ペア$x$、$x'$について、遷移$x\rightarrow x'$が逆にも遷移可能である。これは詳細釣り合い条件で保証される。詳細釣り合いとは、状態$x$から状態$x'$への遷移確率は状態$x'$から状態$x$への遷移確率と等しいこと、つまり、$\pi(x)P(x\rightarrow x') = \pi(x')P(x'\rightarrow x)$である。
2. 定常分布の一意性: 定常分布$\pi(x)$は一意でなければならない。

これはマルコフ過程のエルゴード性から保証される。エルゴード性はすべての状態が非周期的（aperiodic）ででありかつ正再帰でなければならない。非周期とは系が固定距離である状態にもどることができないことであり、正再帰とは非ゼロの確率である全ての状態に遷移可能である。

M-H アルゴリズムは遷移確率の構成により、上記の２つの条件を満たすようにマルコフ過程を設計することができる。

導出を詳細釣り合い条件から始める、

$P(x)P(x\rightarrow x') = P(x')P(x'\rightarrow x)$

これは、以下のように書き換えられる。

$\frac{P(x\rightarrow x')}{P(x'\rightarrow x)} = \frac{P(x')}{P(x)}$.

通常の手法として遷移確率を提案確率分布と採択確率分布に分解する。提案分布$\displaystyle g(x\rightarrow x')$は$x$が与えれたときの状態$x'$の提案する条件付き確率であり、採択確率$\displaystyle A(x\rightarrow x')$は$x$が与えれたときの状態$x'$を採択する条件付き確率である。

$P(x\rightarrow x') = g(x\rightarrow x') A(x\rightarrow x')$

この関係を以前の式に代入して以下の式を得る。

$\frac{A(x\rightarrow x')}{A(x'\rightarrow x)} = \frac{P(x')}{P(x)}\frac{g(x'\rightarrow x)}{g(x\rightarrow x')}$ .

次のステップとして以前の式を満たす採択率を選ぶことが必要である。よくある選択として、メトロポリス選択が知られ以下の式で得られる。この値はアルゴリズムの実装に必要な値である。

$A(x\rightarrow x') = \min\left(1,\frac{P(x')}{P(x)}\frac{g(x'\rightarrow x)}{g(x\rightarrow x')}\right)$

実装の観点からはMetropolis–Hastings アルゴリズムは以下のステップから成り立っている。

1. 初期化： ランダムに$x$を設定する
2. $\displaystyle g(x\rightarrow x')$に従い$x'$を生成する
3. $\displaystyle A(x\rightarrow x')$に従い採択し$x'$に遷移する。採択されない場合は、$x' = x$となり値を更新しない。
4. $T$回以下であれば２に戻る
5. 値を保存する。２に戻る。

サンプルを適切に集めるためには、$T$は提案分布や採択率とが別に決められ、ステップ４においてサンプルが相関していないことが必要である。マルコフ過程の自己相関時間の時間のオーダーによる。

一般的にこのパラメータの決定は簡単ではないことは重要な点である。問題に対して適切にパラメータは決定されるべきである。分布に関する知識が全くない場合には一様分布が提案分布として選ばれることもある。この場合、状態$x$と状態</math>x'</math>はいつも相関しないために$T$の値は$1$に設定される。

アルゴリズムの手順

現時点のサンプル値は $x^t\,$ であるとする。 メトロポリス・ヘイスティングス法では、次のサンプル$x'$は 確率密度関数 $Q(x'; x^t)\,$ に従い生成される。 そして以下の値を計算する。

$a = a_1 a_2\,$

ここで、

$a_1 = \frac{P(x')}{P(x^t)} \,\!$

はサンプルの候補 $x'\,$ とその一つ前のサンプル $x^t\,$ の尤度比であり、

$a_2 = \frac{Q( x^t; x' )}{Q(x';x^t)}$

は２つの提案分布($x^t\,$ から $x'\,$ へとその逆方向)の比率である。 提案分布が状態に関して対称の場合は $\displaystyle a_2$ は 1 となる。

次に、以下のルールをもとづき新しい状態 $x^{t+1}\,$ が選択される。

$\displaystyle a \geq 1$ の場合:

$\displaystyle x^{t+1} = x'$

$\displaystyle a < 1$ の場合:

$\displaystyle a$ の確率で $\displaystyle x^{t+1} = x'$

$\displaystyle 1 - a$ の確率で $\displaystyle x^{t+1} = x^t$

マルコフ連鎖はランダムな初期値 $\displaystyle x^0$ から開始され、初期値が「忘れられる」まで、多くの試行を繰り返す。この間の標本は棄てられ、burn-in(機械や回路を通電した直後の安定しない状態からの比喩、初期通電)とよばれる。 burn-in'後の値 $x$ の集合は、分布 $P(x)$ からのサンプルを表す。

M-Hアルゴリズムは提案分布 $Q(x'; x^t)$の形が、直接サンプリングが困難である目標分布 $\displaystyle P(x)$ の形と類似している場合、つまり $Q(x'; x^t) \approx P(x') \,\!$ のときにうまくアルゴリズムが動作する。 しかし、多くの場合目標分布の形は未知である。

ガウス分布の提案分布 $\displaystyle Q$ が用いられる場合は、分散パラメータの $\displaystyle \sigma^2$ が burn-in 期間のうちに調整される必要がある。これは通常、採択率を計算することで行われる。採択率とは $\displaystyle N$個のサンプルのうち採択されたサンプルの割合である。要求される採択率は目標分布によって異なる。理論的には、一次元ガウス分布を目標分布とすると理想的な採択率は約50% であり、N次元ガウス分布の目標分布では約 23% であることが知られている。

$\displaystyle \sigma^2$ が小さすぎるとサンプリング列は低速混合をおこす。 つまり採択率が高くなり標本空間の移動が遅くなる。 そのため $\displaystyle P(x)$ への収束が遅くなる。 逆に $\displaystyle \sigma^2$ が大きすぎると、 採択率が低すぎサンプルが確率密度の低い所に移動してしまい $\displaystyle a_1$ が非常に小さくなってしまう。 このため $\displaystyle P(x)$ への収束が遅くなる。 したがって、提案分布のパラメータを調整し採択率を適切にする必要がある。

関連記事

• 焼きなまし法
• メルセンヌ・ツイスタ：もっとも実用的な乱数生成方法
• 詳細釣り合い
• 正再帰性
• マルコフ連鎖モンテカルロ法

M-H アルゴリズムの一例として

• ギブスサンプリング；条件付き分布から１変数ごとのサンプリング
• メトロポリス法
• スライスサンプリング
• ハイブリッドモンテカルロ
• MTM アルゴリズム

MCMCではない方法

• 逆関数法　求めたい分布の確率密度関数の逆関数が得られる場合
• 棄却サンプリング　求めたい分布に比例する関数が得られる場合
• 適応的棄却サンプリング　適応的に包絡関数を決め棄却を少なくする方法
• 重点サンプリング　重要度重みを考え、期待値を得るための方法
• SIR sampling-importance-resampling

モンテカルロEMアルゴリズム：EステップをサンプリングにしたEMアルゴリズム。

• メトロポリス光輸送：応用として

関連書籍

• Bernd A. Berg. Markov Chain Monte Carlo Simulations and Their Statistical Analysis. Singapore, World Scientific 2004.
• Siddhartha Chib and Edward Greenberg: "Understanding the Metropolis–Hastings Algorithm". American Statistician, 49(4), 327–335, 1995
• W.K. Hastings. "Monte Carlo Sampling Methods Using Markov Chains and Their Applications", Biometrika, 57(1):97-109, 1970. JSTOR doi:10.2307/2334940
• N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller. "Equations of State Calculations by Fast Computing Machines". Journal of Chemical Physics, 21(6):1087-1092, 1953. Equation of State Calculations by Fast Computing Machines | Browse - Journal of Chemical Physics

外部リンク

• MathWorld article on Simulated Annealing
• Metropolis-Hastings algorithm on xβ