メトロポリス・ヘイスティングスアルゴリズム

統計学および統計物理学において、メトロポリス・ヘイスティングス法は、直接サンプリングが難しい確率分布からランダムサンプルのシーケンスを取得するためのマルコフ連鎖モンテカルロ(MCMC) 法です。新しいサンプルは 2 つのステップでシーケンスに追加されます。まず、前のサンプルに基づいて新しいサンプルが提案され、次に、その時点での確率分布の値に応じて、提案されたサンプルがシーケンスに追加されるか拒否されます。結果のシーケンスは、分布を近似する (例:ヒストグラムを生成する) ため、または積分を計算する(例:期待値) ために使用できます。

メトロポリス・ヘイスティングスやその他の MCMC アルゴリズムは、一般的に多次元分布からのサンプリングに使用され、特に次元数が多い場合に使用されます。1 次元分布の場合、通常、分布から独立したサンプルを直接返すことができる他の方法 (適応拒否サンプリングなど) があり、これらはMCMC 方法に固有の自己相関サンプルの問題がありません。

このアルゴリズムは、1953年にアリアナ・W・ローゼンブルース、マーシャル・ローゼンブルース、オーガスタ・H・テラー、エドワード・テラーと共著した「高速コンピューティングマシンによる状態方程式の計算」という論文の第一共著者であるニコラス・メトロポリスにちなんで名付けられました。長年、このアルゴリズムは単にメトロポリスアルゴリズムとして知られていました。^{[1] [2]}この論文では対称提案分布の場合のアルゴリズムが提案されましたが、1970年にWKヘイスティングスがそれをより一般的なケースに拡張しました。^[3]この一般化された方法は最終的に両方の名前で識別されましたが、「メトロポリス-ヘイスティングスアルゴリズム」という用語が最初に使用されたのは明らかではありません。

メトロポリス アルゴリズムの開発の功績については、議論がある。この手法の計算面に精通していたメトロポリスは、スタニスワフ ウラムと共著した以前の論文で「モンテ カルロ」という用語を作り出し、1952 年の実験で使用されたMANIAC Iコンピュータを設計および構築した理論部門のグループを率いた。しかし、2003 年以前には、このアルゴリズムの開発に関する詳細な説明はなかった。マーシャル ローゼンブルースは、死の直前に、 1953 年の出版 50 周年を記念する LANL での 2003 年の会議に出席した。この会議で、ローゼンブルースは「統計力学のためのモンテ カルロ アルゴリズムの起源」と題するプレゼンテーションで、このアルゴリズムとその開発について説明した。^{[4]グベルナティスは、50 周年記念会議について 2005 年の雑誌記事[5]} でさらに歴史的に説明している。ローゼンブルース氏は、自分と妻のアリアナがこの作業を行ったこと、そしてメトロポリスはコンピューターの使用時間を提供する以外、開発に何ら関与していないことを明確にしている。

これはエドワード・テラーの記述と矛盾する。テラーは回想録の中で、1953年の論文の5人の著者は「何日も（何晩も）」一緒に作業したと述べている。^[6]対照的に、ローゼンブルースの詳細な記述では、テラーが「詳細な運動学に従う代わりに、統計力学を利用して集団平均を取る」という重要だが初期の提案をしたとしている。ローゼンブルースによれば、これが一般化モンテカルロ法について考えるきっかけとなり、ジョン・フォン・ノイマンとよく話し合ったという。アリアナ・ローゼンブルースは（2003年にグベルナティスに）オーガスタ・テラーがコンピューターの仕事を始めたが、アリアナ自身がそれを引き継いでコードを一から書いたと語った。死の直前に記録された口述歴史の中で、^[7]ローゼンブルースは再び、テラーが最初の問題を提起し、テラー自身がそれを解決し、アリアナがコンピューターをプログラミングしたとしている。

