3

我写了一个小的 C 代码来做随机游走大都会，我在 R 中调用它。当我运行它时，R 冻结。我不确定代码的哪一部分不正确。我遵循这个Peng and 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 回答