[Rcpp-devel] for loop and random draws

[Denis Haine]

Sorry, I missed sharing that part comparing the two functions, vector 
vs. matrix without the regression.

// [[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_e1(int iter, NumericMatrix obs_mat, NumericVector 
draws) {
   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);

   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 demonstration
   }
   return e_mat;
}

So these two give the same output:

> 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)
> t0 <- define_e(iter, obs_mat, draws)
> set.seed(1234)
> t1 <- define_e1(iter, obs_mat, draws)
> all.equal(t0, t1)
[1] TRUE

But adding the regression to define_e1() (as define_e2() in previous 
email) does not give the same. The first iteration in the loop is 
identical in both, but not the subsequent ones. That's why I was 
thinking it could be linked to the seed, but (sorry for my limited 
knowledge of C++) it would be strange to manually "re-seed" at each 
iteration. Moreover, I believe the seed set in R is transferred to Rcpp 
and so does not need to be taken specifically care of.

Denis

Le 2024-10-01 11:14, Dirk Eddelbuettel a écrit :

> On 1 October 2024 at 09:59, Dirk Eddelbuettel wrote:
> | This is not reproducible. You show two functions, they both take 
> arguments,
> | you do not supply those. We cannot run this, so we cannot really help 
> you.
> 
> I missed the lines at the bottom past your signature, so I was wrong. 
> My bad.
> 
> I would still encourage you to _simplify_. You have a pair of 
> moderately
> complex functions, and they do not reproduce the same results after two 
> runs
> despite re-seeding. So I would trim elements from them one by one to 
> see
> which removal leads to identical results. Or, inversely, start from 
> twice
> drawing from an RNG after reseed and assert first that the draws are
> identical. Then add your processing.
> 
> Dirk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20241001/2c4cadd1/attachment.htm>