上記の広告は１ヶ月以上更新のないブログに表示されています。新しい記事を書く事で広告が消せます。

以前に乱数発生ライブラリrandlibについて ここに書いたのですが、そのライブラリの乱数のシードを変更する
方法が分かったので、メモしておきます。

シードを変更しないと、毎回同じ乱数が生成されてしまうので

シードの変更にはsetall(long iseed1,long iseed2)という関数を使います。
これですべての乱数発生パラメータに引数の値を設定するのだとか。
(iseed1とiseed2の違いについては、マニュアルには載っていませんでした)

他の関数でvoid setsd(long iseed1,long iseed2)というのもあるのですが、
この関数をコールするにはその前にsetall()を実行していないといけないらしいので、
setall()のみをつかえばいいかな・・・と。

実際には↓こんな感じで使っています。

#include "randlib.h"
#include <stdio.h>
#include <time.h>

int main( void )
{

    setall( (long)time(NULL), (long)time(NULL) );


    double      value1 = genunf( 0.0, 1.0 );

    printf( "TIME = %ld VALUE = %f\n", (long)time(NULL), value1 );

    return 0;
}

これだと、実行するたびに乱数が異なります
（1秒間に何回もプログラムを起動した場合はおなじですけど）

なんだかあまり逆カイ2乗分布について載っているサイトが少ないので、
「Bayesian Data Analysis(Second Edition)」(Andrew Gelman et al.)から逆カイ2乗分布についてメモ。

昨日はrandom walk chainでしたが、今日はindependence chainです。
どこが違うのかというと、昨日のは候補の発生は現在の点+αで行っていましたが、
今日のは候補の発生は現在の点とは独立して行います。

使用する例は昨日のと同じものを使います。ただし候補を発生させる提案分布は以下のものを使います。

　　　　　　　　　　

【メトロポリスヘイスティングス法(independence chainアルゴリズム)の手順】



この手順を基にして作ったプログラムと結果が以下になります。
ちなみにaccept率は77.5%でした。割と高いっすね。

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "randlib.h"

#define PI			3.14159265358979324

// 平均値の初期パラメータ
double		priMean = 0.0;
// 分散の初期パラメータ
double		priVariance = 1000;

// t分布に従う乱数を生成する
double tRandom( double df )
{
	return gennor( 0.0, 1.0 ) / sqrt( genchi( df ) / df );
}


// 標準正規分布の確率密度を求める
double stdNormalFunction( double x )
{
	double	a = 1 / sqrt( 2 * PI );
	double	b = exp( - 0.5 * ( x * x ) );

	return a * b;
}


// ガンマ関数の値を求める
double gammln( double xx ) 
{
	double		x;
	double		y;
	double		tmp;
	double		ser;

	static double	cof[6] = {
		76.18009172947146,
		-86.50532032941677,
		24.01409824083091,
		-1.231739572450155,
		0.1208650973866179e-2,
		-0.5395239384953e-5 };
	
	int			j;

	y = xx;
	x = xx;
	tmp = x + 5.5;
	tmp -= ( x + 0.5 ) * log( tmp );
	ser = 1.000000000190015;
	for( j = 0; j <= 5; j ++ ) {
		ser += cof[j] / ++ y;
	}

	return - tmp + log( 2.5066282746310005 * ser / x );
}


// ベータ関数の値を求める
double betaFunction( double a, double b )
{
	return exp( gammln( a ) + gammln( b ) - gammln( a + b ) );
}


// t分布の確率密度を求める
double tFunction( double x, double df )
{
	double	a = pow( ( 1.0 + ( x * x ) / df ), - 0.5 * ( df + 1 ) );
	double	b = sqrt( df ) * betaFunction( 0.5, 0.5 * df );

	return a / b;
}