メトロポリス-ヘイスティングス アルゴリズムは、確率密度 を持つ任意の確率分布からサンプルを抽出できます。ただし、密度に比例する関数がわかっていて、 の値が計算できる必要があります。は密度と正確に等しいのではなく、密度に比例するだけでよいという要件により、メトロポリス-ヘイスティングス アルゴリズムは特に便利です。これは、実際には非常に困難なことが多い密度の正規化係数を計算する必要がないためです。 ${\displaystyle P(x)}$${\displaystyle f(x)}$ ${\displaystyle P}$${\displaystyle f(x)}$${\displaystyle f(x)}$

メトロポリス-ヘイスティングス アルゴリズムは、サンプル値のシーケンスを生成する際に、生成されるサンプル値が増えるにつれて、値の分布が目的の分布に近づくようにします。これらのサンプル値は、次のサンプルの分布が現在のサンプル値のみに依存するように反復的に生成され、サンプルのシーケンスはマルコフ連鎖になります。具体的には、各反復で、アルゴリズムは現在のサンプル値に基づいて次のサンプル値の候補を提案します。次に、ある確率で、候補が受け入れられるか (この場合、候補値は次の反復で使用されます)、または拒否されます (この場合、候補値は破棄され、次の反復で現在の値が再利用されます)。受け入れられる確率は、目的の分布に関して現在のサンプル値と候補サンプル値の関数の値を比較することによって決定されます。 ${\displaystyle f(x)}$

新しい候補を提案するために使用される方法は、前のサンプル が与えられた場合の、新しく提案されたサンプル の確率分布( と書かれることもある) によって特徴付けられます。これは、提案密度、提案関数、またはジャンピング分布と呼ばれます。 の一般的な選択はを中心とするガウス分布であり、これにより、 に近い点が次に訪問される可能性が高くなり、サンプルのシーケンスがガウスランダムウォークになります。Metropolis ら (1953) による最初の論文では、は からの最大距離に制限された一様分布であると示唆されました。ハミルトンモンテカルロ、ランジュバンモンテカルロ、または前処理されたクランクニコルソンなどのより複雑な提案関数も可能です。 ${\displaystyle g(x\mid y)}$${\displaystyle Q(x\mid y)}$${\displaystyle x}$${\displaystyle y}$${\displaystyle g(x\mid y)}$${\displaystyle y}$${\displaystyle y}$${\displaystyle g(x\mid y)}$${\displaystyle y}$

説明のために、提案関数が対称であるメトロポリス-ヘイスティングス アルゴリズムの特殊なケースであるメトロポリス アルゴリズムについて以下に説明します。

メトロポリスアルゴリズム（対称提案分布）

を目的の確率密度関数（目標分布とも呼ばれる）^[a]に比例する関数とします。 ${\displaystyle f(x)}$${\displaystyle P(x)}$

