[Rcpp-devel] RNGScope

Dirk Eddelbuettel edd at debian.org
Tue Aug 14 17:07:10 CEST 2012


Hi Rodney,

On 14 August 2012 at 09:17, Rodney Sparapani wrote:
| This is my first post so please be gentle ;o)  I have been through

Welcome -- the list is generally very nice and helpful. 

| all of the docs that I can find (but I haven't looked at the source
| code directly).  Here is the behavior that I see:
| 
|  > require(inline)
| Loading required package: inline
|  > source("helper.R")
|  > source("corrMH.R")
|  >
|  > set.seed(66)
|  >
|  > for(i in 1:3) print(corrMH(0.2, 0.5, 100, 0))
| [1] -0.07071726
| [1] -8.085111
| [1] -4.004675
|  >
|  > corrMH.Rcpp <- cxxfunction(
| +  signature(a="numeric", b="numeric", N="numeric", x="numeric"),
| +  body=cat.file("corrMH.Rcpp.cpp"), plugin="Rcpp")
|  >
|  > set.seed(66)
|  >
|  > for(i in 1:3) print(corrMH.Rcpp(0.2, 0.5, 100, 0))
| [1] -0.07071726
| [1] -0.07071726
| [1] -0.07071726
| 
| Notice that the Rcpp version of the function returns the same value
| each time (only 3 shown, but this is true no matter how many times
| that you do it).  I guess this is because we follow the RNGScope state
| of the RNG, but we don't update it.  However, updating is really needed
| in this case.  I hope I am making this clear.  Please let me know if
| you have any ideas.  Thanks

Can you please try to reduces this to a minimal example which does not
pull in a number of other files?  

All the documentation says (SHOUTS actually) to use RNGScope so that get/put
rng state are called.  

Minimally reproducible example I had hanging around in /tmp anyway:

    library(inline)

    f <- cxxfunction(signature(), plugin="Rcpp", body='
         Rcpp::RNGScope tmp;
         return rnorm(3);
    ')

Now if we build that:

    R> library(inline)
    R> f <- cxxfunction(signature(), plugin="Rcpp", body='
    +    Rcpp::RNGScope tmp;
    +    return rnorm(3);
    + ')
    R> 
    R> set.seed(42); f()
    [1]  1.370958 -0.564698  0.363128
    R> set.seed(42); rnorm(3)
    [1]  1.370958 -0.564698  0.363128
    R> 

we saw that it perfectly matches the R draws.  If your boil your code down to
such a bare minimum you should get the same behaviour irresptive what
distrbution you pull from.

Hth, Dirk


| 
| Attachments seem to be allowed so I'm including my code...
| 
| Red Hat Enterprise Linux Server release 6.2 (Santiago)
| g++ (GCC) 4.4.6 20110731 (Red Hat 4.4.6-3)
| R version 2.14.2 (2012-02-29)
| 
|  > installed.packages()["Rcpp", ]
|                       Package                      LibPath
|                        "Rcpp" "/opt/local/lib64/R/library"
|                       Version                     Priority
|                      "0.9.10"                           NA
|                       Depends                      Imports
|      "R (>= 2.12.0), methods"                           NA
|                     LinkingTo                     Suggests
|                            NA  "RUnit, inline, rbenchmark"
|                      Enhances                      OS_type
|                            NA                           NA
|                       License                        Built
|                  "GPL (>= 2)"                     "2.14.2"
| -- 
| Rodney Sparapani, PhD  Center for Patient Care and Outcomes Research
| Sr. Biostatistician               http://www.mcw.edu/pcor
| 4 wheels good, 2 wheels better!   Medical College of Wisconsin (MCW)
| WWLD?:  What Would Lombardi Do?   Milwaukee, WI, USA
| 
| ----------------------------------------------------------------------
| 
| Rcpp::NumericVector _a(a), _b(b), _N(N), _x(x),
|   _y=rcauchy(1, _x[0], pow(_a[0], -0.5)), _q=runif(1, 0., 1.);
| 
| RNGScope scope; // use the same RNG state that R uses
| 
| double x2=pow(_x[0], 2.), y2=pow(_y[0], 2.), 
|   p=pow((1.+y2)/(1.+x2), 0.5*_N[0]-1.5) *
|       exp(_b[0]*(_y[0]*sqrt(1.+y2)-_x[0]*sqrt(1.+x2))-0.5*_a[0]*(y2-x2));
| 
| if (_q[0]<p) return(_y);
| else return(x);
| 
| ----------------------------------------------------------------------
| #library(MASS)
| library(compiler)
| 
| # returns rho/sqrt(1-rho^2) rather than rho
| # a useful transformation as discussed in
| # Johnson, Klotz and Balakrishnan
| # Continuous Univariate Distributions
| # x, which is the current value, is assumed
| # to have come from the same transformation
| 
| corrMH <- cmpfun(function(a, b, N, x) {
|   p <- NaN
| 
|   while (is.nan(p)) {
|     y <- rcauchy(1, location=x, scale=a^(-0.5))
| 
|     p <- ((1+y^2)/(1+x^2))^(0.5*N-1.5) *
|       exp(b*(y*sqrt(1+y^2)-x*sqrt(1+x^2))-0.5*a*(y^2-x^2))
|   }
|          
|   if (runif(1)<p) return(y)
|   else return(x)
| })
| 
| 
| ----------------------------------------------------------------------
| 
| cat.char <- function(char) {
|   # arrays of strings are very annoying
|   # cat.char converts an array of strings into one long string with line-feeds
| 
|   cum <- NULL
|   
|   for(i in 1:length(char)) {
|     if(i==1) cum <- sprintf("%s\n", char[1])
|     else cum <- sprintf("%s%s\n", cum, char[i])
|   }
|   
|   return(cum)
| }
| 
| cat.file <- function(file) {
|   # cat.file returns the contents of a file as a string
|   
|   content <- NULL
|   f1 <- file.info(file)
|   
|   if(!is.na(f1$size) & !f1$isdir & f1$size>0) {
|     f1 <- file(file, open="r")
| 
|     if(isOpen(f1)) {
|       content <- readLines(f1)
| 
|       close(f1)
|     }
|   }
| 
|   return(cat.char(content))
| }
| 
| ## open.file <- function(file) {
|   
| ##   content <- NULL
| ##   f1 <- NULL
|   
| ##   f1 <- file(file, open="r")
| 
| ##   if(isOpen(f1)) {
| ##     content <- readLines(f1)
| 
| ##     close(f1)
| ##   }
| 
| ##   return(content)
| ## }
| 
| 
| ----------------------------------------------------------------------
| require(inline)
| source("helper.R")
| source("corrMH.R")
| 
| set.seed(66)
| 
| for(i in 1:30) print(corrMH(0.2, 0.5, 100, 0))
| 
| corrMH.Rcpp <- cxxfunction(
|  signature(a="numeric", b="numeric", N="numeric", x="numeric"), 
|  body=cat.file("corrMH.Rcpp.cpp"), plugin="Rcpp")
| 
| set.seed(66)
| 
| for(i in 1:30) print(corrMH.Rcpp(0.2, 0.5, 100, 0))
| 
| ----------------------------------------------------------------------
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com  


More information about the Rcpp-devel mailing list