3

ランダム ウォーク メトロポリスを実行する小さな C コードを書きました。これを R で呼び出します。実行すると、R がフリーズします。コードのどの部分が間違っているかわかりません。このPeng と Leeuw のチュートリアル(6 ページ) に従っています。免責事項として: 私は C の経験があまりなく、C++ の基本的な知識しかありません。

#----C code --------
#include <R.h>
#include <Rmath.h>

void mcmc(int *niter, double *mean, double *sd, double *lo_bound, 
          double *hi_bound, double *normal)
{
    int i, j;
    double x, x1, h, p;
    x = runif(-5, 5);
    for(i=0; i < *niter; i++) {
        x1 = runif(*lo_bound, *hi_bound);
        while((x1 + x) > 5 || (x1 + x) < -5)
            x1 = runif(*lo_bound, *hi_bound);
        h = dnorm(x+x1, *mean, *sd, 0)/dnorm(x, *mean, *sd, 0);
        if(h >= 1)
            h = 1;
        p = runif(0, 1);
        if(p < h)
            x += x1;
        normal[i] = x;
    }
}


#-----R code ---------
foo_C<-function(mean, sd, lo_bound, hi_bound, niter)
{
    result <- .C("mcmc",  as.integer(niter), as.double(mean), as.double(sd), 
                 as.double(lo_bound), as.double(hi_bound), normal=double(niter))
    result[["normal"]]
}

それをコンパイルした後：

dyn.load("foo_C.so")
foo_C(0, 1, -0.5, 0.5, 100)

フォローアップ: ループwhileに問題があります。しかし、問題の根源は関数 に関係しているようですrunif。これは、下限と上限の間の確率変数を生成することになっています。しかし、関数が実際に行うことは、上限値 (5) または下限値 (-5) のいずれかをランダムに選択することのようです。

4

1 に答える 1

4

R の乱数生成ルーチンを呼び出す前に、 R 拡張機能の記述、セクション6.3 乱数の生成と呼び出しの指示に従う必要があります。GetRNGstate();また、終了時に電話する必要がありPutRNGstate();ます。

コードが機能し始めた理由は、 C 関数set.seedを呼び出す前に R セッションで呼び出したことが原因である可能性があります。mcmc

したがって、C コードは次のようになります。

#include <R.h>
#include <Rmath.h>

void mcmc(int *niter, double *mean, double *sd, double *lo_bound, 
          double *hi_bound, double *normal)
{
    int i;
    double x, x1, h, p;
    GetRNGstate();
    x = runif(-5.0, 5.0);
    for(i=0; i < *niter; i++) {
        x1 = runif(*lo_bound, *hi_bound);
        while((x1 + x) > 5.0 || (x1 + x) < -5.0) {
            x1 = runif(*lo_bound, *hi_bound);
            //R_CheckUserInterrupt();
        }
        h = dnorm(x+x1, *mean, *sd, 0)/dnorm(x, *mean, *sd, 0);
        if(h >= 1)
            h = 1;
        p = runif(0, 1);
        if(p < h)
            x += x1;
        normal[i] = x;
    }
    PutRNGstate();
}
于 2013-04-09T12:14:49.330 に答える