C 코드를 호출하면 R이 동결됩니다.

Question

랜덤 워크 메트로폴리스를 수행하기 위해 작은 C 코드를 작성했습니다.이 코드는 실행시 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) 

FOLLOW UP :는 while 루프는 문제가있는 곳. 그러나이 문제의 근본 원인은 lower bound와 upper bound 사이의 무작위 변수를 생성하는 함수 runif과 관련이있는 것 같습니다. 그러나 함수가 실제로하는 것은 무작위로 상한값 (5) 또는 하한값 (-5) 중 하나를 선택하는 것입니다.

Answer 1

Writing R Extensions, 6.3 Random number generation 섹션의 지침에 따라 R의 난수 생성 루틴을 호출하기 전에 GetRNGstate();으로 전화해야합니다. 완료되면 PutRNGstate();으로 전화해야합니다.

코드가 작동하기 시작한 이유는 mcmc C 함수를 호출하기 전에 R 세션에서 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(); 
}