Merge pull request #16 from osorensen/rng-fix

Native R and Rcpp random number generation

Author: osorensen

R/compute_sequentially.R

@@ -90,3 +90,4 @@ compute_sequentially <- function(
 
   ret <- run_smc(input_timeseries, hyperparameters, smc_options)
 }
+

src/particle.cpp

@@ -11,11 +11,13 @@ StaticParameters::StaticParameters(const vec& alpha, const umat& rho, const vec&
   alpha { alpha }, rho { rho }, tau { tau } {}
 
 StaticParameters::StaticParameters(const Prior& prior) :
-  alpha { randg(prior.n_clusters, distr_param(prior.alpha_shape, 1 / prior.alpha_rate)) },
+  alpha { Rcpp::rgamma(prior.n_clusters, prior.alpha_shape, 1 / prior.alpha_rate) },
   rho { umat(prior.n_items, prior.n_clusters) },
-  tau { normalise(randg(prior.n_clusters, distr_param(prior.cluster_concentration, 1)), 1) }
+  tau { normalise(Rcpp::as<vec>(Rcpp::rgamma(prior.n_clusters, prior.cluster_concentration, 1)), 1) }
   {
-    rho.each_col([&prior](uvec& a){ a = shuffle(regspace<uvec>(1, prior.n_items)); });
+    rho.each_col([&prior](uvec& a){
+      a = Rcpp::as<uvec>(Rcpp::sample(prior.n_items, prior.n_items, false));
+      });
   }
 
 Particle::Particle(const Options& options, const StaticParameters& parameters,

src/rejuvenate.cpp

@@ -7,7 +7,7 @@
 using namespace arma;
 
 uvec leap_and_shift(const uvec& current_rho, unsigned int cluster, const Prior& prior) {
-  unsigned int u = randi(distr_param(0, prior.n_items - 1));
+  unsigned int u = Rcpp::sample(prior.n_items, 1, false)[0] - 1;
   uvec intermediate_rho = current_rho;
   uvec rho_proposal = intermediate_rho;
   int rho_u = intermediate_rho(u);
@@ -16,7 +16,7 @@ uvec leap_and_shift(const uvec& current_rho, unsigned int cluster, const Prior&
     uvec(rho_u)
   );
 
-  unsigned int index = randi(distr_param(0, support.size() - 1));
+  unsigned int index = Rcpp::sample(support.size(), 1, false)[0] - 1;
   intermediate_rho(u) = support(index);
   for(size_t i{}; i < intermediate_rho.size(); i++) {
     if(current_rho(i) == current_rho(u)) {
@@ -85,7 +85,7 @@ bool Particle::rejuvenate(
   int proposed_particle_filter = Rcpp::sample(probs.size(), 1, false, probs, false)[0];
 
   bool accepted{};
-  if(log_ratio > log(randu())) {
+  if(log_ratio > log(R::runif(0, 1))) {
     this->parameters = StaticParameters{alpha_proposal, rho_proposal, parameters.tau};
     this->conditioned_particle_filter = proposed_particle_filter;
     this->log_incremental_likelihood = proposal_particle.log_incremental_likelihood;

src/resampler.cpp

@@ -26,7 +26,7 @@ ivec count_between_intervals(const vec& cumprob, const vec& u) {
 
 ivec stratsys(int n_samples, vec probs, bool stratified) {
   vec u(n_samples);
-  vec rn = stratified ? randu(n_samples) : vec(n_samples, fill::value(randu()));
+  vec rn = stratified ? Rcpp::runif(n_samples) : vec(n_samples, fill::value(R::runif(0, 1)));
 
   for(size_t i{}; i < n_samples; i++) u(i) = (i + rn(i)) / n_samples;
   return count_between_intervals(cumsum(probs), u);

src/sample_latent_rankings.cpp

@@ -8,6 +8,11 @@
 using namespace arma;
 namespace fs = std::filesystem;
 
+uvec shuffle_rcpp(const uvec& values_in) {
+  ivec inds = Rcpp::sample(values_in.size(), values_in.size(), false) - 1;
+  return values_in(conv_to<uvec>::from(inds));
+}
+
 LatentRankingProposal sample_latent_rankings(
     const std::unique_ptr<Data>& data, unsigned int t, const Prior& prior,
     std::string latent_rank_proposal,
@@ -45,7 +50,7 @@ LatentRankingProposal sample_latent_rankings(
       uvec tmp = ndit->second.observation;
 
       if(latent_rank_proposal == "uniform") {
-        tmp(ndit->second.available_items) = shuffle(ndit->second.available_rankings);
+        tmp(ndit->second.available_items) = shuffle_rcpp(ndit->second.available_rankings);
         proposal.proposal = join_horiz(proposal.proposal, tmp);
         proposal.log_probability = join_vert(
           proposal.log_probability, vec{-lgamma(ndit->second.available_rankings.size() + 1.0)});
@@ -55,7 +60,7 @@ LatentRankingProposal sample_latent_rankings(
         }
         double logprob{0};
 
-        uvec available_items_shuffled = shuffle(ndit->second.available_items);
+        uvec available_items_shuffled = shuffle_rcpp(ndit->second.available_items);
         uvec available_rankings = ndit->second.available_rankings;
 
         while(available_items_shuffled.size() > 1) {
@@ -135,7 +140,7 @@ LatentRankingProposal sample_latent_rankings(
       Rcpp::stop("No files.");
     }
 
-    int random_index = randi(distr_param(0, file_count - 1));
+    int random_index = Rcpp::sample(file_count, 1, false)[0] - 1;
 
     fs::path random_file_path;
     size_t current_index = 0;

tests/testthat/test-compute_sequentially_complete.R

@@ -40,5 +40,5 @@ test_that("compute_sequentially works with complete data", {
   )
   alpha_hat <- weighted.mean(x = as.numeric(mod$alpha), w = mod$importance_weights)
   expect_gt(alpha_hat, .99)
-  expect_lt(alpha_hat, 1.05)
+  expect_lt(alpha_hat, 1.06)
 })

tests/testthat/test-compute_sequentially_mixtures.R

@@ -1,4 +1,4 @@
-test_that("s", {
+test_that("Mixture models work", {
   set.seed(2)
   mod <- compute_sequentially(
     mixtures[1:50, ],
@@ -19,10 +19,10 @@ test_that("s", {
     tau[, i] <- tau[perm$permutations[i, ], i]
   }
 
-  expect_gt(weighted.mean(alpha[1, ], mod$importance_weights), 2)
-  expect_lt(weighted.mean(alpha[1, ], mod$importance_weights), 3)
-  expect_gt(weighted.mean(alpha[2, ], mod$importance_weights), .85)
-  expect_lt(weighted.mean(alpha[2, ], mod$importance_weights), 1.15)
+  expect_gt(weighted.mean(alpha[1, ], mod$importance_weights), 1)
+  expect_lt(weighted.mean(alpha[1, ], mod$importance_weights), 1.1)
+  expect_gt(weighted.mean(alpha[2, ], mod$importance_weights), 2)
+  expect_lt(weighted.mean(alpha[2, ], mod$importance_weights), 2.5)
 
   expect_gt(weighted.mean(tau[1, ], mod$importance_weights), .4)
   expect_lt(weighted.mean(tau[1, ], mod$importance_weights), .6)