int main( void )
{
	// t分布の自由度
	double			df = 5.0;
	// t分布の平均
	double			mean = 5.0;
	// t分布の分散
	double			variance = 4.0;
	// 観測値数
	int				dataNum = 1000;
	// 観測値格納域
	double			y[dataNum];

	int				i;
	int				j;
	double			ran;
	// サンプリングした平均値
	double			mu;
	// 平均値の候補
	double			candiMean;
	// 観測値の平均
	double			xbar = 0.0;
	// サンプリング回数
	int				replaceNum = 10000;
	// burn-in回数
	int				burnIn = 1000;
	// 現状と候補の確率比
	double			r;
	// 平均のサンプリング結果格納域
	double			mean_list[replaceNum];
	// accept回数
	int				acceptNum = 0;
	double			nowStdErr = 0.0;
	double			candiStdErr = 0.0;
	
	// 観測値をt分布(v=5,mu=5,sigma=4)から1000個の乱数発生で生成する
	for( i = 0; i < dataNum; i ++ ) {
		y[i] = mean + sqrt( variance ) * tRandom( df );
		xbar += y[i];
	}
	xbar /= (double)dataNum;

	// 平均値の初期値は0とする
	mu = 0;
	
	// サンプリング開始
	for( i = - burnIn; i < replaceNum; i ++ ) {
		// 候補を求める
		candiMean = xbar + sqrt( df * variance / (double)dataNum ) * 
					gennor( 0.0, 1.0 );

		// 候補をacceptする確率を求める
		nowStdErr = 0.0;
		candiStdErr = 0.0;
		for( j = 0; j < dataNum; j ++ ) {
			candiStdErr += log( 1 / sqrt( variance ) * 
				tFunction( ( y[j] - candiMean ) / sqrt( variance ), df ) );
			nowStdErr += log( 1 / sqrt( variance ) *
				tFunction( ( y[j] - mu ) / sqrt( variance ), df ) );
		}

		r = ( 1 / sqrt( priVariance ) ) * 
			stdNormalFunction( ( candiMean - priMean ) / sqrt( priVariance ) ) /
			( ( 1 / sqrt( priVariance ) ) *
			stdNormalFunction( ( mu - priMean ) / sqrt( priVariance ) ) ) *
			exp( candiStdErr - nowStdErr ) *
			1 / sqrt( df * variance / dataNum ) * 
			stdNormalFunction( ( mu - xbar ) / sqrt( df * variance / dataNum ) ) /
			( 1 / sqrt( df * variance / dataNum ) *
			stdNormalFunction( ( candiMean - xbar ) / sqrt( df * variance / dataNum ) ) );

		// accept or reject判定
		ran = gennor( 0.0, 1.0 );
		if( r >= ran ) {
			// 候補をacceptする
			mu = candiMean;
			acceptNum ++;
		}

		if( i >= 0 ) {
			mean_list[i] = mu;
		}
	}

	// 結果出力
	printf( "sampling number\tmean\n" );
	mu = 0.0;
	for( i = 0; i < replaceNum; i ++ ) {
		mu += mean_list[i];
		printf( "%d\t\t%f\n", i + 1, mean_list[i] );
	}
	mu /= (double)replaceNum;
	printf( "posterior mean = %f, accept rate = %f\n",
			mu,
			(double)acceptNum / ( (double)burnIn + (double)replaceNum ) );

	return 0;
}

今日はメトロポリス・ヘイスティングス法について。
その中でもrandom walk chain(酔歩連鎖)についてです。
だんだんめんどくさくなってきたので、今日はいきなり例を出していきます。



【メトロポリスヘイスティングス法(random walk chainアルゴリズム)の手順】



この手法の長所はアルゴリズムの構築が簡単な点ですが、短所は提案分布のチューニングが
必要な点です。現在の点と候補の距離（分散などのちらばり）が小さくなるように提案分布を設定すると
候補のaccept率は上がりますが状態空間の一部ばかりをサンプリングしてしまいます。
逆にばらつきが大きいとaccept率が低下してサンプリング効率は悪くなってしまいます。

この手法のプログラムは以下のようになります↓

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "randlib.h"

#define PI		3.14159265358979324


// 平均値の初期パラメータ
double		priMean = 0.0;
// 分散の初期パラメータ
double		priVariance = 1000;

// t分布に従う乱数を生成する
double tRandom( double df )
{

	return gennor( 0.0, 1.0 ) / sqrt( genchi( df ) / df );
}


// 標準正規分布の確率密度を求める
double stdNormalFunction( double x )
{
	double	a = 1 / sqrt( 2 * PI );
	double	b = exp( - 0.5 * ( x * x ) );

	return a * b;
}


