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

数学や物理において、メトロポリス・ヘイスティングス法(もしくは M-H アルゴリズム)(メトロポリス・ヘイスティングスほう、Metropolis-Hastings algorithm)　はマルコフ連鎖モンテカルロ法の一つで、直接的に乱数の生成が難しい確率分布に対し、その確率分布に収束するマルコフ連鎖を生成する手法である。生成されたマルコフ連鎖は、確率分布の近似（ヒストグラム）などの期待値、すなわち積分の近似計算に用いられる。

歴史

このアルゴリズムは1953年にボルツマン分布のための特殊形で発表したニコラス・メトロポリスらによって提案され.^[1] 、1970年に W.K. ヘイスティングスによって一般化された^[2]。ギブスサンプリング法は棄却されることのないM-H アルゴリズムと捉えることが出来て、調整パラメータが少ない点で構成が容易であるが、応用範囲は狭まる。

アルゴリズムの説明

M-H アルゴリズムを用いると、確率分布${\displaystyle \displaystyle P(x)}$ の確率密度関数もしくは確率関数（煩雑なので、以下では「確率密度関数」に統一する）に比例する関数が計算できるならば、 いかなる確率分布 ${\displaystyle \displaystyle P(x)}$ からでも、基本的には乱数列を得ることができる。

統計学や統計物理学では、しばしば興味のある確率密度関数の定数倍しかわからないことがある。定数倍部分がわからなくても乱数の生成が可能である点はM-H アルゴリズムの大きな長所である。

M-H アルゴリズムは乱数列を生成する。 乱数列の長さが長ければ長いほど、その乱数列は目標の分布P(x)を近似することになる。 乱数は反復計算によって生成される。次の乱数が従う確率分布は現在の乱数の値にのみ依存する。すなわち、乱数列はマルコフ連鎖である。

乱数列を生成する反復計算方法がアルゴリズムの重要な点である。反復計算は、次の状態となる候補を計算することと、それを採択、もしくは棄却する手続きで構成される。 採択とは候補となった状態が次の反復計算の際に使用されることで、棄却とは次の反復計算でも現在の状態が使用されることである。採択、もしくは棄却を決定する採択確率は、P(x)に関する現在の状態と次の状態の確率密度関数を比較して決定される。

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

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

${\displaystyle f(x)}$を目標分布${\displaystyle P(x)}$の確率密度関数に比例する関数とする。

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

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

棄却法など、確率分布から直接独立に乱数を生成する方法と異なり、M-H アルゴリズムなどのマルコフ連鎖モンテカルロ法はいくつかの短所がある。

• 乱数列が相関していること。長い乱数列を生成しても、近接した乱数は相関をもち、確率分布を正しく反映したものではない。相関が高ければ、目標となる確率分布を近似するのに、多くの反復計算が必要になる。相関を小さくする工夫はいくつかある。ひとつは、あらかじめ決めた整数${\displaystyle n}$に対し、${\displaystyle n}$個飛ばしで乱数を集める方法である。整数${\displaystyle n}$の値は、しばしば乱数列を用いて計算した自己相関関数をもとに決められる。しかし整数${\displaystyle n}$の値が大きすぎると、その分反復計算が増えてしまう。もうひとつは、次の状態の候補を現在時点より、適度に遠くへ提案する方法である。しかしあまり遠くへ提案しすぎれば棄却確率が増加してかえって相関が小さくなる。
• 反復計算で生成されたマルコフ連鎖は、ゆるい十分条件のもとで、目標の分布に収束されることが保証されている。初期値が確率密度関数の小さい領域にあると、目標の分布に近づくのが遅れ、目標の分布の近似に大きなバイアスを生む危険がある。その対策として、反復計算の最初の方の部分を切り捨てる、burn-inがしばしば有効である。

確率分布の近似に使われる基本的な棄却法は、その棄却確率が次元数の関数として指数関数的に増加する、次元の呪いの影響を受ける。M-Hアルゴリズムなどマルコフ連鎖モンテカルロ法では、次元の影響が同程度には起きないことが多い。そのため、確率分布が定義されている状態空間の次元が高い場合、マルコフ連鎖モンテカルロ法を用いることは自然である。そのためマルコフ連鎖モンテカルロ法は、多くの分野で使用されている階層ベイズモデルや高次元な統計モデルの事後分布の近似方法としてよく選ばれている。

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

