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

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

歴史

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

アルゴリズムの説明

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

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

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

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

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

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

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

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

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

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

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

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

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

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

定式化

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

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

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

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

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

$P(x)P(x'|x)=P(x')P(x|x')$

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

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

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

$P(x'|x)=Q(x'|x)A(x',x)$

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

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

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

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

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

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

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

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

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

一般的にこのパラメータの決定は簡単ではないことは重要な点である。問題に対して適切にパラメータは決定されるべきである。分布に関する知識が全くない場合には一様分布が提案分布として選ばれることもある。この場合、状態 $x$ と状態 $x'$ はいつも相関しないために $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)$ への収束が遅くなる。したがって、提案分布のパラメータを調整し採択率を適切にする必要がある。

脚注

[脚注の使い方]

^ Metropolis, N.; Rosenbluth, A.W.; Rosenbluth, M.N.; Teller, A.H.; Teller, E. (1953). “Equation of State Calculations by Fast Computing Machines”. Journal of Chemical Physics 21 (6): 1087–1092. Bibcode: 1953JChPh..21.1087M. doi:10.1063/1.1699114.
^ Hastings, W.K. (1970). “Monte Carlo Sampling Methods Using Markov Chains and Their Applications”. Biometrika 57 (1): 97–109. Bibcode: 1970Bimka..57...97H. doi:10.1093/biomet/57.1.97. JSTOR 2334940. Zbl 0219.65008.
^ ^a ^b Robert, Christian; Casella, George (2004). Monte Carlo Statistical Methods. Springer. ISBN 0387212396
^ Kulik, Alexei (2017). Ergodic Behavior of Markov Processes: With Applications to Limit Theorems. De Gruyter. ISBN 978-3110458701
^ Newman, M. E. J.; Barkema, G. T. (1999). Monte Carlo Methods in Statistical Physics. USA: Oxford University Press. ISBN 0198517971

外部リンク

[metropolis-1] Metropolis, N.; Rosenbluth, A.W.; Rosenbluth, M.N.; Teller, A.H.; Teller, E. (1953). “Equation of State Calculations by Fast Computing Machines”. Journal of Chemical Physics 21 (6): 1087–1092. Bibcode: 1953JChPh..21.1087M. doi:10.1063/1.1699114.

[Hastings-2] Hastings, W.K. (1970). “Monte Carlo Sampling Methods Using Markov Chains and Their Applications”. Biometrika 57 (1): 97–109. Bibcode: 1970Bimka..57...97H. doi:10.1093/biomet/57.1.97. JSTOR 2334940. Zbl 0219.65008.

[Roberts_Casella-3] Robert, Christian; Casella, George (2004). Monte Carlo Statistical Methods. Springer. ISBN 0387212396

[Kulik-4] Kulik, Alexei (2017). Ergodic Behavior of Markov Processes: With Applications to Limit Theorems. De Gruyter. ISBN 978-3110458701

[Newman_Barkema-5] Newman, M. E. J.; Barkema, G. T. (1999). Monte Carlo Methods in Statistical Physics. USA: Oxford University Press. ISBN 0198517971

[1]

[2]

[3]

[4]

[5]

表話編歴確率分布
離散単変量で有限台	ベンフォードベルヌーイベータ二項（英語版）二項 categorical（英語版）超幾何ポワソン二項ラーデマッハ（英語版）離散一様ジップジップ–マンデルブロー（英語版）
離散単変量で無限台	ベータ負二項（英語版）ボレル（英語版）コンウェイ–マクスウェル–ポワソン（英語版）離散位相型（英語版）ドラポルト（英語版）拡張負二項（英語版）ガウス–クズミン幾何対数（英語版）負の二項放物フラクタル（英語版）ポワソンスケラム（英語版）ユール–サイモン（英語版）ゼータ（英語版）
連続単変量で有界区間に台を持つ	逆正弦（英語版） ARGUS（英語版）バルディング–ニコルス（英語版）ベイツ（英語版）ベータ beta rectangular（英語版）アーウィン–ホール（英語版）クマラスワミー（英語版）ロジット-正規（英語版）非中心ベータ（英語版） raised cosine（英語版） reciprocal（英語版）三角 U-quadratic（英語版）一様ウィグナー半円
連続単変量で半無限区間に台を持つ	ベニーニ（英語版）ベンクタンダー第一種（英語版）ベンクタンダー第二種（英語版）第2種ベータ Burr（英語版）カイ二乗カイ（英語版） Dagum（英語版）デービス（英語版）指数-対数（英語版）アーラン指数 F folded normal（英語版） Flory–Schulz（英語版）フレシェガンマ gamma/Gompertz（英語版）一般逆ガウス（英語版） Gompertz（英語版） half-logistic（英語版） half-normal（英語版） Hotelling's T-squared（英語版）超アーラン（英語版）超指数（英語版） hypoexponential（英語版）逆カイ二乗（英語版） scaled inverse chi-squared（英語版）逆ガウス逆ガンマコルモゴロフレヴィ対数コーシー対数ラプラス（英語版）対数ロジスティック（英語版）対数正規ロマックス（英語版）行列指数（英語版）マクスウェル–ボルツマンマクスウェル–ユットナー（英語版）ミッタク-レフラー（英語版）仲上（英語版）非心カイ二乗パレート位相型（英語版） poly-Weibull（英語版）レイリー relativistic Breit–Wigner（英語版）ライス（英語版） shifted Gompertz（英語版）切断正規タイプ2ガンベル（英語版）ワイブル離散ワイブル（英語版）ウィルクスのラムダ（英語版）
連続単変量で実数直線全体に台を持つ	コーシー（ローレンツ、ブライト・ウィグナー）指数冪（英語版）フィッシャーの z（英語版）ガウスの q（英語版）一般正規（英語版）一般化双曲型幾何安定（英語版）ガンベルホルツマルク（英語版）双曲線正割ジョンソンの S_U（英語版）ランダウラプラス非対称ラプラス（英語版）ロジスティック非心 t 正規 (ガウス) 正規逆ガウス（英語版）歪正規（英語版）スラッシュ安定スチューデントの t タイプ1ガンベル（英語版）トレイシー–ウィダム（英語版）分散ガンマ（英語版）フォークト
連続単変量でタイプの変わる台を持つ	一般極値一般パレート（英語版）マルチェンコ–パストゥール（英語版） q-指数（英語版） q-ガウス q-ワイブル（英語版） shifted log-logistic（英語版）トゥーキーのラムダ（英語版）
混連続-離散単変量	rectified Gaussian（英語版）
多変量 (結合)	【離散】エウェンズ（英語版）多項ディリクレ多項（英語版）負多項（英語版）【連続】ディリクレ一般ディリクレ（英語版）多変量正規多変量安定（英語版）多変量 t（英語版）正規逆ガンマ（英語版）正規ガンマ（英語版）【行列値】逆行列ガンマ（英語版）逆ウィッシャート（英語版）行列正規（英語版）行列 t（英語版）行列ガンマ（英語版）正規逆ウィッシャート（英語版）正規ウィッシャート（英語版）ウィッシャート
方向	【単変量 (円周) 方向】円周一様（英語版）単変数フォン・ミーゼス wrapped 正規（英語版） wrapped コーシー（英語版） wrapped 指数（英語版） wrapped 非対称ラプラス（英語版） wrapped レヴィ（英語版）【二変量 (球面)】ケント（英語版）【二変量 (トロイダル)】二変数フォン・ミーゼス（英語版）【多変量】フォン・ミーゼス–フィッシャー（英語版）ビンガム（英語版）
退化と特異	【退化】ディラックのデルタ関数【特異】カントール
族	円周（英語版）混合ポワソン（英語版）楕円（英語版）指数自然指数（英語版）位置尺度（英語版）最大エントロピー（英語版）混合（英語版）ピアソン（英語版）トウィーディ（英語版） wrapped（英語版）
サンプリング法（英語版）	逆関数サンプリング法マルコフ連鎖モンテカルロ法（メトロポリス・ヘイスティングス法・ギブスサンプリング・スライスサンプリング）粒子フィルタボックス＝ミュラー法棄却サンプリング（英語版）ジッグラト法（英語版）マルサグリア法（英語版）
一覧（英語版）カテゴリ

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

歴史

アルゴリズムの説明

定式化

アルゴリズムの手順

関連記事

関連書籍

脚注

外部リンク