MCMCプログラミング
１つの乱数生成法から他の分布に従う乱数を得る方法の１つにマルコフ連鎖モンテカルロ（Markov
chain Monte Carlo; MCMC）と呼ばれているものがある。この方法では、１つの値
から次の値
をの値のみに依存しての値を決めるというマルコフ連鎖の実現値として・・・、、・・・の値を得る。このときMCMCの方法を選ぶことによって、これらの実現値の値の集合
において、実現値の時間順序を無視して実現値の分布をみると、目的とする分布の実現値とみなすことができる。MCMCの代表的な方法としてＭＨアルゴリズム（Metropolis-Hasting
algorithm）と呼ばれているものがある。これは、以下の手順のものである。
MHアルゴリズム：
（１）
を時刻
における値とするとき、
に従う分布から値を１つ選び
で表わす。すなわち
である。
（２）（１）における値に対して
における値
を次式
で与える。確率
で与えることなどアルゴリズムにおける確率的な手順は、後に示すサンプルプログラムのように一様乱数を用いて行うことができる。ここで
である。
（３）
と更新して（１）に戻る。
詳しくは、「岡本安晴『計量心理学』、2006、培風館」の「付録Ｉ マルコフ連鎖モンテカルロ」、あるいは「岡本安晴『統計学を学ぶための数学入門[上]』、2008、培風館」の「10.2 MCMC (Markov chainMonte Carlo)」などを用意したので参照されたい。ただし、岡本（2006）では、
を
と表わしている。
ベイズ的データ分析（Bayesian Data Analysis）では、target
distribution 
として事後分布をとることにより、MCMCによって事後分布の分析が行われる。MCMCはPC上で簡単に実現できるため、ベイズ的分析の実用性が飛躍的に高まった。ベイズ的データ分析については、「岡本安晴『統計学を学ぶための数学入門[下]』2009、培風館」の「9.12 確率モデルによるデータ分析」における「ベイズ的分析」、「岡本安晴『データ分析のための統計学入門』2009、おうふう」の「6.4節 ベイズ的方法」、「岡本安晴『計量心理学』2006、培風館」の「付録Ｈ ベイズ的アプローチ」で説明している。
本ホームページでは、MCMCの簡単な例を２つ用意した。１つは、
として一様分布をとり、目的とする分布を正規分布とするものである。他のものは、逆にを正規分布として、目的とする分布を一様分布とするものである。
１．一様分布から正規分布を生成するプログラムmcmc
正規分布を生成する方法は有名なものがあるが、プログラムmcmcでは敢えてMCMCによって生成してみた。MCMCに関わる個所のコードは以下のようである。
for (int i = 0; i
<= N; i++){
ux
= (rn.uni() - 0.5) * 5 + xv;
uy
= (rn.uni() - 0.5) * 5 + yv;
a
= normal2var( ux, uy ) / normal2var( xv, yv );
if (rn.uni() < a){
xv
= ux; yv = uy;
cnt++;
}
g->FillEllipse(
bshW, Math::Round(x0 + xv * pu), Math::Round(y0 + yv * pu), 3, 3 );
if ((i % 1000) == 0)
this->Text = i.ToString() + "/" + N.ToString() + "...normal distribution by MCMC/Kiss." +
"
a-level = " + (cnt/N).ToString();
}
xvとyvの値を中心とする一様乱数をuxとuyに設定してから、確率aでxvとyvの値をuxとuyに更新し、座標（xv, yv）の点を描画している。実行形式のプログラムmcmc.exeは、ファイル名mcmc.exeをクリックしてダウロード・実行することができる（.NET Framework3.5対応）。プログラム・ソース・ファイルなどをまとめたものは圧縮ファイルmcmcPrg.zipとして用意したので、ファイル名mcmcPrg.zipをクリックしてダウンロードすることができる。ダウンロードしたファイルを、マウスの右ボタンでクリックして表示されるメニュから「解凍」あるいは「展開」を選んで解凍すれば、Visual C++などから開くことができる。
プログラムmcmc.exeを実行すると図１のようなウィンドウが表示される。
図１ 実行開始時のフォーム
適当なボタンのクリックで、ボタンに表示された乱数の生成が行われ、
を座標値とする点が描画される。例えば、「WH-mcmc」ボタンをクリックするとWichmann-Hillの方法による一様乱数によって２次元正規分布に従う乱数の対が生成される。
２．正規分布から一様乱数を生成するプログラムUniFromNormal
一様乱数は直接生成されるものであるが、MCMCの説明としてとして正規分布をとって、一様分布をMCMCによって生成してみたプログラムがUniFromNormalである。MCMCに関わる個所のコードは以下のようになっている。
for (int i = 0; i
< N; i++){
z
= rn.normal()*0.8 + xt;
a
= uni_distri( z );
if (rn.uni() < a){
xt
= z;
cnt++;
}
ck_uni[ck_cat(xt)]++;
}
正規乱数をxtを中心として生成して、zに格納している。確率aでxtの値をzに更新している。更新値は関数ck_catによって分類して、カウントを格納する配列ck_uniの対応する要素の値を１つ増やしている。このプログラムのソースファイルなどはUniFromNormal.zipとしてまとめたので、ファイル名UniFromNormal.zipをクリックしてダウンロードできる。ダウンロードしたファイルは、マウスの右ボタンクリックで表示されるメニュから「解凍」あるいは「展開」を選らぶと解凍して、Visual C++まどから開くことができる。実行形式のプログラムUniFromNormal.exe（.NET Framework3.5対応）は、ファイル名UniFromNormal.exeをクリックして実行できる。プログラムを実行すると、コンソール・ウィンドウが図２のように表示される。
図２ コンソール・ウィンドウの表示
区間
を10等分したときの各部分区間にある生成乱数の比率が、配列ck_uniの要素の値として表示されている。どれも1/10とみなせる値である。
なお、Visual C++プログラミングの解説書として＜岡本安晴「大学生のための心理学VC++プログラミング入門」勁草書房＞を用意した。