提案密度 ${\displaystyle \displaystyle Q(x'|x_{t})}$ が ${\displaystyle \displaystyle x_{t}}$ に一切依存しないことも可能であり、その場合はアルゴリズムは「独立型メトロポリス・ヘイスティングス法」という。 ふさわしい提案分布を持った独立型メトロポリス・ヘイスティングス法は高い精度が期待されるが、目標となる確率分布の事前の知識を必要とするし、次元の呪いを強く受ける危険がある。

定式化

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

ここでは簡単のため、離散の状態空間を考えることにする。マルコフ連鎖は、２つの状態間の遷移確率${\displaystyle P(x'|x)}$によって一意に定義される。 次の２つの条件が満たされるとき、マルコフ連鎖は定常分布に収束する^[3]。このとき、マルコフ連鎖はエルゴード性をもつという。

1. 定常分布の存在：定常である確率分布${\displaystyle \pi (x)}$が存在しなければならない。一つの十分条件として、詳細釣り合い条件がある。詳細釣り合い条件とは、状態${\displaystyle x}$が${\displaystyle \pi (x)}$からの乱数であるとき、状態${\displaystyle x}$から状態${\displaystyle x'}$への遷移確率が状態${\displaystyle x'}$から状態${\displaystyle x}$への遷移確率と等しいこと、つまり、${\displaystyle \pi (x)P(x'|x)=\pi (x')P(x|x')}$となることである。
2. 定常分布の一意性: 定常分布${\displaystyle \pi (x)}$は一意でなければならない。十分条件の一つは、${\displaystyle P(x'|x)}$がすべての${\displaystyle x,x'}$について正になることである^[4]。

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

詳細釣り合い条件を確認しよう。${\displaystyle \pi (x)=P(x)}$として

${\displaystyle P(x)P(x'|x)=P(x')P(x|x')}$

が成り立つ必要がある。これは、以下のように書き換えられる。

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

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

${\displaystyle P(x'|x)=Q(x'|x)A(x',x)}$

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

${\displaystyle {\frac {A(x',x)}{A(x,x')}}={\frac {P(x')}{P(x)}}{\frac {Q(x|x')}{Q(x'|x)}}}$ .

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

${\displaystyle A(x',x)=\min \left(1,{\frac {P(x')}{P(x)}}{\frac {Q(x|x')}{Q(x'|x)}}\right)}$

この式が前の式を満たすことは、${\displaystyle A(x',x)}$ /${\displaystyle A(x,x')}$か${\displaystyle A(x,x')}$/${\displaystyle A(x',x)}$ の少なくとも片方 が1以上になることから確認できる。

また、これは${\displaystyle A(x',x)}$と${\displaystyle A(x,x')}$を一般性を失うことなく入れ替えることができるからである。

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

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

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

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

アルゴリズムの手順

現時点のサンプル値は ${\displaystyle x^{t}\,}$ であるとする。 メトロポリス・ヘイスティングス法では、次のサンプル${\displaystyle x'}$は 確率密度関数 ${\displaystyle Q(x'|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 {Q(x^{t}|x')}{Q(x'|x^{t})}}}$

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

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

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

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

${\displaystyle \displaystyle a<1}$ の場合:

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

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

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

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

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

${\displaystyle \displaystyle \sigma ^{2}}$ が小さすぎるとサンプリング列は低速混合をおこす。 つまり採択率が高くなり標本空間の移動が遅くなる。 そのため ${\displaystyle \displaystyle P(x)}$ への収束が遅くなる。 逆に ${\displaystyle \displaystyle \sigma ^{2}}$ が大きすぎると、 採択率が低すぎサンプルが確率密度の低い所に移動してしまい ${\displaystyle \displaystyle a_{1}}$ が非常に小さくなってしまう。 このため ${\displaystyle \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

脚注

外部リンク

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