// ガンマ関数の値を求める
double gammln( double xx )
{
	double		x;
	double		y;
	double		tmp;
	double		ser;

	static double cof[6] = {
		 76.18009172947146,
		-86.50532032941677,
		 24.01409824083091,
		 -1.231739572450155,
		  0.1208650973866179e-2,
		 -0.5395239384953e-5 };

	int			j;

	y = xx;
	x = xx;
	tmp = x + 5.5;
	tmp -= ( x + 0.5 ) * log( tmp );
	ser = 1.000000000190015;
	for( j = 0; j <= 5; j ++ ) {
		ser += cof[j] / ++ y;
	}

	return - tmp + log( 2.5066282746310005 * ser / x );
}

// ベータ関数の値を求める
double betaFunction( double a, double b )
{
	return exp( gammln( a ) + gammln( b ) - gammln( a + b ) );
}


// t分布の確率密度を求める
double tFunction( double x, double df )
{
	double	a = pow( ( 1.0 + ( x * x ) / df ), - 0.5 * ( df + 1 ) );
	double	b = sqrt( df ) * betaFunction( 0.5, 0.5 * df );

	return a / b;
}


int main( void )
{
	// t分布の自由度
	double		df = 5.0;
	// t分布の平均
	double		mean = 5.0;
	// t分布の分散
	double		variance = 4.0;
	// 観測値数
	int			dataNum = 1000;
	// 観測値格納域
	double		y[dataNum];

	int			i;
	int			j;
	double		ran;
	// サンプリングした平均値
	double		mu;
	// 平均値の候補
	double		candiMean;
	// サンプリング回数
	int			replaceNum = 10000;
	// burn-in回数
	int			burnIn = 1000;
	// 提案分布のパラメータ
	double		tausq = 0.2;
	// 現状と候補の確率比
	double		r;
	// 平均のサンプリング結果格納域
	double		mean_list[replaceNum];
	// accept回数
	int			acceptNum = 0;
	double		nowStdErr = 0.0;
	double		candiStdErr = 0.0;

	// 観測値をt(v=5,mu=5,sigma=4)から1000個の乱数発生で生成する
	for( i = 0; i < dataNum; i ++ ) { 
		y[i] = mean + sqrt( variance ) * tRandom( df );
	}

	// 平均値の初期値をN(0,1000)から発生させる
	mu = gennor( priMean, sqrt( priVariance ) );
	for( i = 0; i < dataNum; i ++ ) {
		nowStdErr += log( 1 / sqrt( variance ) * 
				tFunction( ( y[i] - mu ) / sqrt( variance ), df ) );	
	}

	// サンプリング開始
	for( i = - burnIn; i < replaceNum; i ++ ) {
		// 候補を求める
		candiMean = mu + gennor( 0.0, sqrt( tausq ) );

		// 現状と候補の事後確率密度の比を求める	
		candiStdErr = 0.0;
		for( j = 0; j < dataNum; j ++ ) {
			candiStdErr += log( 1 / sqrt( variance ) *
				tFunction( ( y[j] - candiMean ) / sqrt( variance ), df ) );
		}
		r = ( 1 / sqrt( priVariance ) ) * 
			stdNormalFunction( ( candiMean - priMean ) / sqrt( priVariance ) ) /
			( ( 1 / sqrt( priVariance ) ) * 
			stdNormalFunction( ( mu - priMean ) / sqrt( priVariance ) ) ) *
			exp( candiStdErr - nowStdErr ); 

		// accept or reject判定
		while( 1 ) {
			ran = gennor( 0.0, 1.0 );
			if( ran != 0.0 && ran != 1.0 ) {
				break;
			}
		}
		if( r >= ran || r >= 1.0 ) {
			// 候補をacceptする
			mu = candiMean;
			nowStdErr = candiStdErr;
			acceptNum ++;
		}

		if( i >= 0 ) {
			mean_list[i] = mu;
		}
	}

	// 結果出力
	printf( "sampling number\tmean\n" );
	mu = 0.0;
	for( i = 0; i < replaceNum; i ++ ) {
		mu += mean_list[i];
		printf( "%d\t\t%f\n", i + 1, mean_list[i] );
	}
	mu /= (double)replaceNum;
	printf( "posterior mean = %f, accept rate = %f\n", 
			mu, 
			(double)acceptNum / ( (double)burnIn + (double)replaceNum ) );

	return 0;
}

このプログラムで提案密度のパラメータを0.2としたときのサンプリングの様子をプロットしてみました↓
このときの事後平均の結果は5.06で、accept率は71.2%でした。



今度はパラメータを0.001に変更してサンプリングしてみました↓
このときのaccept率は83.4%と高くなりましたが、真の値にたどり着くまでのサンプリングの動きが緩やか
になっています。



