MCMCの簡単なC++プログラム例
MCMC（Markov chain Monte Carlo）の簡単な例が＜J. K. Kruschke (2011). “Doing Bayesian Data Analysis: A tutorial with R and BUGS”, Chap. 7＞で説明されている。この説明に従ったプログラム例をC++/Linux(Ubuntu)で作成してみた。ターゲットとする確率分布は
であり、提案分布は
である。このとき、現在の値
から提案分布で与えられるサンプル値
に値を更新する確率
は
で与えられる（岡本安晴（2006）「計量心理学」培風館、付録Ｉ「マルコフ連鎖モンテカルロ」参照）。
上のことをプログラムで表すと以下のようになる。
=======================================================================================================
#include <cmath>
#include <iostream>
#include "RNP2to191.h" // 64ビットLinux(Ubuntu12)用の乱数生成用クラス型 RNP2to191
が宣言されている
// 32ビットLinuxのときは、RNP2to191_32.hを用いる
using namespace std;
using namespace rnp2to191;
RNP2to191 rn; // 乱数生成用オブジェクト
int proposal( int icurrent ){
if
(rn.uni() < 0.5) return icurrent
- 1;
else
return icurrent + 1;
}
double TargetProb( int
ipos ){
if
((ipos < 1) || (ipos
> 7)){
return
0.0;
}else{
return
ipos / 28.0;
}
}
int main(){
const int nsteps
= 28000;
int mcmc_pos[nsteps
+ 1];
mcmc_pos[0] = 4; // 初期値
int count[] = { 0, 0, 0, 0, 0, 0, 0 }; // 各状態をとった回数の初期化
for
(int istep = 0; istep < nsteps; istep++){
int proposed = proposal( mcmc_pos[istep] );
double
a = min( 1.0, TargetProb( proposed ) / TargetProb( mcmc_pos[istep] ) );
if
(rn.uni() < a){
mcmc_pos[istep+1] = proposed;
}else{
mcmc_pos[istep+1] = mcmc_pos[istep];
}
count[mcmc_pos[istep+1]-1]++;
cout << (istep+1) << ":";
for
(int i = 0; i < mcmc_pos[istep+1]; i++)
cout << " ";
cout << mcmc_pos[istep+1]
<< endl;
}
cout << endl;
for
(int i = 0; i < 7; i++){
cout << "Position-" << (i+1) <<
": p*28 = "
<<
(28.0 * count[i] / double(nsteps))
<< " I";
for
(int j = 0; j < count[i]
/ 300; j++)
cout << "x";
cout << endl;
}
cout << endl;
return
0;
}
=======================================================================================================
乱数のクラス型 RNP2to191 は、ヘッダーファイル RNP2to191.h （64ビットOS用、32ビットOS用はRNP2to191_32.h）で宣言されているが、この乱数についてはウェブサイト＜疑似乱数の生成＞で説明している。
上のプログラムを実行すると図１のような結果となる。
図１
上の出力において
「28000: 6」
とあるのは、28000ステップにおける値が６であることを表している。MCMCによる各ステップにおける値の遷移が視覚的に表示されている。
「Position-1: p*28 = 1.081 Ixxx」
とあるのは、ステップにおける値が１であった比率（確率）を28倍すると1.081になることを表している。28倍しているのは、式（１）との対応を分かり易くするためである。Position-iのiの値が１から７に変わるに伴って、p*28 の値がほぼ１から７まで変化していることがわかる。すなわち、上のＭＣＭＣのプログラムによって式（２）の分布に合うサンプルが得られたことが示されている。「Ixxx」は横方向の棒グラフである。
プログラムのソースコードファイルなどは圧縮ファイルMCMCSamplePrg64.zip（64ビットOS用）あるいはMCMCSamplePrg32.zip（32ビットOS用）として用意した。クリックしてダウロード解凍すれば、自由に利用できる。解凍したファイルには、makefileも含まれているのでmakeコマンドで実行形式ファイルを作ることもできる。