[Rcpp-devel] for loop and random draws

[Denis Haine]

Hi,

I have two different results within a for loop if a random draw is 
followed by an additional operation, i.e.:

// [[Rcpp::export]]
NumericMatrix define_e(int iter, NumericMatrix obs_mat, NumericVector 
draws) {
   Rcpp::NumericVector d = obs_mat(_, 1);
   Rcpp::NumericVector pr = no_init(d.length());
   Rcpp::NumericMatrix e(d.length(), iter);
   for (int i = 0; i < iter; i++) {
     pr = obs_mat(_, 0) * draws(i);
     e(_, i) = cpprbinom(d.length(), 1, pr);
   }
   return e;
}

// [[Rcpp::export]]
NumericMatrix define_e2(int iter, NumericMatrix obs_mat, NumericVector 
draws) {
   Environment pkg = Environment::namespace_env("fastglm");
   Function f = pkg["fastglm"];
   Function s = pkg["summary.fastglm"];

   Rcpp::NumericVector d = obs_mat(_, 1);
   Rcpp::NumericVector pr = no_init(d.length());
   Rcpp::NumericVector e = no_init(d.length());
   Rcpp::NumericMatrix e_mat(d.length(), iter);

   Rcpp::NumericMatrix ematrix(d.length(), 2);
   Rcpp::NumericVector v(d.length(), 1);
   ematrix(_, 0) = v;
   Rcpp::List mod_pois;
   Rcpp::NumericVector mod_coef = no_init(d.length());
   double coef;
   for (int i = 0; i < iter; i++) {
     pr = obs_mat(_, 0) * draws(i);
     e = cpprbinom(d.length(), 1, pr);
     e_mat(_, i) = e; // save e as matrix for sake of comparison with 
define_e function

     ematrix(_, 1) = e;
     if (all(is_na(e)).is_true()) {
       coef = NA_REAL;
       } else {
       mod_pois = f(Named("x", ematrix), Named("y", d), Named("family", 
"poisson"));
       mod_coef = mod_pois["coefficients"];
       coef = {mod_coef[1]};
     }
   }
   return e_mat;
}

define_e saves the random draws in a matrix (cpprbinom is from 
https://stackoverflow.com/questions/29430726/vectorised-rcpp-random-binomial-draws), 
while define_e2 uses a vector to be provided to a regression model. 
Without the regression model in define_e2(), the two functions output 
the same values for e(), but not after. I absolutely have no clue why 
that is (maybe some "wiggling" with RNG seed??) and how to get the same 
output as with define_e but without saving as a matrix (as it can 
quickly become unmanageable for memory with N rows X number of 
iterations X 8k). Any help is welcome! Thanks,

Denis

> draws <- rbeta(200, 100, 60)
> obs_mat <- cbind(c(rep(1, 40), rep(0, 60)), c(rep(1, 20), rep(0, 80)))
> iter <- length(draws)
> set.seed(1234)
> t1 <- define_e(iter, obs_mat, draws)
> set.seed(1234)
> t2 <- define_e2(iter, obs_mat, draws)
> all.equal(t1, t2)
[1] "Mean relative difference: 2.270964"
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20241001/926b10ed/attachment.htm>