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

数学や物理において、メトロポリス・ヘイスティングス法(もしくは M-H アルゴリズム)(メトロポリス・ヘイスティングスほう、"Metropolis-Hastings algorithm")　はマルコフ連鎖モンテカルロ法の一つで、直接的に乱数の生成が難しい確率分布に対し、その確率分布に収束するマルコフ連鎖を生成する手法である。生成されたマルコフ連鎖は、確率分布の近似（ヒストグラム）などの期待値、すなわち積分の近似計算に用いられる。
歴史.
このアルゴリズムは1953年にボルツマン分布のための特殊形で発表したニコラス・メトロポリスらによって提案され. 、1970年に W.K. ヘイスティングスによって一般化された。ギブスサンプリング法は棄却されることのないM-H アルゴリズムと捉えることが出来て、調整パラメータが少ない点で構成が容易であるが、応用範囲は狭まる。
アルゴリズムの説明.
M-H アルゴリズムを用いると、確率分布formula_1 の確率密度関数もしくは確率関数（煩雑なので、以下では「確率密度関数」に統一する）に比例する関数が計算できるならば、
いかなる確率分布 formula_1 からでも、基本的には乱数列を得ることができる。
統計学や統計物理学では、しばしば興味のある確率密度関数の定数倍しかわからないことがある。定数倍部分がわからなくても乱数の生成が可能である点はM-H アルゴリズムの大きな長所である。
M-H アルゴリズムは乱数列を生成する。
乱数列の長さが長ければ長いほど、その乱数列は目標の分布P(x)を近似することになる。
乱数は反復計算によって生成される。次の乱数が従う確率分布は現在の乱数の値にのみ依存する。すなわち、乱数列はマルコフ連鎖である。
乱数列を生成する反復計算方法がアルゴリズムの重要な点である。反復計算は、次の状態となる候補を計算することと、それを採択、もしくは棄却する手続きで構成される。
採択とは候補となった状態が次の反復計算の際に使用されることで、棄却とは次の反復計算でも現在の状態が使用されることである。採択、もしくは棄却を決定する採択確率は、P(x)に関する現在の状態と次の状態の確率密度関数を比較して決定される。
以下では簡単のために、M-H アルゴリズムの最も基本的な例である、メトロポリスアルゴリズムを示す。
メトロポリスアルゴリズム
formula_3を目標分布formula_4の確率密度関数に比例する関数とする。
この反復計算により、乱数はある状態にとどまったり、動いたりを繰り返し、ランダムに状態空間を動きまわる。採択率は、現在の状態と次の状態の候補の確率密度関数を比べることで、生成された候補がどれだけ採択されうるかを表す。現在の状態よりも次の候補の確率密度関数が高ければ、必ず次の状態を採択する。しかし次の候補の確率密度関数が高くない場合、ある確率で移動せずに候補を棄却する。このため、高い確率密度関数の領域からは多くの乱数を、また低い確率密度関数の領域は少ない乱数を生成することになる。直感的にはこれがアルゴリズムの仕組みであり、目的の確率分布に従った乱数を近似的に生成する方法である。
棄却法など、確率分布から直接独立に乱数を生成する方法と異なり、M-H アルゴリズムなどのマルコフ連鎖モンテカルロ法はいくつかの短所がある。
確率分布の近似に使われる基本的な棄却法は、その棄却確率が次元数の関数として指数関数的に増加する、次元の呪いの影響を受ける。M-Hアルゴリズムなどマルコフ連鎖モンテカルロ法では、次元の影響が同程度には起きないことが多い。そのため、確率分布が定義されている状態空間の次元が高い場合、マルコフ連鎖モンテカルロ法を用いることは自然である。そのためマルコフ連鎖モンテカルロ法は、多くの分野で使用されている階層ベイズモデルや高次元な統計モデルの事後分布の近似方法としてよく選ばれている。
次元が高い場合には、個々の次元ごとに異なった振る舞いをすることや、遅い混合を避けるためにすべての次元に関して適切なジャンプの大きさを決定することが問題となるため、提案分布を適切に選択することが自体が困難である。
そのような状況でしばしば取られる代替案としてはギブスサンプリングがある。ギブスサンプリングは、すべての次元から一度にサンプルするのではなく、個々の次元に関してサンプリングを行う。
これは多くの典型な階層モデルにあるように、少数の変数が他の変数を条件付けている場合には有効な方法である。
個々の変数は他の変数に関して条件づけされて１度にサンプリングされる。
他には適応的棄却サンプリング、一次元M-H ステップ、スライスサンプリングが考えられる。
提案密度 formula_32 が formula_33 に一切依存しないことも可能であり、その場合はアルゴリズムは「独立型メトロポリス・ヘイスティングス法」という。
ふさわしい提案分布を持った独立型メトロポリス・ヘイスティングス法は高い精度が期待されるが、目標となる確率分布の事前の知識を必要とするし、次元の呪いを強く受ける危険がある。
定式化.
M-H アルゴリズムは目標確率分布formula_4に従ったサンプルの生成を行うことが目的である。
これを達成するために、漸近的に唯一の定常分布π(x)に収束するマルコフ連鎖を用いる。
ここでは簡単のため、離散の状態空間を考えることにする。マルコフ連鎖は、２つの状態間の遷移確率formula_35によって一意に定義される。
次の２つの条件が満たされるとき、マルコフ連鎖は定常分布に収束する。このとき、マルコフ連鎖はエルゴード性をもつという。
M-H アルゴリズムは遷移確率の構成により、上記の２つの条件を満たすようにマルコフ過程を設計することができる。
詳細釣り合い条件を確認しよう。formula_47として
formula_48
が成り立つ必要がある。これは、以下のように書き換えられる。
formula_49.
通常の手法として遷移確率を提案確率分布と採択確率分布に分解する。提案分布formula_50はformula_9が与えられたときの状態formula_40を提案する条件付き確率であり、採択確率formula_53はformula_9が与えられたときの状態formula_40を採択する条件付き確率である。
formula_56 
この関係を以前の式に代入して以下の式を得る。
formula_57 .
次のステップとして、この式を満たす採択率を選ぶことが必要である。よくある選択として、メトロポリス選択が知られ以下の式で得られる。この値はアルゴリズムの実装に必要な値である。
formula_58
この式が前の式を満たすことは、formula_59 /formula_60かformula_61/formula_59 の少なくとも片方 が1以上になることから確認できる。
また、これはformula_59とformula_61を一般性を失うことなく入れ替えることができるからである。
実装の観点からはMetropolis–Hastings アルゴリズムは以下のステップから成り立っている。
サンプルを適切に集めるためには、formula_71は提案分布や採択率とが別に決められ、ステップ４においてサンプルが相関していないことが必要である。マルコフ過程の自己相関時間の時間のオーダーによる。
一般的にこのパラメータの決定は簡単ではないことは重要な点である。問題に対して適切にパラメータは決定されるべきである。分布に関する知識が全くない場合には一様分布が提案分布として選ばれることもある。この場合、状態formula_9と状態formula_40はいつも相関しないためにformula_71の値はformula_76に設定される。
アルゴリズムの手順.
現時点のサンプル値は formula_77 であるとする。
メトロポリス・ヘイスティングス法では、次のサンプルformula_40は
確率密度関数 formula_79 に従い生成される。
そして以下の値を計算する。
ここで、
はサンプルの候補 formula_82 とその一つ前のサンプル formula_77 の尤度比であり、
は２つの提案分布(formula_77 から formula_82 へとその逆方向)の比率である。
提案分布が状態に関して対称の場合は formula_87 は 1 となる。
次に、以下のルールにもとづき新しい状態 formula_88 が選択される。
formula_89 の場合:
formula_91 の場合:
マルコフ連鎖はランダムな初期値 formula_96 から開始され、初期値が「忘れられる」まで、多くの試行を繰り返す。この間の標本は棄てられ、"burn-in"(機械や回路を通電した直後の安定しない状態からの比喩、初期通電)とよばれる。
burn-in'後の値 formula_9 の集合は、分布 formula_4 からのサンプルを表す。
M-Hアルゴリズムは提案分布 formula_99の形が、直接サンプリングが困難である目標分布 formula_1 の形と類似している場合、つまり formula_101 のときにうまくアルゴリズムが動作する。
しかし、多くの場合目標分布の形は未知である。
ガウス分布の提案分布 formula_102 が用いられる場合は、分散パラメータの formula_103 が burn-in 期間のうちに調整される必要がある。これは通常、"採択率"を計算することで行われる。採択率とは formula_104個のサンプルのうち採択されたサンプルの割合である。要求される採択率は目標分布によって異なる。理論的には、一次元ガウス分布を目標分布とすると理想的な採択率は約50% であり、N次元ガウス分布の目標分布では約 23% であることが知られている。
formula_103 が小さすぎるとサンプリング列は低速混合をおこす。
つまり採択率が高くなり標本空間の移動が遅くなる。
そのため formula_1 への収束が遅くなる。
逆に formula_103 が大きすぎると、
採択率が低すぎサンプルが確率密度の低い所に移動してしまい formula_108 が非常に小さくなってしまう。
このため formula_1 への収束が遅くなる。
したがって、提案分布のパラメータを調整し"採択率"を適切にする必要がある。
関連記事.
M-H アルゴリズムの一例として
MCMCではない方法
モンテカルロEMアルゴリズム：EステップをサンプリングにしたEMアルゴリズム。