MCMCの簡単な例
Linux用のサンプルプログラムのウェブサイトも用意しました。
MCMC（Markov chain Monte Carlo）の簡単な例が＜J. K. Kruschke (2011). “Doing Bayesian Data Analysis: A tutorial with R and BUGS”, Chap. 7＞で説明されている。この説明に従ったプログラム例をVisual C++/CLR 2012 コンソールアプリケーションとして作成してみた。ターゲットとする確率分布は
であり、提案分布は
である。このとき、現在の値
から提案分布で与えられるサンプル値
に値を更新する確率
は
で与えられる（岡本安晴（2006）「計量心理学」培風館、付録Ｉ「マルコフ連鎖モンテカルロ」参照）。
上のことをプログラムで表すと以下のようになる。
#include "stdafx.h"
#include "RNP2to191.h" // 乱数生成用クラス型 RNP2to191 が宣言されている
using namespace System;
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(array<System::String ^> ^args)
{
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 = Math::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]++;
Console::Write((istep+1).ToString() + ":");
for (int i = 0; i < mcmc_pos[istep+1]; i++) Console::Write(" ");
Console::WriteLine(mcmc_pos[istep+1].ToString());
}
Console::WriteLine();
for (int i = 0; i < 7; i++){
Console::Write("Position-" + (i+1).ToString() + ": p*28 = " + (28.0 * count[i] / double(nsteps)).ToString("F3") + " I");
for (int j = 0; j < count[i] / 300; j++)
Console::Write("x");
Console::WriteLine();
}
Console::WriteLine();
Console::Write("Enterキーを押して終了。");
Console::ReadLine();
return 0;
}
乱数のクラス型 RNP2to191 は、ヘッダーファイル RNP2to191.h で宣言されているが、この乱数についてはウェブサイト＜疑似乱数の生成＞で説明している。
上のプログラムを実行すると図１のような結果となる。
図１
上の出力において
「28000: 6」
とあるのは、28000ステップにおける値が６であることを表している。MCMCによる各ステップにおける値の遷移が視覚的に表示されている。
「Position-1: p*28 = 1.081 Ixxx」
とあるのは、ステップにおける値が１であった比率（確率）を28倍すると1.081になることを表している。28倍しているのは、式（１）との対応を分かり易くするためである。Position-iのiの値が１から７に変わるに伴って、p*28 の値がほぼ１から７まで変化していることがわかる。すなわち、上のＭＣＭＣのプログラムによって式（２）の分布に合うサンプルが得られたことが示されている。「Ixxx」は横方向の棒グラフである。
プログラムのソースコードファイルなどは圧縮ファイルConsoleApplication.zipとして用意した。クリックしてダウロード解凍すれば、Visual Studio 2012で開くことができる。Visual Studio 2012の無料版が用意されているが、ウェブサイト＜Windowsフォームアプリケーションの作成＞で説明している。Visual C++/CLRでのコンソールアプリケーションの作成については＜VC++入門:コンソールアプリケーションの作成＞で説明している。
実行形式のファイルも用意しているので、プログラムファイル名ConsoleApplication.exeをクリックしてダウンロード実行することができる。ただし、.NET Framework 4.5がインストールされていないとエラーになるので注意。