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

| 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

Am Sonntag, den 04.12.2011, 20:10 -0600 schrieb Dirk Eddelbuettel:

On 5 December 2011 at 00:16, Paul Menzel wrote:

| 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.

Thanks!

|         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:

[?]

Indeed. Thank you very much.

| 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.

Understood.

Thank you very much again and I am sorry for not finding that out
myself. I took the post from [1] as an example and did not think about
the `set.seed(42)` call.

Now looking for `RNGScope` it is even in section 3.3 of the `Rcpp-FAQ`
(`vignette("Rcpp-FAQ")`) and as you wrote in the header file.


Thanks,

Paul


[1] http://lists.r-forge.r-project.org/pipermail/rcpp-devel/2011-September/002805.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 198 bytes
Desc: This is a digitally signed message part
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20111205/92254d4e/attachment.pgp>