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

[Dirk Eddelbuettel]

On 16 April 2012 at 12:43, Marie Auger-Methe wrote:
| 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 

Super!

| 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?

I am still travelling and have not had a chance to look at this, but I
suspect that maybe something odd happen with PutRNGState / GetRNGState on the
side of rvm (and RNGScope deals with this for us).  There is no reason why it
shouldn't work, maybe we can work out a fix to be passed to CircStats. Else
we can always provide a suitable rvm in Rcpp too I suppose...

Dirk

 
| 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
| >
| 

-- 
R/Finance 2012 Conference on May 11 and 12, 2012 at UIC in Chicago, IL
See agenda, registration details and more at http://www.RinFinance.com