This is my first post so please be gentle ;o) I have been through
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
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@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