今度はパラメータを4.0にしてサンプリングしてみました↓
このときのaccept率は70.4%で、サンプリング幅が大きくなっています。

いよいよGibbs Samplingについてのメモです。
教科書などでは「ギブスサンプラー」と書かれている方が多いのですが、
私はGibbs Samplingで習ったのでここではこれで通します。



【Gibbs Samplingの手順】

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "randlib.h"








int main( void )
{

	// 母集団の平均値
	double			trueMean = 5.0;
	// 母集団の分散
	double			trueVar = 1.0;
	// 観測値数
	int				dataNum = 1000;
	// 観測値格納域
	double			y[dataNum];
	// 観測値の平均
	double			xbar = 0.0;
	// 観測値の分散
	double			xvar = 0.0;

	// 事前分布/条件付き事後分布のパラメータ
	double			mu = 0.0;
	double			m = 1.0;
	double			alpha = 0.02;
	double			beta = 0.02;

	int				i;
	// サンプリング回数
	int				replaceNum = 10000;
	// burn-in回数
	int				burnIn = 1000;
	// 平均のサンプリング結果格納域
	double			mean_list[replaceNum];
	// 分散のサンプリング結果格納域
	double			variance_list[replaceNum];

	// 観測値をN(5,1)から100個の乱数発生で生成する
	for( i = 0; i < dataNum; i++ ) {
		y[i] = gennor( trueMean, sqrt( trueVar ) );
	}

	// 観測値の平均を求める
	for( i = 0; i < dataNum; i ++ ) {
		xbar += y[i];
	}
	xbar /= dataNum;

	// 観測値の2乗差を求める
	for( i = 0; i < dataNum; i ++ ) {
		xvar += ( y[i] - xbar ) * ( y[i] - xbar );
	}

	// 分散の初期値を事前分布IG(alpha/2, beta/2)から生成する
	double			estVar = 1 / gengam( alpha / 2, beta / 2 );

	// 平均の初期値を事前分布N(0,σ 2/m)から生成する	
	double			estMean = gennor( mu, estVar / m );

	// 各パラメータを更新
	beta += xvar + ( dataNum * m * pow( ( xbar - mu ), 2.0 ) ) /
			( dataNum + m );
	mu = ( m * mu + dataNum * xbar ) / ( m + dataNum );
	alpha += dataNum;
	m += dataNum;

	// ギブスサンプリング開始
	for( i = -burnIn; i < replaceNum; i ++ ) {

		// 平均のサンプリングを条件付き事後分布から求める
		estMean = gennor( mu, estVar / m );

		// 分散のサンプリングを条件付き事後分布から求める
		estVar = 1 / 
				 gengam( ( alpha + 1.0 ) / 2, 
						( m * pow( ( estMean - mu ), 2.0 ) + beta ) / 2 ); 

		if( i >= 0 ) {
			mean_list[i] = estMean;
			variance_list[i] = estVar;
		}
	}

	// 結果出力
	printf( "sampling number\tmean\tvariance\n" );
	estMean = 0.0;
	estVar = 0.0;
	for( i = 0; i < replaceNum; i ++ ) {
		estMean += mean_list[i];
		estVar += variance_list[i];
		printf( "%d\t\t%f\t\t%f\n", i + 1, mean_list[i], variance_list[i] );
	}
	estMean /= (double)replaceNum;
	estVar /= (double)replaceNum;
	printf( "posterior mean : mean = %f, variance = %f\n", 
			estMean,
			estVar );

	return 0;
}

事後平均の結果は平均が4.965364で分散が1.037286で、そこそこ正確に推定することができました。
したのプロットは上が平均のサンプリングの様子で下が分散のサンプリングの様子です。

引き続きマルコフ連鎖モンテカルロ法です。
学部生の頃はまさか自分がMCMC法をこんな勉強することになるとは夢にも思わなかったので
（なんかよくわからんけど文献に時々出てくるスゴイmethodみたいな認識でした）、
なんだかちょっと感慨深けです

今日はA-R（acceptance-rejection）法をやります。自分がある密度関数 から乱数を生成したいけど、分布関数が非常に複雑で乱数発生が困難な時、正規分布などの簡単な確率分布の密度関数 を使って目的の確率分布の乱数を生成する方法です。

