レプリカ交換モンテカルロ法（パラレル・テンパリング）による混合ガウス分布に従う乱数の生成

マルコフ連鎖モンテカルロ法（メトロポリス法）による混合ガウス分布に従う乱数の生成 - My Life as a Mock Quant
でやった内容の欠点とそれを補うためにもうちょっと高尚な手法である拡張アンサンブル法の１つ「レプリカ交換モンテカルロ法」を用いてやりましたよというお話。
シミュレーション条件・パラメーター設定に関しては適当なんで要考察。自分が扱いたい問題に応じてこれはやらないといけない。

通常のマルコフ連鎖モンテカルロ法(MCMC)での問題点

マルコフ連鎖モンテカルロ法（メトロポリス法）による混合ガウス分布に従う乱数の生成 - My Life as a Mock Quant
でやった通常のいわゆるマルコフ連鎖モンテカルロ法を用いた「混合ガウス分布」からの乱数生成方法だと、２つのガウシアンのピーク（平均値）が離れているとうまくいかない。
実際に平均値を(-3,-3),(3,3)という次候補選択の際に使用している「移動の幅」*1に比べて極端に離れているケースの場合、MCMCの特性である「局所的な探索・移動を繰り返しながらの乱数生成」は遠く離れたもう１つの混合ガウス分布のピークへ行くことができず、MCMC自体がうまくいかない。実際にコードを書いて走らせてみると↓のように片ピークからのサンプルしか得ることができない。従って、これだと「混合ガウス分布」からのサンプリングとしては不適切となる。

レプリカ交換モンテカルロ法とは

上述のように今回扱ったような「多峰性のある確率分布（ピークが複数ある確率分布）」の場合、MCMCからのサンプリングが正しいものとならないケースが存在し、ここではそれを回避するために拡張アンサンブルモンテカルロ法(extended ensemble Monte Carlo, generalized ensemble Monte Carlo)と呼ばれる手法群の１つであるレプリカ交換モンテカルロ法を用いる*2。拡張アンサンブルモンテカルロ法のアイディアとしてはパラメータの異なるいろいろな分布をまとめた分布の集まり（分布族）をシミュレーション計算の対象とし、全体の構成要素である各分布をうまく混合することで問題点を回避、実際に欲しい対象とする確率分布からのサンプリングを行う。

レプリカ交換モンテカルロ法(パラレル・テンパリング)では、異なるパラメータ*3の値 $\theta_k, k \in \left[1,\cdots,K\right]$ を持つ分布 $P(x|\theta_k)$ *4、レプリカ系を複数まとめた同時分布
$P(x_1,\cdots,x_K) = \prod_{k=1}^{K} P(x_k|\theta_k)$
からマルコフ連鎖モンテカルロ法を用いてサンプリングする。

上記の同時分布を定常分布にするような遷移として次の２ステップを考え、これらを交互に実行する。

通常のMCMCの実行： $K$ 個の各分布に対してメトロポリス法等のMCMCを実行する
レプリカの交換：適当なステップ数ごとに、ランダムに選んだ $k \in \left[1,\cdots,K\right]$ について状態 $x_k$ と $x_{k+1}$ を確率 $\min(1,r)$ で交換する（レプリカ $x_{.}$ 同士に対してメトロポリス法を適用）。

ここで、レプリカ同士の交換にはメトロポリス法のアイディアを使っているので、
$r = \frac{P(x_1,\cdots,x_{k+1},x_k,\cdots,x_K)}{P(x_1,\cdots,x_k,x_{k+1},\cdots,x_K)} = \frac{P(x_{k+1}|\theta_k)P(x_{k}|\theta_{k+1})}{P(x_k|\theta_k)P(x_{k+1}|\theta_{k+1})}$
となる。これにより分布全体が詳細釣り合いの条件を満たすので、全体として見たとき、 $P(x_1,\cdots,x_K)$ が定常分布となるようにマルコフ連鎖が設計できている。同時分布全体のシミュレーションを行いながら、自分が欲しいパラメータ値のマルコフ連鎖のみを取得すれば本来の所望の分布からのサンプリングとなる。