[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