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

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

メトロポリス・ヘイスティングス法をはじめとするMCMCアルゴリズムは、一般的に多次元分布からのサンプリング、特に次元数が多い場合に用いられます。一方、一次元分布の場合、分布から直接独立したサンプルを取得できる他の手法（例えば、適応型棄却サンプリング）が通常存在し、これらの手法はMCMC法に内在する自己相関サンプルの問題から解放されます。

このアルゴリズムは、1953 年にArianna W. Rosenbluth、Marshall Rosenbluth、Augusta H. Teller、Edward Tellerと共に「高速計算機による状態方程式の計算」というタイトルの論文の最初の共著者であるNicholas Metropolisにちなんで名付けられました。長年にわたり、このアルゴリズムは単にメトロポリス アルゴリズムとして知られていました。^{[ 1 ] [ 2 ]}この論文では、対称的な提案分布の場合にこのアルゴリズムが提案されましたが、1970 年にWK Hastings がそれをより一般的なケースに拡張しました。^{[ 3 ]}一般化された方法は最終的に両方の名前で識別されるようになりましたが、「メトロポリス-ヘイスティングス アルゴリズム」という用語が最初に使用された時期は不明です。

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

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

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

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

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

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

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

させて${\displaystyle f(x)}$目的の確率密度関数に比例する関数である${\displaystyle P(x)}$（別名：目標分布）^{[ a ]}

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})}$。
   • 合格率を計算する${\displaystyle \alpha =f(x')/f(x_{t})}$これは、候補者を受け入れるか拒否するかを決定するために使用されます。^{[ b ]} f はPの密度に比例するため、次のようになります。${\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次元メトロポリス-ヘイスティングスステップ、またはスライスサンプリングなどがあります。

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

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

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. 定常分布の一意性：定常分布${\displaystyle \pi (x)}$一意でなければならない。これはマルコフ過程のエルゴード性によって保証されており、すべての状態は(1)非周期的であること（システムが一定の間隔で同じ状態に戻らないこと）、および(2)正の再帰性を持つこと（同じ状態に戻るまでの期待ステップ数が有限であること）が必要である。

メトロポリス・ヘイスティングスアルゴリズムは、上記の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)}}.}$

導出の次のステップは、上記の条件を満たす受容率を選択することです。よく用いられる選択肢の一つは、メトロポリス選択です。

${\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}$。

指定された条件が満たされる場合、保存された状態の経験的分布${\displaystyle x_{0},\ldots ,x_{T}}$近づく${\displaystyle P(x)}$反復回数（${\displaystyle T}$効果的に推定するために必要なもの${\displaystyle P(x)}$関係性を含む要因の数に依存します${\displaystyle P(x)}$提案分布と推定の望ましい精度。^{[ 13 ]}離散状態空間上の分布の場合、マルコフ過程の自己相関時間のオーダーでなければならない。 ^{[ 14 ]メトロポリス-ヘイスティングスの収束理論の分かりやすい説明は[ 15 ]}に示されている。

一般的な問題では、どの分布が${\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 )}}$

そして、推定すると${\displaystyle P(E)}$指標関数の期待値を推定することで達成できる。${\displaystyle A_{E}(x)\equiv \mathbf {1} _{E}(x)}$、これは 1 のとき${\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 方向の提案密度の比率です (${\displaystyle x_{t}}$に${\displaystyle x'}$（逆もまた同様です）。提案密度が対称であれば、これは 1 に等しくなります。すると新しい状態は${\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)}$。

提案密度が目標分布の形状と一致する場合に、このアルゴリズムは最も効果的に機能します。${\displaystyle P(x)}$直接サンプリングが難しい、つまり${\displaystyle g(x'\mid x_{t})\approx P(x')}$ガウス型提案密度の場合${\displaystyle g}$分散パラメータが使用される${\displaystyle \sigma ^{2}}$バーンイン期間中に調整する必要があります。これは通常、受理率を計算することによって行われます。受理率とは、最後の期間のウィンドウ内で提案されたサンプルのうち受理されたサンプルの割合です。${\displaystyle N}$サンプル。望ましい受入率は目標分布に依存しますが、理論的には、一次元ガウス分布の理想的な受入率は約 50% であり、${\displaystyle N}$次元ガウス目標分布。^{[ 16 ]}これらのガイドラインは、ベルンシュタイン-フォン・ミーゼスの定理を用いて確立できる多変量正規分布に従うことが多い、十分に規則的なベイズ事後分布からサンプリングする場合にうまく機能します。^{[ 17 ]}

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

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}$（条件付き）提案確率。

• 遺伝的アルゴリズム
• 平均場粒子法
• 大都市軽便輸送
• 複数回挑戦できるメトロポリス
• 平行焼戻し
• 逐次モンテカルロ法
• シミュレーテッドアニーリング

• Berg, Bernd A (2004年10月).マルコフ連鎖モンテカルロシミュレーションとその統計分析：WebベースのFortranコード付き. World Scientific . doi : 10.1142/5602 . ISBN 978-981-238-935-0。
• Chib, Siddhartha; Greenberg, Edward (1995年11月) 「メトロポリス・ヘイスティングスアルゴリズムの理解」『アメリカ統計学者』49 (4):327. doi : 10.2307/2684568 .
• Minh, David DL; Minh, Do Le (Paul) (2015-02-07). "Hastingsアルゴリズムの理解" . Communications in Statistics - Simulation and Computation . 44 (2): 332– 349. arXiv : 1408.4438 . doi : 10.1080/03610918.2013.777455 . ISSN 0361-0918 . 
• ボルスタッド、ウィリアム・M. (2010).計算ベイズ統計学の理解. 計算統計学シリーズ. ホーボーケン、ニュージャージー州: ワイリー. ISBN 978-0-470-04609-8。