簡単なギブズ サンプリング プログラム例
ギブズサンプリング（Gibbs Sampling）の基本的な考え方の説明が、J. K. Kruschke (2011).”Doing Bayesian Data Analysis: A tutorial with R and BUGS”の第８章にある。これをVisual C++/CLR 2012のコンソールアプリケーションとして作成してみた。基本的な部分は以下に示されるように簡単なコードで表すことができる。
for (int t = 1; t <= NSamples; t++){
theta1[t] = rn.beta( z1 + a1, N1 - z1 + b1 );
theta2[t] = rn.beta( z2 + a2, N2 - z2 + b2 );
thetaDiff[t] = theta1[t] - theta2[t];
}
プログラムのソースコードなどは圧縮ファイルDemoProgram.zipとしてまとめたので、クリックしてダウンロード・解凍すればVisual Studio 2012で開くことができる。
上のコードで表されるギブズサンプリングは次のデータに対してベイズ的分析を行うものである。
２枚のコイン、コイン１とコイン２、があり、それぞれのコインを投げたときの表が出る確率を
および
とおく。
の事前分布を
とおく。コイン１およびコイン２をそれぞれ
および
回投げたとき、表の出た回数がそれぞれ
および
回であったとする（このデータをDで表す）。２項分布の仮定とコイン１とコイン２の結果の独立性の仮定の下で事後確率は
で与えらえる。したがって、条件付き事後確率は
となる。これより、ギブズサンプリングが上のコードで表されることがわかる。thetaDiff[t]は、確率の差
の事後分布を求めるためのものである。
上のコードを基にして、ギブズサンプリングの様子、およびの事後分布をヒストグラムで表すプログラムを以下のように作成した。
#include "stdafx.h"
#include "RNP2to191.h"
using namespace System;
using namespace rnp2to191;
int main(array<System::String ^> ^args)
{
RNP2to191 rn;
const int NSamples = 10000;
//
// Parameter values of prior distributions
//
double a1 = 3, b1 = 3, a2 = 3, b2 = 3;
//
// Data
//
int N1 = 7, z1 = 5, N2 = 7, z2 = 2;
double theta1[NSamples + 1], theta2[NSamples + 1], thetaDiff[NSamples + 1];
theta1[0] = theta2[0] = 0.5;
thetaDiff[0] = theta1[0] - theta2[0];
for (int t = 1; t <= NSamples; t++){
//
// Gibbs sampling
//
theta1[t] = rn.beta( z1 + a1, N1 - z1 + b1 );
theta2[t] = rn.beta( z2 + a2, N2 - z2 + b2 );
thetaDiff[t] = theta1[t] - theta2[t];
//
// Displaying the sampled values
//
int pos1 = int(theta1[t] * 50), pos2 = int(theta2[t] * 50);
Console::Write(t.ToString()->PadLeft(6) + ": ");
if (pos1 < pos2){
Console::WriteLine(String("1").PadLeft(pos1) + String("2").PadLeft(pos2 - pos1));
}else if (pos1 > pos2){
Console::WriteLine(String("2").PadLeft(pos2) + String("1").PadLeft(pos1 - pos2));
}else{
Console::WriteLine(String("X").PadLeft(pos1));
}
}
Console::WriteLine();
double sum = 0.0;
int n_greater0;
int freq[27];
for (int i = 0; i < 27; i++) freq[i] = 0;
for (int t = 1; t <= NSamples; t++){
//
// Sum of thetaDiff = theta1 - theta2
//
sum += thetaDiff[t];
//
// Number of thetaDiff > 0
//
if (thetaDiff[t] > 0.0) n_greater0++;
//
// Counting values for a histogram
//
int pos = int((thetaDiff[t] + 0.4 + 0.025) / 0.05);
if (pos < 0.0){
freq[0]++;
}else if(pos > 24){
freq[26]++;
}else{
freq[pos+1]++;
}
}
double mean = sum / NSamples;
Console::WriteLine("Mean = " + mean.ToString());
Console::WriteLine("Prob( diff. > 0.0) = " + (double(n_greater0) / NSamples).ToString("F5"));
Console::WriteLine();
Console::WriteLine("Histogram of theta1 - theta2");
int nmax = 0;
for (int i = 0; i < 26; i++)
if (nmax < freq[i]) nmax = freq[i];
for (int ix = 1; ix <= 25; ix++){
Console::Write((-0.4 + (ix - 1)*0.05).ToString("F2")->PadLeft(6) + " I");
for (int i = 0; i < int(freq[ix] * 50.0 / nmax); i++)
Console::Write("X");
Console::WriteLine();
}
Console::WriteLine();
Console::WriteLine("Press Enter key to exit.");
Console::ReadLine();
return 0;
}
上のプログラムを実行すると図１のようになる。
図１
図１におけるヒストグラムは、Kruschke (2011)のFigure 8.8 に対応するものである。ヒストグラムの前に
「9999: 2 1」
などとあるのは、9999回目のサンプル値の位置を横方向に図示するものである。２はtheta2[9999]の値、１はtheta1[9999]の値の位置を表している。２とか１でなく、Xが記されているのときは、１と２が重なっていることを表す。
図１の結果は、実行形式のプログラムDemoProgram.exeをダウンロード・実行することによっても得られる。