【A-R法の手順】
１．確率密度関数 から乱数 を生成する。

２．区間(0,1)から乱数 を生成する。

３．乱数 のaccept/reject判定基準として以下のrを計算する。

　　　　　　　　

ここでcは であるすべてのxについて が成り立つような正の定数である。

４． ならば　はrejectされ、それ以外ならばacceptされる。

この手法の長所としては　の基準化定数を求めなくて済むこと（メトロポリス法も求めなくてよい）だそうですが、短所としては に近似するを見つけなくてはいけないことと、cの値を大きくしすぎるとサンプリング効率が悪くなってしまうことです。そう考えると、ちょっとめんどーな方法と言えなくもないなー

以下のソースは自由度1のカイ2乗分布を用いてガンマ分布(2,1)からサンプリングを発生させています。cの値は3とし、各確率密度関数は以下の通りになります。

　　　　　　　　

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "randlib.h"

#define PI	3.14159265358979324

// 自由度1のχ 2乗分布の確率密度を求める
double chisqFunction( double x )
{
	return ( pow( x, -0.5 ) * exp( -0.5 * x ) ) / sqrt( 2 * PI );
}



// ガンマ分布(2,1)の確率密度を求める
double gammaFunction( double x )
{
	return x * exp( - x );
}


int main( void )
{
	// 反復回数
	int			replaceNum = 10000;
	// Cの値
	const double	C = 3.0;

	// xの値
	double		x_list[replaceNum];
	// サンプリングで得られたガンマ分布の乱数
	double		gamma_list[replaceNum];
	// c*χ２乗分布の乱数
	double		chisq_list[replaceNum];

	int			i;
	double		x;
	double		gamma;
	double		chisq;
	double		r;
	double		ran;
	int			cnt = 0;

	for( i = 0; i < replaceNum; i ++ ) {
		// 自由度1のχ 2乗分布から乱数を発生させる
		x = genchi( 1 );

		// (0,1)区間から一様乱数を生成する
		while( 1 ) {
			ran = genunf( 0.0, 1.0 );
			if( ran != 0.0 && ran != 1.0 ) {
				break;
			}
		}

		// 確率変数xのガンマ分布の確率密度を求める
		gamma = gammaFunction( x );

		// 確率変数xのχ 2乗分布の確率密度を求める
		chisq = chisqFunction( x );

		r = gamma / ( C * chisq );

		if( ran <= r ) {
			// 乱数xをacceptする
			x_list[cnt] = x;
			gamma_list[cnt] = gamma;
			chisq_list[cnt] = C * chisq;
			cnt ++;
		}
	}

	printf( "x\t\tgamma\t\tchisq\n" );
	for( i = 0; i < cnt; i ++ ) {
		printf( "%f\t%f\t%f\n", x_list[i], gamma_list[i], chisq_list[i] );
	}

	printf( "%d  RATE = %f\n", cnt, ( (double)cnt / (double)replaceNum ) );

	return 0;
}