1. 初期化:サンプルの最初の観測点となる任意の点を選択し、提案​​関数 を選択します。このセクションでは、は対称であると仮定します。つまり、 を満たす必要があります。${\displaystyle x_{t}}$${\displaystyle g(x\mid y)}$${\displaystyle g}$${\displaystyle g(x\mid y)=g(y\mid x)}$
2. 各反復tについて:
   • 分布から選択して次のサンプルの候補を提案します。${\displaystyle x'}$${\displaystyle g(x'\mid x_{t})}$
   • ^候補[b]を受け入れるか拒否するかを決定するために使用される受け入れ率 を計算します。f はPの密度に比例するため、 となります。${\displaystyle \alpha =f(x')/f(x_{t})}$${\displaystyle \alpha =f(x')/f(x_{t})=P(x')/P(x_{t})}$
   • 承認または拒否:
     • 一様乱数を生成します。${\displaystyle u\in [0,1]}$
     • の場合、 を設定して候補を受け入れます。${\displaystyle u\leq \alpha }$${\displaystyle x_{t+1}=x'}$
     • の場合、候補を拒否し、代わりに を設定します。${\displaystyle u>\alpha }$${\displaystyle x_{t+1}=x_{t}}$

このアルゴリズムは、サンプル空間内をランダムに移動しようとすることで進行し、移動を受け入れることもあれば、その場に留まることもあります。特定のポイントでの は、アルゴリズムによってポイントに費やされた反復に比例します。受け入れ率は、密度が である分布に従って、現在のサンプルに対して提案された新しいサンプルがどの程度ありそうかを示していることに注意してください。既存のポイントよりも可能性の高いポイント (つまり、に対応するの高密度領域内のポイント) に移動しようとすると、常にその移動を受け入れます。ただし、可能性の低いポイントに移動しようとすると、移動を拒否することがあり、確率の相対的な低下が大きいほど、新しいポイントを拒否する可能性が高くなります。したがって、 の高密度領域にとどまり (そしてそこから多数のサンプルを返し) 、低密度領域を訪れるのはたまにだけになる傾向があります。直感的に、これがこのアルゴリズムが機能し、密度 で目的の分布に従うサンプルを返す理由です。 ${\displaystyle P(x)}$${\displaystyle x}$${\displaystyle \alpha}$${\displaystyle P(x)}$${\displaystyle P(x)}$${\displaystyle \alpha >1\geq u}$${\displaystyle P(x)}$${\displaystyle P(x)}$

分布から独立したサンプルを直接生成する 適応的除去サンプリング^[8]のようなアルゴリズムと比較すると、メトロポリス-ヘイスティングスや他のMCMCアルゴリズムにはいくつかの欠点がある。

• サンプルは自己相関しています。長期的には正しく に従いますが、近くのサンプルのセットは互いに相関し、分布を正しく反映しません。つまり、有効なサンプル サイズは実際に取得されたサンプルの数よりも大幅に少なくなる可能性があり、大きな誤差につながります。${\displaystyle P(x)}$
• マルコフ連鎖は最終的には望ましい分布に収束しますが、特に開始点が低密度領域にある場合、初期のサンプルは非常に異なる分布に従う可能性があります。その結果、通常はバーンイン期間が必要になり、^[9]初期サンプル数を破棄します。

一方、ほとんどの単純な棄却サンプリング法は、「次元の呪い」に悩まされており、棄却の確率は次元の数に応じて指数関数的に増加します。メトロポリス-ヘイスティング法は、他の MCMC 法と同様に、この問題がそれほど発生しないため、サンプリングする分布の次元数が多い場合に利用できる唯一のソリューションとなることがよくあります。その結果、MCMC 法は、階層型ベイズ モデルや、現在多くの分野で使用されているその他の高次元統計モデルからサンプルを生成するための選択肢となることがよくあります。

多変量分布では、上で説明した古典的なメトロポリス-ヘイスティングス アルゴリズムでは、新しい多次元サンプル ポイントを選択します。次元の数が多い場合、使用する適切なジャンプ分布を見つけるのは困難です。これは、個々の次元が非常に異なる方法で動作し、ジャンプ幅 (上記を参照) がすべての次元で同時に「ちょうど良い」ものでなければ、混合が過度に遅くなるのを防ぐことができないためです。このような状況でより適切に機能する代替アプローチは、ギブス サンプリングと呼ばれ、一度にすべての次元のサンプルを選択するのではなく、各次元の新しいサンプルを他の次元とは別に選択します。この方法では、潜在的に高次元の空間からサンプリングするという問題は、小さな次元からサンプリングするという一連の問題に縮小されます。^[10]これは、多変量分布が、各変数が他の少数の変数のみに条件付けられている一連の個別のランダム変数で構成されている場合、特に当てはまります。これは、ほとんどの典型的な階層モデルの場合に当てはまります。個々の変数は、一度に 1 つずつサンプリングされ、各変数は他のすべての変数の最新の値に条件付けられます。多変量分布の正確な形に応じて、これらの個々のサンプルを選択するためにさまざまなアルゴリズムを使用できます。いくつかの可能性としては、適応拒否サンプリング法、^[8]適応拒否メトロポリスサンプリングアルゴリズム、^[11]単純な1次元メトロポリス-ヘイスティングスステップ、またはスライスサンプリングがあります。

メトロポリス・ヘイスティングスアルゴリズムの目的は、望ましい分布に従って状態の集合を生成することである。これを実現するために、アルゴリズムはマルコフ過程を使用し、漸近的にとなる唯一の定常分布に到達する。^[12]${\displaystyle P(x)}$ ${\displaystyle \pi (x)}$${\displaystyle \pi (x)=P(x)}$

マルコフ過程は遷移確率、つまり任意の状態から他の任意の状態へ遷移する確率によって一意に定義される。次の2つの条件が満たされる場合、マルコフ過程は一意に定常分布を持つ。 ^[12]${\displaystyle P(x'\mid x)}$${\displaystyle x}$${\displaystyle x'}$${\displaystyle \pi (x)}$

1. 定常分布の存在: 定常分布が存在する必要があります。十分な条件ではあるが必要条件ではないのが詳細なバランスであり、これは各遷移が可逆的であることを必要とします。つまり、すべての状態ペアについて、状態 にあり状態 に遷移する確率は、状態 にあり状態 に遷移する確率と等しくなければなりません。${\displaystyle \pi (x)}$${\displaystyle x\to x'}$${\displaystyle x,x'}$${\displaystyle x}$${\displaystyle x'}$${\displaystyle x'}$${\displaystyle x}$${\displaystyle \pi (x)P(x'\mid x)=\pi (x')P(x\mid x')}$
2. 定常分布の一意性: 定常分布は一意でなければなりません。これは、マルコフ過程のエルゴード性によって保証されます。エルゴード性では、すべての状態が (1) 非周期的 (システムは固定間隔で同じ状態に戻らない) であり、(2) 正の再帰的 (同じ状態に戻るための期待ステップ数が有限) である必要があります。${\displaystyle \pi (x)}$

メトロポリス・ヘイスティングスアルゴリズムでは、上記の 2 つの条件を満たすマルコフ過程を (遷移確率を構築することによって) 設計し、その定常分布がとなるように選択します。アルゴリズムの導出は、詳細なバランスの条件から始まります。 ${\displaystyle \pi (x)}$${\displaystyle P(x)}$

${\displaystyle P(x'\mid x)P(x)=P(x\mid x')P(x'),}$

これを書き直すと

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

このアプローチでは、遷移を提案と承認-拒否の 2 つのサブステップに分離します。提案分布は、与えられた状態を提案する条件付き確率であり、承認分布は、提案された状態を受け入れる確率です。遷移確率は、それらの積として表すことができます。 ${\displaystyle g(x'\mid x)}$${\displaystyle x'}$${\displaystyle x}$${\displaystyle A(x',x)}$${\displaystyle x'}$

${\displaystyle P(x'\mid x)=g(x'\mid x)A(x',x).}$

この関係を前の式に代入すると、

${\displaystyle {\frac {A(x',x)}{A(x,x')}}={\frac {P(x')}{P(x)}}{\frac {g(x\mid x')}{g(x'\mid x)}}.}$

導出の次のステップは、上記の条件を満たす受け入れ率を選択することです。一般的な選択肢の 1 つは、メトロポリスの選択です。

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

このメトロポリスの受入比率については、またはのいずれの場合も、条件は満たされます。 ${\displaystyle A}$${\displaystyle A(x',x)=1}$${\displaystyle A(x,x')=1}$

したがって、メトロポリス-ヘイスティングスアルゴリズムは次のように記述できます。

1. 初期化
   1. 初期状態を選択します。${\displaystyle x_{0}}$
   2. セット。${\displaystyle t=0}$
2. 繰り返し
   1. に従ってランダムな候補状態を生成します。${\displaystyle x'}$${\displaystyle g(x'\mid x_{t})}$
   2. 受け入れ確率を計算します。${\displaystyle A(x',x_{t})=\min \left(1,{\frac {P(x')}{P(x_{t})}}{\frac {g(x_{t}\mid x')}{g(x'\mid x_{t})}}\right)}$
   3. 承認または拒否:
      1. 一様乱数を生成する。${\displaystyle u\in [0,1]}$
      2. の場合、新しい状態を受け入れてを設定します。${\displaystyle u\leq A(x',x_{t})}$${\displaystyle x_{t+1}=x'}$
      3. の場合、新しい状態を拒否し、古い状態をコピーして転送します。${\displaystyle u>A(x',x_{t})}$${\displaystyle x_{t+1}=x_{t}}$
   4. 増分: 設定します。${\displaystyle t=t+1}$

指定された条件が満たされれば、保存された状態の経験的分布はに近づきます。効果的に推定するために必要な反復回数 ( ) は、と提案分布の関係や推定の望ましい精度など、要因の数によって異なります。 ^[13]離散状態空間上の分布の場合、マルコフ過程の自己相関時間 のオーダーでなければなりません。 ^[14]${\displaystyle x_{0},\ldots ,x_{T}}$${\displaystyle P(x)}$${\displaystyle T}$${\displaystyle P(x)}$${\displaystyle P(x)}$

一般的な問題では、どの分布を使用するか、または適切な推定に必要な反復回数は明確ではないことに注意することが重要です。どちらも方法の自由なパラメータであり、特定の問題に合わせて調整する必要があります。 ${\displaystyle g(x'\mid x)}$

メトロポリス・ヘイスティングス法の一般的な用途は積分を計算することである。具体的には、空間と、、上の確率分布を考える。メトロポリス・ヘイスティングス法は、次の形式の積分を推定することができる。 ${\displaystyle \Omega \subset \mathbb {R} }$${\displaystyle P(x)}$${\displaystyle \Omega }$${\displaystyle x\in \Omega }$

${\displaystyle P(E)=\int _{\Omega }A(x)P(x)\,dx,}$

ここで、 は（測定可能な）関心関数です。 ${\displaystyle A(x)}$

たとえば、統計量 とその確率分布（周辺分布）を考えてみましょう。目標はの裾におけるを推定することだとします。正式には、は次のように記述できます 。${\displaystyle E(x)}$${\displaystyle P(E)}$${\displaystyle P(E)}$${\displaystyle E}$${\displaystyle P(E)}$${\displaystyle P(E)}$

${\displaystyle P(E)=\int _{\Omega }P(E\mid x)P(x)\,dx=\int _{\Omega }\delta {\big (}E-E(x){\big )}P(x)\,dx=E{\big (}P(E\mid X){\big )}}$

したがって、推定は、指標関数の期待値を推定することによって実行できます。この値は、 のときは 1、それ以外のときは 0 です。は の裾にあるため、の裾にがある状態を描く確率はに比例しますが、これは定義により小さい値です。ここで、メトロポリス・ヘイスティングス アルゴリズムを使用して、（まれな）状態をより確実にサンプリングし、裾の推定に使用するサンプル数を増やすことができます。これは、たとえば、 の状態（など）を優先するサンプリング分布を使用することによって実行できます。 ${\displaystyle P(E)}$ ${\displaystyle A_{E}(x)\equiv \mathbf {1} _{E}(x)}$${\displaystyle E(x)\in [E,E+\Delta E]}$${\displaystyle E}$${\displaystyle P(E)}$${\displaystyle x}$${\displaystyle E(x)}$${\displaystyle P(E)}$${\displaystyle P(E)}$${\displaystyle P(E)}$${\displaystyle \pi (x)}$${\displaystyle \pi (x)\propto e^{aE}}$${\displaystyle a>0}$

最も最近サンプリングされた値が であると仮定します。メトロポリス・ヘイスティングスアルゴリズムに従うために、次に確率密度を持つ新しい提案状態を描き、値を計算し ${\displaystyle x_{t}}$${\displaystyle x'}$${\displaystyle g(x'\mid x_{t})}$

${\displaystyle a=a_{1}a_{2},}$

どこ

${\displaystyle a_{1}={\frac {P(x')}{P(x_{t})}}}$

は提案されたサンプルと前のサンプルの間の確率比（例えばベイズ事後確率）であり、 ${\displaystyle x'}$${\displaystyle x_{t}}$

${\displaystyle a_{2}={\frac {g(x_{t}\mid x')}{g(x'\mid x_{t})}}}$

は、2 つの方向 (からへ、およびその逆)の提案密度の比率です。提案密度が対称の場合、この値は 1 になります。新しい状態は、次の規則に従って選択されます。 ${\displaystyle x_{t}}$${\displaystyle x'}$${\displaystyle x_{t+1}}$

もし${\displaystyle a\geq 1{:}}$

${\displaystyle x_{t+1}=x',}$

それ以外：

${\displaystyle x_{t+1}={\begin{cases}x'&{\text{with probability }}a,\\x_{t}&{\text{with probability }}1-a.\end{cases}}}$

マルコフ連鎖は任意の初期値 から開始され、この初期状態が「忘れられる」までアルゴリズムは多くの反復で実行されます。破棄されるこれらのサンプルはバーンインと呼ばれます。 の受け入れられた値の残りのセットは、分布 からのサンプルを表します。 ${\displaystyle x_{0}}$${\displaystyle x}$${\displaystyle P(x)}$

このアルゴリズムは、提案密度がターゲット分布の形状と一致する場合、つまり直接サンプリングが難しい場合に最もよく機能します。ガウス提案密度を使用する場合、分散パラメータはバーンイン期間中に調整する必要があります。これは通常、受け入れ率を計算することによって行われます。受け入れ率は、最後のサンプルのウィンドウで受け入れられる提案サンプルの割合です。望ましい受け入れ率はターゲット分布によって異なりますが、1次元ガウス分布の理想的な受け入れ率は約50％で、次元ガウスターゲット分布の場合は約23％に低下することが理論的に示されています。^{[15]これらのガイドラインは、十分に規則的なベイズ事後分布からサンプリングする場合にうまく機能します。これは、} Bernstein–von Misesの定理を使用して確立できる多変量正規分布に従うことが多いためです。^[16]${\displaystyle P(x)}$${\displaystyle g(x'\mid x_{t})\approx P(x')}$${\displaystyle g}$${\displaystyle \sigma ^{2}}$${\displaystyle N}$${\displaystyle N}$

が小さすぎる場合、チェーンの混合はゆっくりになります（つまり、受け入れ率は高くなりますが、後続のサンプルは空間内をゆっくりと移動し、チェーンは にゆっくりと収束します）。一方、が大きすぎる場合、提案は確率密度のはるかに低い領域に着地する可能性が高くなるため受け入れ率は非常に低くなり、 は非常に小さくなり、チェーンはやはり非常にゆっくりと収束します。通常、提案の分布を調整して、前の段落で述べた理論的な推定値に沿って、アルゴリズムが全サンプルの約 30% を受け入れるようにします。 ${\displaystyle \sigma ^{2}}$${\displaystyle P(x)}$${\displaystyle \sigma ^{2}}$${\displaystyle a_{1}}$

MCMC は、統計モデルの事後分布からサンプルを抽出するために使用できます。受け入れ確率は次のように与えられます。 ここで、は尤度、は事前確率密度、は (条件付き) 提案確率です。 ${\displaystyle P_{acc}(\theta _{i}\to \theta ^{*})=\min \left(1,{\frac {{\mathcal {L}}(y|\theta ^{*})P(\theta ^{*})}{{\mathcal {L}}(y|\theta _{i})P(\theta _{i})}}{\frac {Q(\theta _{i}|\theta ^{*})}{Q(\theta ^{*}|\theta _{i})}}\right),}$${\displaystyle {\mathcal {L}}}$${\displaystyle P(\theta )}$${\displaystyle Q}$

• 遺伝的アルゴリズム
• 平均場粒子法
• メトロポリスライトトランスポート
• 複数回トライするメトロポリス
• 平行焼戻し
• シーケンシャルモンテカルロ
• シミュレーテッドアニーリング

• Bernd A. Berg .マルコフ連鎖モンテカルロシミュレーションとその統計分析シンガポール、World Scientific、2004年。
• チブ、シッダールタ、グリーンバーグ、エドワード（1995）。「メトロポリス・ヘイスティングスアルゴリズムの理解」アメリカ統計学者、49（4）、327-335。
• David DL Minh および Do Le Minh。「ヘイスティングス アルゴリズムの理解」Communications in Statistics - Simulation and Computation、44:2 332–349、2015 年
• ボルスタッド、ウィリアム M. (2010)計算ベイズ統計の理解、ジョン ワイリー アンド サンズ ISBN 0-470-04609-0