提案分布 Q はランダムウォーク の粒子が次に移動する候補点を提案する。 数学 や物理 において、メトロポリス・ヘイスティングス法 (もしくは M-H アルゴリズム )(メトロポリス・ヘイスティングスほう、Metropolis-Hastings algorithm ) はマルコフ連鎖モンテカルロ法 の一つで、直接的に乱数の生成が難しい確率分布 に対し、その確率分布に収束するマルコフ連鎖 を生成する手法である。生成されたマルコフ連鎖は、確率分布の近似（ヒストグラム ）などの期待値、すなわち積分の近似計算に用いられる。
このアルゴリズムは1953年にボルツマン分布 のための特殊形で発表したニコラス・メトロポリス らによって提案され.[ 1] [ 2] ギブスサンプリング 法は棄却されることのないM-H アルゴリズム と捉えることが出来て、調整パラメータが少ない点で構成が容易であるが、応用範囲は狭まる。
M-H アルゴリズム を用いると、確率分布
P
(
x
)
{\displaystyle \displaystyle P(x)}
確率分布
P
(
x
)
{\displaystyle \displaystyle P(x)}
統計学や統計物理学では、しばしば興味のある確率密度関数の定数倍しかわからないことがある。定数倍部分がわからなくても乱数の生成が可能である点はM-H アルゴリズム の大きな長所である。
M-H アルゴリズムは乱数列を生成する。
乱数列の長さが長ければ長いほど、その乱数列は目標の分布P(x)を近似することになる。
乱数は反復計算によって生成される。次の乱数が従う確率分布は現在の乱数の値にのみ依存する。すなわち、乱数列はマルコフ連鎖である。
乱数列を生成する反復計算方法がアルゴリズムの重要な点である。反復計算は、次の状態となる候補を計算することと、それを採択、もしくは棄却する手続きで構成される。
採択とは候補となった状態が次の反復計算の際に使用されることで、棄却とは次の反復計算でも現在の状態が使用されることである。採択、もしくは棄却を決定する採択確率は、P(x)に関する現在の状態と次の状態の確率密度関数を比較して決定される。
以下では簡単のために、M-H アルゴリズムの最も基本的な例である、メトロポリスアルゴリズムを示す。
メトロポリスアルゴリズム
f
(
x
)
{\displaystyle f(x)}
P
(
x
)
{\displaystyle P(x)}
初期化: 初期値
x
0
{\displaystyle x_{0}}
Q
(
x
|
y
)
{\displaystyle Q(x|y)}
関数
Q
(
x
|
y
)
{\displaystyle Q(x|y)}
y
{\displaystyle y}
x
{\displaystyle x}
Q
{\displaystyle Q}
Q
(
x
|
y
)
=
Q
(
y
|
x
)
{\displaystyle Q(x|y)=Q(y|x)}
Q
(
x
|
y
)
{\displaystyle Q(x|y)}
y
{\displaystyle y}
ガウス分布 があり、この場合、
y
{\displaystyle y}
Q
{\displaystyle Q}
提案分布 と呼ばれる。
t
{\displaystyle t}
確率分布
Q
(
x
∗
|
x
t
)
{\displaystyle Q(x^{*}|x_{t})}
x
∗
{\displaystyle x^{*}}
採択率
α
=
f
(
x
∗
)
f
(
x
t
)
{\displaystyle \alpha ={\frac {f(x^{*})}{f(x_{t})}}}
P
{\displaystyle P}
f
{\displaystyle f}
α
=
f
(
x
∗
)
f
(
x
t
)
=
P
(
x
∗
)
P
(
x
t
)
{\displaystyle \alpha ={\frac {f(x^{*})}{f(x_{t})}}={\frac {P(x^{*})}{P(x_{t})}}}
α
≥
1
{\displaystyle \alpha \geq 1}
x
∗
{\displaystyle x^{*}}
x
t
+
1
=
x
∗
{\displaystyle x_{t+1}=x^{*}}
α
{\displaystyle \alpha }
x
t
+
1
=
x
t
{\displaystyle x_{t+1}=x_{t}}
この反復計算により、乱数はある状態にとどまったり、動いたりを繰り返し、ランダムに状態空間を動きまわる。採択率は、現在の状態と次の状態の候補の確率密度関数を比べることで、生成された候補がどれだけ採択されうるかを表す。現在の状態よりも次の候補の確率密度関数が高ければ、必ず次の状態を採択する。しかし次の候補の確率密度関数が高くない場合、ある確率で移動せずに候補を棄却する。このため、高い確率密度関数の領域からは多くの乱数を、また低い確率密度関数の領域は少ない乱数を生成することになる。直感的にはこれがアルゴリズムの仕組みであり、目的の確率分布に従った乱数を近似的に生成する方法である。
棄却法など、確率分布から直接独立に乱数を生成する方法と異なり、M-H アルゴリズムなどのマルコフ連鎖モンテカルロ法はいくつかの短所がある。
乱数列が相関していること。長い乱数列を生成しても、近接した乱数は相関をもち、確率分布を正しく反映したものではない。相関が高ければ、目標となる確率分布を近似するのに、多くの反復計算が必要になる。相関を小さくする工夫はいくつかある。ひとつは、あらかじめ決めた整数
n
{\displaystyle n}
n
{\displaystyle n}
n
{\displaystyle n}
n
{\displaystyle n}
反復計算で生成されたマルコフ連鎖は、ゆるい十分条件のもとで、目標の分布に収束されることが保証されている。初期値が確率密度関数の小さい領域にあると、目標の分布に近づくのが遅れ、目標の分布の近似に大きなバイアスを生む危険がある。その対策として、反復計算の最初の方の部分を切り捨てる、burn-inがしばしば有効である。 確率分布の近似に使われる基本的な棄却法は、その棄却確率が次元数の関数として指数関数的に増加する、次元の呪い の影響を受ける。M-Hアルゴリズムなどマルコフ連鎖モンテカルロ法では、次元の影響が同程度には起きないことが多い。そのため、確率分布が定義されている状態空間の次元が高い場合、マルコフ連鎖モンテカルロ法を用いることは自然である。そのためマルコフ連鎖モンテカルロ法は、多くの分野で使用されている階層ベイズモデルや高次元な統計モデルの事後分布の近似方法としてよく選ばれている。
次元が高い場合には、個々の次元ごとに異なった振る舞いをすることや、遅い混合を避けるためにすべての次元に関して適切なジャンプの大きさを決定することが問題となるため、提案分布を適切に選択することが自体が困難である。
そのような状況でしばしば取られる代替案としてはギブスサンプリング がある。ギブスサンプリングは、すべての次元から一度にサンプルするのではなく、個々の次元に関してサンプリングを行う。
これは多くの典型な階層モデルにあるように、少数の変数が他の変数を条件付けている場合には有効な方法である。
個々の変数は他の変数に関して条件づけされて１度にサンプリングされる。
他には適応的棄却サンプリング 、一次元M-H ステップ、スライスサンプリング が考えられる。
提案密度
Q
(
x
′
|
x
t
)
{\displaystyle \displaystyle Q(x'|x_{t})}
x
t
{\displaystyle \displaystyle x_{t}}
M-H アルゴリズムは目標確率分布
P
(
x
)
{\displaystyle P(x)}
[ 3]
ここでは簡単のため、離散の状態空間を考えることにする。マルコフ連鎖は、２つの状態間の遷移確率
P
(
x
′
|
x
)
{\displaystyle P(x'|x)}
[ 3]
定常分布の存在 ：定常である確率分布
π
(
x
)
{\displaystyle \pi (x)}
x
{\displaystyle x}
π
(
x
)
{\displaystyle \pi (x)}
x
{\displaystyle x}
x
′
{\displaystyle x'}
x
′
{\displaystyle x'}
x
{\displaystyle x}
π
(
x
)
P
(
x
′
|
x
)
=
π
(
x
′
)
P
(
x
|
x
′
)
{\displaystyle \pi (x)P(x'|x)=\pi (x')P(x|x')}
定常分布の一意性 : 定常分布
π
(
x
)
{\displaystyle \pi (x)}
P
(
x
′
|
x
)
{\displaystyle P(x'|x)}
x
,
x
′
{\displaystyle x,x'}
[ 4] M-H アルゴリズムは遷移確率の構成により、上記の２つの条件を満たすようにマルコフ過程を設計することができる。
詳細釣り合い条件を確認しよう。
π
(
x
)
=
P
(
x
)
{\displaystyle \pi (x)=P(x)}
P
(
x
)
P
(
x
′
|
x
)
=
P
(
x
′
)
P
(
x
|
x
′
)
{\displaystyle P(x)P(x'|x)=P(x')P(x|x')}
が成り立つ必要がある。これは、以下のように書き換えられる。
P
(
x
′
|
x
)
P
(
x
|
x
′
)
=
P
(
x
′
)
P
(
x
)
{\displaystyle {\frac {P(x'|x)}{P(x|x')}}={\frac {P(x')}{P(x)}}}
通常の手法として遷移確率を提案確率分布と採択確率分布に分解する。提案分布
Q
(
x
′
|
x
)
{\displaystyle \displaystyle Q(x'|x)}
x
{\displaystyle x}
x
′
{\displaystyle x'}
採択確率
A
(
x
′
,
x
)
{\displaystyle \displaystyle A(x',x)}
x
{\displaystyle x}
x
′
{\displaystyle x'}
P
(
x
′
|
x
)
=
Q
(
x
′
|
x
)
A
(
x
′
,
x
)
{\displaystyle P(x'|x)=Q(x'|x)A(x',x)}
この関係を以前の式に代入して以下の式を得る。
A
(
x
′
,
x
)
A
(
x
,
x
′
)
=
P
(
x
′
)
P
(
x
)
Q
(
x
|
x
′
)
Q
(
x
′
|
x
)
{\displaystyle {\frac {A(x',x)}{A(x,x')}}={\frac {P(x')}{P(x)}}{\frac {Q(x|x')}{Q(x'|x)}}}
次のステップとして、この式を満たす採択率を選ぶことが必要である。よくある選択として、メトロポリス選択が知られ以下の式で得られる。この値はアルゴリズムの実装に必要な値である。
A
(
x
′
,
x
)
=
min
(
1
,
P
(
x
′
)
P
(
x
)
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)}
この式が前の式を満たすことは、
A
(
x
′
,
x
)
{\displaystyle A(x',x)}
A
(
x
,
x
′
)
{\displaystyle A(x,x')}
A
(
x
,
x
′
)
{\displaystyle A(x,x')}
A
(
x
′
,
x
)
{\displaystyle A(x',x)}
また、これは
A
(
x
′
,
x
)
{\displaystyle A(x',x)}
A
(
x
,
x
′
)
{\displaystyle A(x,x')}
実装の観点からはMetropolis–Hastings アルゴリズムは以下のステップから成り立っている。
初期化： ランダムに
x
{\displaystyle x}
Q
(
x
′
|
x
)
{\displaystyle \displaystyle Q(x'|x)}
x
′
{\displaystyle x'}
A
(
x
′
,
x
)
{\displaystyle \displaystyle A(x',x)}
x
′
{\displaystyle x'}
x
′
=
x
{\displaystyle x'=x}
T
{\displaystyle T}
値を保存する。２に戻る。 サンプルを適切に集めるためには、
T
{\displaystyle T}
[ 5]
一般的にこのパラメータの決定は簡単ではないことは重要な点である。問題に対して適切にパラメータは決定されるべきである。分布に関する知識が全くない場合には一様分布が提案分布として選ばれることもある。この場合、状態
x
{\displaystyle x}
x
′
{\displaystyle x'}
T
{\displaystyle T}
1
{\displaystyle 1}
現時点のサンプル値は
x
t
{\displaystyle x^{t}\,}
x
′
{\displaystyle x'}
Q
(
x
′
|
x
t
)
{\displaystyle Q(x'|x^{t})\,}
a
=
a
1
a
2
{\displaystyle a=a_{1}a_{2}\,}
ここで、
a
1
=
P
(
x
′
)
P
(
x
t
)
{\displaystyle a_{1}={\frac {P(x')}{P(x^{t})}}\,\!}
はサンプルの候補
x
′
{\displaystyle x'\,}
x
t
{\displaystyle x^{t}\,}
a
2
=
Q
(
x
t
|
x
′
)
Q
(
x
′
|
x
t
)
{\displaystyle a_{2}={\frac {Q(x^{t}|x')}{Q(x'|x^{t})}}}
は２つの提案分布(
x
t
{\displaystyle x^{t}\,}
x
′
{\displaystyle x'\,}
a
2
{\displaystyle \displaystyle a_{2}}
次に、以下のルールにもとづき新しい状態
x
t
+
1
{\displaystyle x^{t+1}\,}
a
≥
1
{\displaystyle \displaystyle a\geq 1}
x
t
+
1
=
x
′
{\displaystyle \displaystyle x^{t+1}=x'}
a
<
1
{\displaystyle \displaystyle a<1}
a
{\displaystyle \displaystyle a}
x
t
+
1
=
x
′
{\displaystyle \displaystyle x^{t+1}=x'}
1
−
a
{\displaystyle \displaystyle 1-a}
x
t
+
1
=
x
t
{\displaystyle \displaystyle x^{t+1}=x^{t}}
マルコフ連鎖はランダムな初期値
x
0
{\displaystyle \displaystyle x^{0}}
burn-in (機械や回路を通電した直後の安定しない状態からの比喩、初期通電)とよばれる。
burn-in'後の値
x
{\displaystyle x}
P
(
x
)
{\displaystyle P(x)}
サンプル を表す。
Q
(
x
′
;
x
t
)
{\displaystyle Q(x';x^{t})}
P
(
x
)
{\displaystyle \displaystyle P(x)}
Q
(
x
′
|
x
t
)
≈
P
(
x
′
)
{\displaystyle Q(x'|x^{t})\approx P(x')\,\!}
ガウス分布 の提案分布
Q
{\displaystyle \displaystyle Q}
σ
2
{\displaystyle \displaystyle \sigma ^{2}}
採択率 を計算することで行われる。採択率とは
N
{\displaystyle \displaystyle N}
σ
2
{\displaystyle \displaystyle \sigma ^{2}}
低速混合 をおこす。
つまり採択率が高くなり標本空間の移動が遅くなる。
そのため
P
(
x
)
{\displaystyle \displaystyle P(x)}
σ
2
{\displaystyle \displaystyle \sigma ^{2}}
a
1
{\displaystyle \displaystyle a_{1}}
P
(
x
)
{\displaystyle \displaystyle P(x)}
採択率 を適切にする必要がある。
M-H アルゴリズム の一例として
MCMCではない方法
モンテカルロ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
^ 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
離散単変量で 離散単変量で 連続単変量で 連続単変量で 連続単変量で 連続単変量で 混連続-離散単変量 多変量 (結合) 方向 退化 と特異 族 サンプリング法 （英語版 ）
