How to setup setseed for one's Rcpp function using RcppZiggurat samplers?

Tomasz 0
Viewed 35

I am developing some sampling algorithms based on RcppZiggurat. My objective is to use the samplers in cpp code but not having to come up with new ways of setting the seed but to rely on existing functions, such as RcppZiggurat::zsetseed.

This are my two cpp functions using the sampler zigg.norm():

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

#include <Ziggurat.h>
// [[Rcpp::depends(RcppZiggurat)]]

using namespace Rcpp;
using namespace arma;

static Ziggurat::Ziggurat::Ziggurat zigg;

// [[Rcpp::export]]
vec zzrnorm(int n) {
  vec x(n);
  for (int i=0; i<n; i++) {
    x(i) = zigg.norm();
  }
  return x;
}

// [[Rcpp::export]]
vec zzrnorm_ss(int n, unsigned long int set_seed) {
  zigg.setSeed(set_seed);
  vec x(n);
  for (int i=0; i<n; i++) {
    x(i) = zigg.norm();
  }
  return x;
}

And here is the R code and outputs that I am studying and cannot understand the outputs:

> Rcpp::sourceCpp("so_Ziggurat.cpp")
> RcppZiggurat::zsetseed(1)
> RcppZiggurat::zrnorm(3)
[1] -1.4808858  0.9991196 -0.4898021

> RcppZiggurat::zsetseed(1)
> RcppZiggurat::zrnorm(3)
[1] -1.4808858  0.9991196 -0.4898021

> RcppZiggurat::zsetseed(1)

> zzrnorm(3)
           [,1]
[1,] -1.1409054
[2,]  0.8759027
[3,] -0.1969075

> RcppZiggurat::zsetseed(1)

> zzrnorm(3)
           [,1]
[1,] -0.1245053
[2,]  0.2695581
[3,] -0.4185105

> zzrnorm_ss(3, 1)
           [,1]
[1,] -1.4808858
[2,]  0.9991196
[3,] -0.4898021

> zzrnorm_ss(3, 1)
           [,1]
[1,] -1.4808858
[2,]  0.9991196
[3,] -0.4898021

The first two generators are from the RcppZiggyrat and use RcppZiggurat::zsetseed to set seed and RcppZiggurat::zrnorm to generate random numbers. They are reproducible, of course.

The second set of outputs is something I don't understand. I use the package provided RcppZiggurat::zsetseed to set seed and my function zzrnorm to sample random numbers. My function uses the same class of samplers for which RcppZiggurat::zsetseed sets the seed based on

static Ziggurat::Ziggurat::Ziggurat zigg;

But RcppZiggurat::zsetseed does not set the seed, and the generated numbers differ despite the fact that RcppZiggurat::zsetseed uses the same source code: zigg.setSeed(set_seed); to set the seed.

Now, if I include RcppZiggurat::zsetseed within my function as in the third example then the seed is set properly and the results are reproducible.

So, why isn't that the case in the second set of outputs?

I would greatly appreciate some suggestions as this would be my preferred way of working with the samplers in R.

rcpp r c++
1 Answers

Nice puzzle. I think the key is that your two calls in the middle call your call local instance which has been initialized with the default seed. Witness an Rscript call on the command-line with defaults getting the same six values:

$ Rscript -e 'RcppZiggurat::zrnorm(6)'       
[1] -1.140905  0.875903 -0.196907 -0.124505  0.269558 -0.418510  
$

So to make the calls comparable, you need to set seeds. A simpler version follow. Here we call the instance from the package twice, with seeding, and then your instance, also with seeding. Unsurprisingly, the results are the same.

Code

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]

#include <Ziggurat.h>
// [[Rcpp::depends(RcppZiggurat)]]

static Ziggurat::Ziggurat::Ziggurat zigg;

// [[Rcpp::export]]
arma::rowvec zzrnorm_ss(int n, unsigned long int set_seed) {
  zigg.setSeed(set_seed);
  arma::rowvec x(n);
  for (int i=0; i<n; i++) {
    x(i) = zigg.norm();
  }
  return x;
}

/*** R
RcppZiggurat::zsetseed(1)
v1 <- RcppZiggurat::zrnorm(3)
RcppZiggurat::zsetseed(1)
v2 <- RcppZiggurat::zrnorm(3)
v3 <- zzrnorm_ss(3, 1)
v4 <- zzrnorm_ss(3, 1)
rbind(v1, v2, v3, v4)
*/

Output

> Rcpp::sourceCpp("~/git/stackoverflow/73679283/answer.cpp")

> RcppZiggurat::zsetseed(1)

> v1 <- RcppZiggurat::zrnorm(3)

> RcppZiggurat::zsetseed(1)

> v2 <- RcppZiggurat::zrnorm(3)

> v3 <- zzrnorm_ss(3, 1)

> v4 <- zzrnorm_ss(3, 1)

> rbind(v1, v2, v3, v4)
       [,1]    [,2]      [,3]
v1 -1.48089 0.99912 -0.489802
v2 -1.48089 0.99912 -0.489802
   -1.48089 0.99912 -0.489802
   -1.48089 0.99912 -0.489802
>
Dirk Eddelbuettel0
Related