[Rcpp-devel] R’s random number generator not in itialized correctly when called at first

[Dirk Eddelbuettel]

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 math-library 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 (2011-10-31)
|         Copyright (C) 2011 The R Foundation for Statistical Computing
|         ISBN 3-900051-07-0
|         Platform: i486-pc-linux-gnu (32-bit)
|         
|         > 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/pgp-signature [Click mouse-2 to save to a file]
| 
| ----------------------------------------------------------------------
| _______________________________________________
| 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
-- 
"Outside of a dog, a book is a man's best friend. Inside of a dog, it is too
dark to read." -- Groucho Marx