c - 当我调用 C 代码时 R 卡住

Question

我写了一个小的 C 代码来做随机游走大都市，我在 R 中调用它。当我运行它时，R 卡住。我不确定代码的哪一部分不正确。我关注这个Peng and Leeuw tutorial (第 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)。

Answer 1

您需要按照 Writing R Extensions 中的说明进行操作, 第 6.3 Random number generation并在调用 R 的随机数生成例程之前调用 GetRNGstate();。您还需要在完成后调用 PutRNGstate();。

您的代码开始工作的原因可能是因为您在调用 mcmc C 函数之前在 R session 中调用了 set.seed。

所以你的 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();
}