[Rcpp-devel] RNGScope, importing function, and if-else statement

Marie Auger-Methe marie.augermethe at gmail.com
Mon Apr 16 13:43:38 CEST 2012


Hi Dirk and list,

As suggested by Dirk I have re-written the function rvm in Rcpp (see 
below) and it now gives the same results as the rvm function in R. So my 
problem is fixed (although it would nice if you could just import an R 
RNG function in Rcpp... maybe a future development ;) ).
I'm still not sure what was the problem with my original method, my only 
clue is that the function rvm itself randomly generates multiple numbers 
using runif?

Thanks for all the help,

Marie

suppressMessages(library(CircStats))
suppressMessages(library(inline))
suppressMessages(library(sugaR))

rvmscr <- '
   // Inputs
   int nc = as<int>(n);
   double meanc = as<double>(mean);
   double kc = as<double>(k);

   // Assign memory
   Rcpp::NumericVector vm(nc);
   double U1;
   double U2;
   double U3;
   double z;
   double f;
   double c;

   double a = 1 + pow((1 + 4 * pow(kc,2)),0.5);
   double b = (a - pow(2*a,0.5))/(2*kc);
   double r = (1 + pow(b,2))/(2 * b);

   int obs = 0;

   RNGScope scope;

   while(obs < nc){
     U1 = as<double>(runif(1, 0, 1));
     z = cos(PI * U1);
     f = (1 + r * z)/(r + z);
     c = kc * (r - f);
     U2 = as<double>(runif(1, 0, 1));
     if(c * (2 - c) - U2 > 0){
       U3 = as<double>(runif(1, 0, 1));
         if(U3 < 0.5){
           vm[obs] = 2 * PI - acos(f) + meanc;
         }else{
           vm[obs] = acos(f) + meanc;
         }
       vm[obs] = fmod(vm[obs], 2*PI);
       obs += 1;
     }else{
       if(log(c/U2) + 1 - c >= 0){
         U3 = as<double>(runif(1, 0, 1));
         if(U3 < 0.5){
           vm[obs] = 2 * PI - acos(f) + meanc;
         }else{
           vm[obs] = acos(f) + meanc;
         }
         vm[obs] = fmod(vm[obs], 2*PI);
         obs += 1;
       }
     }
   }

   return vm;
'

rvmc <- cxxfunction(signature(n ="integer", mean = "numeric", k = 
"numeric"),
                     body = rvmscr, plugin = "Rcpp")

n <- 1000
mu <- 0
k <- 3
set.seed(89)
r1 <- rvmc(n,mu,k)
set.seed(89)
r2 <- rvm(n,mu,k)

max(r1 - r2)
cbind(r1[1:20],r2[1:20])

On 16/04/2012 10:29 AM, Dirk Eddelbuettel wrote:
> Marie,
>
> On 16 April 2012 at 08:56, Marie Auger-Methe wrote:
> | Hi Dirk,
> |
> | Thank you for the prompt reply!
> | By reimplement rvm using Rcpp sugar functions, do you mean re-write the
> | rvm function using sugar functions?
>
> Yes.
>
> (But as I said I have not yet looked at the code for rvm at all...)
>
> Dirk
>



More information about the Rcpp-devel mailing list