以下のグラフは上記A-R法で得られたガンマ分布(2,1)の乱数のプロットです。
ちゃんとガンマ分布(2,1)に従っていることが分かります（　´_ゝ`）



次のグラフはA-R法で得られたガンマ分布(2,1)の乱数とカイ2乗分布のプロットです。
これを見るとxが0.0～1.0あたりのとき2つの分布の近似が良くないためにサンプリング効率が悪くなっています。ちなみに上記プログラムでは10000回反復のうち、acceptされたのは 約3300回でした。

今日は変数の範囲が不等式で制約されている例について、メトロポリス法を適用してみました。
変数が従う分布は前回と同じで、 には以下のような制約が存在し、



その範囲外での確率は0とします。このとき、メトロポリス法は以下のような手順で実行します。
（ちなみに１．～３．までは制約なしのメトロポリスと同じです）

【2変量正規分布の制約ありのメトロポリス法の手順】
１． の任意の初期値を設定する。

２． を原点対照の確率密度 から発生させる。

３．「候補」を以下のように設定する。

　　　　　　　

４． ならば候補はrejectされ、手順２．に進む。acceptされる場合は５．に進む。

５．候補の状態と現在の状態の確率の比を計算する。

　　　　　　　　

６． の一様乱数を生成する。 ならば候補の値を現在の値とし、それ以外はrejectする。

７．２．～６．を繰り返す。

　手順４．はとりあえず候補を選んで、範囲外ならrejectするというアルゴリズムが正しい。範囲外にならないように候補を選ぶといった方法を使ってはいけない。

上記メトロポリス法の例からサンプルの様子を出力してみました↓

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "randlib.h"

// 確率p(x1,x2)を計算する
// 正規化定数は計算しない
double cal_pro( double x1, double x2, double b )
{
	double	val = -0.5 * ( x1 * x1 - 2 * b * x1 * x2 + x2 * x2 );
	double	ans = exp( val );

	return ans;
}



int main( void )
{
	// x1, x2の初期値
	double		x1 = 2.0;
	double		x2 = 7.0;

	// 新しいx1, x2
	double		nx1;
	double		nx2;

	// 反復数
	int			replaceNum = 3200;

	// 固定パラメータ
	double		b = 0.8;
	// ステップ幅
	double		delta = 1.0;
	// 前回の確率
	double		prior_pro = cal_pro( x1, x2, b );
	// 次回の確率
	double		post_pro;
	// x1の格納域
	double		x1_list[replaceNum];
	// x2の格納域
	double		x2_list[replaceNum];

	int			i;
	int			j;
	double		r;
	double		ran;
	double		x1mean;
	double		x2mean;

	printf( "x1mean\t\tx2mean\n" );
	for( i = 0; i < replaceNum; i ++ ) {
		// x1の候補を求める
		// nx1 = x1 + (平均0分散delta^2の正規乱数)
		nx1 = x1 + gennor( 0.0, delta * delta );
		// x2の候補を求める
		// nx2 = x2 + (平均0分散delta^2の正規乱数)
		nx2 = x2 + gennor( 0.0, delta * delta );

		// x1の候補とx2の候補がnx2 >= 0.5*nx1^2+3を満たすかどうか
		// 判定する
		double	rec = 0.5 * ( nx1 * nx1 ) + 3.0;
		if( nx2 >= rec ) {
			// 制約を満たしているので確率の比から
			// 候補のaccept/rejectを判定する

			// 候補の確率を求める
			post_pro = cal_pro( nx1, nx2, b );
			// 現在と候補の確率の比を計算する
			r = post_pro / prior_pro;

			// 候補をacceptするかrejectするか判定する
			ran = (double)rand() / ( (double)RAND_MAX + 1.0 );
			if( ran < r ) {
				// 候補をacceptする
				x1 = nx1;
				x2 = nx2;
			}
		}
		// rejectされた場合はそのまま
		x1_list[i] = x1;
		x2_list[i] = x2;
		x1mean = 0.0;
		x2mean = 0.0;
		for( j = 0; j < i; j ++ ) {
			x1mean += x1_list[j];
			x2mean += x2_list[j];
		}
		x1mean /= (double)( i + 1 );
		x2mean /= (double)( i + 1 );
		printf( "%f\t\t%f\n", x1mean, x2mean );
		prior_pro =post_pro;
	}
	
	printf( "x1\t\tx2\n" );
	for( i = 0; i < replaceNum; i ++ ) {
		printf( "%f\t\t%f\n", x1_list[i], x2_list[i] );
	}

	return 0;
}

上記グラフで、点線部分はです。これを見ると範囲内にちゃんとサンプリングされていることが分かります。
以下のグラフは平均値の動く様子です。平均値は1点に収束されていくのが分かります。

前回は離散変数に対してメトロポリス法を使ったので、今回は連続変数に対してメトロポリス法を やってみます。
2変量正規分布



から各確率変数をメトロポリス法によってサンプリングします。ここで

は正規化定数ですが、候補と現在の確率の比を計算するときに相殺されるので、ここでは計算しません。

【2変量正規分布のメトロポリス法の手順】
１． の任意の初期値を設定する。

２． を原点対称の確率密度から発生させる。

３．「候補」を以下のように設定する。

　　　　　　

４．「候補」の状態と現在の状態の確率の比を計算する。

　　　　　　　

５． の一様乱数を生成する。ならば、「候補」の値を現在の値とし、それ以外は何もしない。

６．２．～５．を繰り返す。

ここで確率密度はいろいろな方法がありますが、私は を適当な「ステップ幅」として、
平均０、分散の正規分布からを独立に選ぶという方法をとりました。

上記のメトロポリス法の例を使って、サンプリングを行うプログラムを作成してみました。
ちなみに、 の初期値は(3.0,9.0)とし、は0.8としました。また、は0.8としました。

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include "randlib.h"

// 確率p(x1,x2)を計算する
// 正規化定数は計算しない
double cal_pro( double x1, double x2, double b )
{
	double	val = -0.5 * ( x1 * x1 - 2 * b * x1 * x2 + x2 * x2 );
	double	ans = exp( val );

	return ans;
}



int main( void )
{
	// x1, x2の初期値
	double		x1 = 3.0;
	double		x2 = 9.0;

	// 新しいx1, x2
	double		nx1;
	double		nx2;

	// 反復数
	int			replaceNum = 800;

	// 固定パラメータ
	double		b = 0.8;
	// ステップ幅
	double		delta = 0.8;
	// 前回の確率
	double		prior_pro = cal_pro( x1, x2, b );
	// 次回の確率
	double		post_pro;

	int			i;
	int			j;
	double		r;
	double		ran;

	printf( "x1\t\tx2\n" );
	for( i = 1; i <= replaceNum; i ++ ) {
		// x1の候補を求める
		// nx1 = x1 + (平均0分散delta^2の正規乱数)
		nx1 = x1 + gennor( 0.0, delta * delta );
		// x2の候補を求める
		// nx2 = x2 + (平均0分散delta^2の正規乱数)
		nx2 = x2 + gennor( 0.0, delta * delta );

		// 候補の確率を求める
		post_pro = cal_pro( nx1, nx2, b );
		// 現在と候補の確率の比を計算する
		r = post_pro / prior_pro;

		// 候補をacceptするかrejectするか判定する
		ran = (double)rand() / ( (double)RAND_MAX + 1.0 );
		if( ran < r ) {
			// 候補をacceptする
			x1 = nx1;
			x2 = nx2;
		}
		// rejectされた場合はそのまま
		printf( "%f\t\t%f\n", x1, x2 );
		prior_pro =post_pro;
	}

	return 0;
}

サンプリング結果をプロットすると、以下のようになりました。

昨日からメトロポリス・ヘイスティングス法を勉強しています
昨日はマルコフ連鎖モンテカルロ法の1つであるメトロポリス法を勉強しました。
今回の教科書は計算統計Ⅱです。以前にGibbsSamplingもこれで勉強したのですが、なかなか分かりやすいと思います。

そもそもマルコフ連鎖モンテカルロ法とは複雑な分布に従う確率変数についてパラメータ(期待値や
分散など）が不明なとき、データが与えられた時のパラメータのもとでのデータの確率とパラメータの
事前分布を使って、パラメータの条件付き事後確率からサンプリングしたデータからパラメータを
求める方法です。

要はちょっとずつ事後分布を変更しながら何回もサンプリングを行い、そのサンプリングから
パラメータを求めるという方法です。

メトロポリス法については3つの離散確率変数がそれぞれ1もしくは-1をとるときを例として
説明してありました。

【メトロポリス法の手順】
１．3つの離散確率変数 における分布が



のとき、変数のどれか1つを選ぶ。

２．選んだ変数の値をで置き換えて、他の2つの変数をそのまました状態を次のサンプリングの「候補」
とする。
ex.　状態でを選んだとすると、候補はとなる。

３．「候補」の確率と現在の状態の確率の比rを求める。



４． となる一様乱数Ｒを生成する。

５． なら「候補」を次の状態に採用する(accept)。 なら現在の状態をそのまま次の状態とする
(reject)。

　この方法では「候補」の確率が現在より大きい場合は候補は必ず採用される。また、「候補」の確率が
現在より小さい場合は確率rで候補は採用される。

　このメトロポリス法を十分に繰り返した後に得られるサンプルを使っての平均や分散を計算する手法をマルコフ連鎖モンテカルロ法という。

　メトロポリス法によって得られたサンプル間には相関が存在するため、分散や精度の高いパラメータの
計算を行う場合は、サンプル相関の補正が必要となる。また、十分に間隔を開けて収集したサンプルは
ほぼ独立のサンプルとみなすことができるので、例えば採取したサンプルのうち10000回に1個ずつを
パラメータの推定に用いれば相関による影響はないとみなすことができる。
　上記のような間隔をあけてサンプルを収集する場合は、rejectされた場合も1回とカウントする。また、このように収集したサンプルの繰り返し数をモンテカルロ・ステップ(MCS)と呼ぶ。

上記のメトロポリス法の例を使って、の平均を求めるプログラムを作成してみました↓

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

// 確率p(x1,x2,x3)を計算する
double cal_pro(
	int		x1,
	int		x2,
	int		x3,
	double	theta )
{
	double		numerator = 
				exp( theta * x1 * x2 + theta * x2 * x3 + theta * x3 * x1 );
	double		denominator = 
				2 * exp(3*theta) + 6 * exp( 0 - theta );
	double		ans =
				numerator / denominator;
	
	return ans;
}


int main( void )
{
	// x1, x2, x3の初期値
	int			x1 = 1;
	int			x2 = 1;
	int			x3 = -1;

	// 新しいx1, x2, x3
	int			nx1;
	int			nx2;
	int			nx3;

	// 繰り返し数
	int			rep = 3;
	// 最大反復数
	int			rep_max = 24000;
	// 固定パラメータ
	double		theta = 0.8;
	// 前回の確率
	double		pri_pro = cal_pro( x1, x2, x3, theta ); 
	// 次回の確率
	double		post_pro;
	// x1の値の格納域
	int			x1_list[rep_max/rep];
	// x1の平均値
	double		mean;

	int			i;
	int			j;
	int			cnt;
	int			int_ran;
	double		real_ran;
	double		r;

	for( i = 1; i <= rep_max; i ++ ) {
		// どの変数の値を変更するか決定する→ 候補
		int_ran = rand();
		if( ( int_ran % 3 ) == 0 ) {
			// x1を変更
			nx1 = 0 - x1;
			nx2 = x2;
			nx3 = x3;
		} else if( ( int_ran % 3 ) == 1 ) {
			// x2を変更
			nx1 = x1;
			nx2 = 0 - x2;
			nx3 = x3;
		} else {
			// x3を変更
			nx1 = x1;
			nx2 = x2;
			nx3 = 0 - x3;
		}

		// 候補の確率を求める
		post_pro = cal_pro( nx1, nx2, nx3, theta );

		// 現在と候補の確率の比を計算する
		r = post_pro / pri_pro;

		// 候補をacceptするかrejectするか判定する
		real_ran = (double)rand() / (double)( RAND_MAX + 1.0 );
		if( real_ran < r ) {
			// 候補をacceptする
			x1 = nx1;
			x2 = nx2;
			x3 = nx3;
		}
		// rejectされた場合はそのまま
		if( ( i % rep ) == 1 ) {
			x1_list[i/rep] = x1;	
			// x1の平均値を計算する
			mean = 0.0;
			cnt = 0;
			for( j = 0; j <= ( i / rep ); j ++ ) {
				mean += x1_list[j];
				cnt ++;
			}
			mean /= cnt; 
			printf( "%f\n", mean );
		}
		pri_pro = post_pro;
	}

	return 0;
}

これでθが0.2の場合と0.8の場合のの平均値の収束具合を比較してみました。上段が0.2で下段が
0.8です。これを見るとθが大きくなると収束が遅くなることがわかります。

正規乱数を生成するプログラムをＣ言語で書きました。
引数のチェックなどは省略してあります。

ちなみに、参考にしたサイトはここ。C++ですが、かなり丁寧に書いてあって、ものすごく参考になりました

あと、ここにも載っていました。こっちはCです。

 

---

 

#include
#include
#include

const int PI = 3.141592653589793

static int BOX_MULLER = 0;
static double X_NORMAL_RAN = 0.0;
static double Y_NORMAL_RAN = 0.0;

double normalRandom( double mean, double var )
{
 int   i;
 double  u1 = 0.0;
 double  u2 = 0.0;

 if( BOX_MULLER == 0 ) {
  BOX_MULLER = 1;
  u1 = (double)rand() / (double)RAND_MAX;
  while( u1 == 0.0 ) {
   u1 = (double)rand() / (double)RAND_MAX;
  }
  u2 = (double)rand() / (double)RAND_MAX;
  X_NORMAL_RAN = sqrt( -2.0 * log( u1 ) ) * cos( 2.0 * PI * u2 );
  Y_NORMAL_RAN = sqrt( -2.0 * log( u1 ) ) * sin( 2.0 * PI * u2 );
  return X_NORMAL_RAN * sqrt( var ) + mean;
 
 } else {
  BOX_MULLER = 0;
  return Y=NORMAL_RAN * sqrt( var ) + mean;
 }
} 
   

---

 

---