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"
_______________________________________________
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
