<html><head><meta http-equiv="Content-Type" content="text/html; charset=UTF-8" /></head><body style='font-size: 10pt; font-family: Verdana,Geneva,sans-serif'>
<div style="font-size: 10pt; font-family: Verdana,Geneva,sans-serif;">
<p>Hi,</p>
<p>I have two different results within a for loop if a random draw is followed by an additional operation, i.e.:</p>
<p>// [[Rcpp::export]]<br />NumericMatrix define_e(int iter, NumericMatrix obs_mat, NumericVector draws) {<br /> Rcpp::NumericVector d = obs_mat(_, 1);<br /> Rcpp::NumericVector pr = no_init(d.length());<br /> Rcpp::NumericMatrix e(d.length(), iter);<br /> for (int i = 0; i < iter; i++) {<br /> pr = obs_mat(_, 0) * draws(i);<br /> e(_, i) = cpprbinom(d.length(), 1, pr);<br /> }<br /> return e;<br />}<br /><br />// [[Rcpp::export]]<br />NumericMatrix define_e2(int iter, NumericMatrix obs_mat, NumericVector draws) {<br /> Environment pkg = Environment::namespace_env("fastglm");<br /> Function f = pkg["fastglm"];<br /> Function s = pkg["summary.fastglm"];<br /><br /> Rcpp::NumericVector d = obs_mat(_, 1);<br /> Rcpp::NumericVector pr = no_init(d.length());<br /> Rcpp::NumericVector e = no_init(d.length());<br /> Rcpp::NumericMatrix e_mat(d.length(), iter);<br /><br /> Rcpp::NumericMatrix ematrix(d.length(), 2);<br /> Rcpp::NumericVector v(d.length(), 1);<br /> ematrix(_, 0) = v;<br /> Rcpp::List mod_pois;<br /> Rcpp::NumericVector mod_coef = no_init(d.length());<br /> double coef;<br /> for (int i = 0; i < iter; i++) {<br /> pr = obs_mat(_, 0) * draws(i);<br /> e = cpprbinom(d.length(), 1, pr);<br /> e_mat(_, i) = e; // save e as matrix for sake of comparison with define_e function<br /><br /> ematrix(_, 1) = e;<br /> if (all(is_na(e)).is_true()) {<br /> coef = NA_REAL;<br /> } else {<br /> mod_pois = f(Named("x", ematrix), Named("y", d), Named("family", "poisson"));<br /> mod_coef = mod_pois["coefficients"];<br /> coef = {mod_coef[1]};<br /> }<br /> }<br /> return e_mat;<br />}</p>
<p>define_e saves the random draws in a matrix (cpprbinom is from <a href="https://stackoverflow.com/questions/29430726/vectorised-rcpp-random-binomial-draws" target="_blank" rel="noopener noreferrer">https://stackoverflow.com/questions/29430726/vectorised-rcpp-random-binomial-draws</a>), 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,</p>
<p>Denis</p>
<p>> draws <- rbeta(200, 100, 60)<br />> obs_mat <- cbind(c(rep(1, 40), rep(0, 60)), c(rep(1, 20), rep(0, 80)))<br />> iter <- length(draws)<br />> set.seed(1234)<br />> t1 <- define_e(iter, obs_mat, draws)<br />> set.seed(1234)<br />> t2 <- define_e2(iter, obs_mat, draws)<br />> all.equal(t1, t2)<br />[1] "Mean relative difference: 2.270964"</p>
<div id="v1_rc_sig"></div>
</div>
</body></html>