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
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
