logo

nikq::cube

2007年1月25日(木) 23:54

メトロポリス法

メトロポリス法への第一ステップ。
いろんな次元の球の体積をMCとメトロポリスで求める。
コード
収束比較グラフ

一応きちんと同じ値に収束してるみたい。

普通にモンテカルロするときはこう。
accept = total  = 0;
for(i = 0; i < sample; i ++)
{
    accept += one_try( random_vector(dim) );  // 通った回数
    total  ++;                                // 総試行回数
}
rate = (double) accept / (double) total;
これをメトロポリス法に移植する。

ユーティリティ関数をまず用意する。
// 完全にランダムなベクトルで初期化
void make_seed( double *r, int dim ){
    for(int i=0;i// ちょっとだけずらす。たまにすごくずらす(エルゴード性確保のため)
void mutate_seed( double *r, double *r2, double p, int dim ){
    int i;
    if( genrand_real1() < 0.1 ) // large mutation
        make_seed( r2, dim );
    else{                       // small mutation
        for(i=0;i

その上で、実際のメトロポリスステップは以下のようになる。

make_seed( seed1 , dim );   // initial bootup
for(i = 0; i < sample; i ++){
    // mutate
    mutate_seed( seed1, seed2, 0.01, dim );
    contribute = one_try( seed2, dim );
    if( genrand_real1() < acceptance() ){
        accept_seed( seed2, seed1, dim );
        accept += contribute;
        total  ++;
    }
}
rate = (double) accept / (double) total;

動いた。

written by nikq [/program] [この記事のURL] [コメントを書く] [コメント(2)]

Comments

『間違えてたらごめんなさい』

突然すいません。
ちょっとメトロポリスに興味があったもので、ソース拝見させていただきましたが、

double acceptance( void ){
double a = T_yx()*f_y() / (T_xy()*f_x());
if( a < 1. )return a;
return 1.;
}

の部分で[return a;]される可能性はないきがするのですが?
実際にプログラムはしらせてみたところ、この[return a;]は実行されていなみたいです。

ですから、単純にすべてのサンプルがアクセプトされて、単純なモンテカルロしている事になっていたりするのでは・・・?

すごく頓珍漢なこといってたらすいません (ノ_・、)

written by あつし

『accept』

確かにそうなんですよね。
このサンプルでは超次元の球の体積を求めてるわけですが、
体積問題って密度一定な積分じゃないですか。
だから、メトロポリス法とモンテカルロ法に
あまり大きな差が出ないんじゃないかな。

範囲外rejectをmutateに追加すれば、
acceptance:return 1.;
でも問題ないです。多分。

それでも一応、
境界付近をランダムウォークしつつ動くので、
高次元の場合には単純モンテカルロよりも
多少サンプリング効率はいいはずです。

written by nikq

TrackBacks

nikq::cube

MySketch 2.7.2 written by 夕雨