[Rcppdevel] R’s random number generator not in itialized correctly when called at first
Dirk Eddelbuettel
edd at debian.org
Mon Dec 5 03:10:47 CET 2011
On 5 December 2011 at 00:16, Paul Menzel wrote:
 Dear Rcpp folks,


 I am running an up to date Debian Sid/unstable system and I am
 simulating how long a random walks stays negative and I am using Rcpp
 and inline.

 library(inline)
 inc < "
 #include <stdio.h>
 #include <stdlib.h>
 "
As an aside, that 'inc' snippet is not needed.
 rwRcpp <cxxfunction(signature(ns="numeric", reps="numeric"),
 includes=inc, plugin="Rcpp", body='
 double n = Rcpp::as<double>(ns);
 double rep = Rcpp::as<double>(reps);

 double x = 0;

 Rcpp::NumericVector s(n + 2);

 for (int i = 0; i < rep; i++) {
 double len = 0;
 double csum = 0;

 do {
Your mistake: no use of RNGScope. Quickly google the docs, and/or the
mailing list archives. In short, somewhere in the file add a line
RNGScope scp;
(or some other variable name instead of scp). Also read 'Writing R
Extensions' on the standalone mathlibrary and the need for keeping state.
It simply is your error, and we plastered the header files with notes about
how you MUSt add this.
Once added, all is well:
edd at max:/tmp$ cat paul.r
library(inline)
rwRcpp <cxxfunction(signature(ns="numeric", reps="numeric"),
plugin="Rcpp", body='
double n = Rcpp::as<double>(ns);
double rep = Rcpp::as<double>(reps);
double x = 0;
Rcpp::NumericVector s(n + 2);
RNGScope scp;
for (int i = 0; i < rep; i++) {
double len = 0;
double csum = 0;
do {
x = rnorm(1, 0, 1)[0];
Rprintf("x = %f\\n", x);
csum = csum + x;
} while ((csum < 0) && (len++ < n));
s[len]++;
}
return s;
')
print(rwRcpp(2,5))
print(rwRcpp(2,5))
edd at max:/tmp$ r paul.r
Loading required package: methods
x = 0.938212
x = 0.007566
x = 0.483957
x = 0.523709
x = 2.352028
x = 0.136084
x = 0.201972
[1] 4 0 0 1
x = 0.767253
x = 0.168308
x = 1.270389
x = 0.255534
x = 1.427180
x = 1.333946
x = 1.474912
x = 0.035768
x = 1.124917
x = 0.863675
x = 0.065182
x = 0.518603
x = 0.163492
[1] 1 0 1 3
edd at max:/tmp$
 x = rnorm(1, 0, 1)[0];
 Rprintf("x = %f\\n", x);

 csum = csum + x;
 } while ((csum < 0) && (len++ < n));

 s[len]++;
 }

 return s;
 ')

 Now the interesting thing is that the random numbers do not seem to be
 random when the function is called at the very beginning.

 $ R
 R version 2.14.0 (20111031)
 Copyright (C) 2011 The R Foundation for Statistical Computing
 ISBN 3900051070
 Platform: i486pclinuxgnu (32bit)

 > source("rw.r")
 > rwRcpp(4,10)
 x = 8.773321
 […]
 x = 8.773321
 [1] 0 0 0 0 0 10
 > rwRcpp(2,5)
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 [1] 0 0 0 5

 Setting the seed value and it starts working.

 > set.seed(1)
 > rwRcpp(2,5)
 x = 0.626454
 x = 0.183643
 x = 0.835629
 x = 1.595281
 x = 0.329508
 x = 0.820468
 x = 0.487429
 x = 0.738325
 x = 0.575781
 [1] 3 0 1 1

 Opening a now R session, I tested also the following.

 > source("rw.r")
 > rwRcpp(2,5)
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 x = 8.773321
 [1] 0 0 0 5
 > rnorm(2)
 [1] 0.9520938 2.0510458
 > rwRcpp(2,5)
 x = 0.438619
 x = 0.606704
 x = 0.658996
 x = 0.082486
 x = 2.124667
 x = 0.127444
 x = 0.447259
 x = 0.363092
 x = 0.822970
 x = 0.387696
 x = 1.454236
 [1] 1 2 0 2

 So it looks like, Rcpp is not initializing(?) the random number
 generator in R? Could that be the cause?
You called it incorrectly, but there is literally no way for us to prevent
this, or else we'd put such a safeguard up.
Dirk

 Thanks,

 Paul

 
 library(inline)
 inc < "
 #include <stdio.h>
 #include <stdlib.h>
 "

 rwRcpp <cxxfunction(signature(ns="numeric", reps="numeric"), includes=inc, plugin="Rcpp", body='
 double n = Rcpp::as<double>(ns);
 double rep = Rcpp::as<double>(reps);

 double x = 0;

 Rcpp::NumericVector s(n + 2);

 for (int i = 0; i < rep; i++) {
 double len = 0;
 double csum = 0;

 do {
 x = rnorm(1, 0, 1)[0];
 Rprintf("x = %f\\n", x);

 csum = csum + x;
 } while ((csum < 0) && (len++ < n));

 s[len]++;
 }

 return s;
 ')
 xapplication/pgpsignature [Click mouse2 to save to a file]

 
 _______________________________________________
 Rcppdevel mailing list
 Rcppdevel at lists.rforge.rproject.org
 https://lists.rforge.rproject.org/cgibin/mailman/listinfo/rcppdevel

"Outside of a dog, a book is a man's best friend. Inside of a dog, it is too
dark to read."  Groucho Marx
More information about the Rcppdevel
mailing list