diff --git a/src/stan/model/model_header.hpp b/src/stan/model/model_header.hpp
index 75ac6c19cbd..34948377f74 100644
--- a/src/stan/model/model_header.hpp
+++ b/src/stan/model/model_header.hpp
@@ -13,9 +13,6 @@
 #include <stan/model/indexing.hpp>
 #include <stan/services/util/create_rng.hpp>
 
-#include <boost/random/additive_combine.hpp>
-#include <boost/random/linear_congruential.hpp>
-
 #include <cmath>
 #include <cstddef>
 #include <fstream>
diff --git a/src/stan/services/experimental/advi/fullrank.hpp b/src/stan/services/experimental/advi/fullrank.hpp
index 5c651773ded..96cf5d05c0b 100644
--- a/src/stan/services/experimental/advi/fullrank.hpp
+++ b/src/stan/services/experimental/advi/fullrank.hpp
@@ -10,7 +10,6 @@
 #include <stan/services/error_codes.hpp>
 #include <stan/io/var_context.hpp>
 #include <stan/variational/advi.hpp>
-#include <boost/random/additive_combine.hpp>
 #include <string>
 #include <vector>
 
diff --git a/src/stan/services/experimental/advi/meanfield.hpp b/src/stan/services/experimental/advi/meanfield.hpp
index a89aa521cd7..6cffe548acf 100644
--- a/src/stan/services/experimental/advi/meanfield.hpp
+++ b/src/stan/services/experimental/advi/meanfield.hpp
@@ -10,7 +10,6 @@
 #include <stan/services/error_codes.hpp>
 #include <stan/io/var_context.hpp>
 #include <stan/variational/advi.hpp>
-#include <boost/random/additive_combine.hpp>
 #include <string>
 #include <vector>
 
diff --git a/src/test/unit/io/random_var_context_test.cpp b/src/test/unit/io/random_var_context_test.cpp
index 18bfeb7f2b0..64ea370f17a 100644
--- a/src/test/unit/io/random_var_context_test.cpp
+++ b/src/test/unit/io/random_var_context_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/io/random_var_context.hpp>
 #include <stan/io/empty_var_context.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
 #include <boost/random/uniform_real_distribution.hpp>
 #include <test/test-models/good/services/test_lp.hpp>
 #include <test/unit/util.hpp>
@@ -95,12 +95,12 @@ class random_var_context : public testing::Test {
   random_var_context()
       : empty_context(),
         model(empty_context, 0, static_cast<std::stringstream*>(0)),
-        rng(0),
+        rng(stan::services::util::create_rng(0, 0)),
         throwing_model() {}
 
   stan::io::empty_var_context empty_context;
   stan_model model;
-  boost::ecuyer1988 rng;
+  stan::rng_t rng;
   test::mock_throwing_model_in_write_array throwing_model;
 };
 
diff --git a/src/test/unit/mcmc/chains_test.cpp b/src/test/unit/mcmc/chains_test.cpp
index 987279de66e..0862b962215 100644
--- a/src/test/unit/mcmc/chains_test.cpp
+++ b/src/test/unit/mcmc/chains_test.cpp
@@ -1,7 +1,6 @@
 #include <stan/mcmc/chains.hpp>
 #include <stan/io/stan_csv_reader.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/additive_combine.hpp>
 #include <set>
 #include <exception>
 #include <utility>
diff --git a/src/test/unit/mcmc/hmc/base_hmc_test.cpp b/src/test/unit/mcmc/hmc/base_hmc_test.cpp
index 00742532075..ce3d8f51f45 100644
--- a/src/test/unit/mcmc/hmc/base_hmc_test.cpp
+++ b/src/test/unit/mcmc/hmc/base_hmc_test.cpp
@@ -1,22 +1,20 @@
 #include <test/unit/mcmc/hmc/mock_hmc.hpp>
 #include <stan/callbacks/stream_writer.hpp>
 #include <stan/mcmc/hmc/base_hmc.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <boost/algorithm/string/split.hpp>
 #include <test/unit/util.hpp>
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 namespace stan {
 namespace mcmc {
 
-class mock_hmc
-    : public base_hmc<mock_model, mock_hamiltonian, mock_integrator, rng_t> {
+class mock_hmc : public base_hmc<mock_model, mock_hamiltonian, mock_integrator,
+                                 stan::rng_t> {
  public:
-  mock_hmc(const mock_model& m, rng_t& rng)
-      : base_hmc<mock_model, mock_hamiltonian, mock_integrator, rng_t>(m, rng) {
-  }
+  mock_hmc(const mock_model& m, stan::rng_t& rng)
+      : base_hmc<mock_model, mock_hamiltonian, mock_integrator, stan::rng_t>(
+          m, rng) {}
 
   sample transition(sample& init_sample, callbacks::logger& logger) {
     this->seed(init_sample.cont_params());
@@ -33,7 +31,7 @@ class mock_hmc
 }  // namespace stan
 
 TEST(McmcBaseHMC, point_construction) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -47,14 +45,14 @@ TEST(McmcBaseHMC, point_construction) {
 }
 
 TEST(McmcBaseHMC, point_access_from_const_hmc) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
   const stan::mcmc::mock_hmc sampler(stan::mcmc::mock_model(2), base_rng);
   EXPECT_EQ(2, sampler.z().q.size());
   EXPECT_EQ(2, sampler.z().g.size());
 }
 
 TEST(McmcBaseHMC, seed) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -70,7 +68,7 @@ TEST(McmcBaseHMC, seed) {
 }
 
 TEST(McmcBaseHMC, set_nominal_stepsize) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -88,7 +86,7 @@ TEST(McmcBaseHMC, set_nominal_stepsize) {
 }
 
 TEST(McmcBaseHMC, set_stepsize_jitter) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -108,7 +106,7 @@ TEST(McmcBaseHMC, set_stepsize_jitter) {
 TEST(McmcBaseHMC, streams) {
   stan::test::capture_std_streams();
 
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
diff --git a/src/test/unit/mcmc/hmc/hamiltonians/base_hamiltonian_test.cpp b/src/test/unit/mcmc/hmc/hamiltonians/base_hamiltonian_test.cpp
index 4926f6daef0..41537c282c0 100644
--- a/src/test/unit/mcmc/hmc/hamiltonians/base_hamiltonian_test.cpp
+++ b/src/test/unit/mcmc/hmc/hamiltonians/base_hamiltonian_test.cpp
@@ -3,12 +3,10 @@
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <stan/callbacks/stream_writer.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <test/unit/util.hpp>
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(BaseHamiltonian, update_potential_gradient) {
   stan::io::empty_var_context data_var_context;
 
@@ -16,7 +14,7 @@ TEST(BaseHamiltonian, update_potential_gradient) {
   funnel_model_namespace::funnel_model model(data_var_context, 0,
                                              &model_output);
 
-  stan::mcmc::mock_hamiltonian<funnel_model_namespace::funnel_model, rng_t>
+  stan::mcmc::mock_hamiltonian<funnel_model_namespace::funnel_model, stan::rng_t>
       metric(model);
   stan::mcmc::ps_point z(11);
   z.q.setOnes();
diff --git a/src/test/unit/mcmc/hmc/hamiltonians/dense_e_metric_test.cpp b/src/test/unit/mcmc/hmc/hamiltonians/dense_e_metric_test.cpp
index 31f5bcbcd75..9f78a430246 100644
--- a/src/test/unit/mcmc/hmc/hamiltonians/dense_e_metric_test.cpp
+++ b/src/test/unit/mcmc/hmc/hamiltonians/dense_e_metric_test.cpp
@@ -1,5 +1,5 @@
 #include <string>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <test/unit/mcmc/hmc/mock_hmc.hpp>
 #include <stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp>
@@ -8,10 +8,8 @@
 #include <test/unit/util.hpp>
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcDenseEMetric, sample_p) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::Matrix2d m(2, 2);
   m(0, 0) = 3.0;
@@ -23,7 +21,7 @@ TEST(McmcDenseEMetric, sample_p) {
 
   stan::mcmc::mock_model model(2);
 
-  stan::mcmc::dense_e_metric<stan::mcmc::mock_model, rng_t> metric(model);
+  stan::mcmc::dense_e_metric<stan::mcmc::mock_model, stan::rng_t> metric(model);
   stan::mcmc::dense_e_point z(2);
   z.set_metric(m_inv);
 
@@ -62,7 +60,6 @@ TEST(McmcDenseEMetric, sample_p) {
 }
 
 TEST(McmcDenseEMetric, gradients) {
-  rng_t base_rng(0);
 
   Eigen::VectorXd q = Eigen::VectorXd::Ones(11);
 
@@ -79,7 +76,7 @@ TEST(McmcDenseEMetric, gradients) {
   funnel_model_namespace::funnel_model model(data_var_context, 0,
                                              &model_output);
 
-  stan::mcmc::dense_e_metric<funnel_model_namespace::funnel_model, rng_t>
+  stan::mcmc::dense_e_metric<funnel_model_namespace::funnel_model, stan::rng_t>
       metric(model);
 
   double epsilon = 1e-6;
@@ -156,7 +153,6 @@ TEST(McmcDenseEMetric, gradients) {
 TEST(McmcDenseEMetric, streams) {
   stan::test::capture_std_streams();
 
-  rng_t base_rng(0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -165,7 +161,8 @@ TEST(McmcDenseEMetric, streams) {
   stan::mcmc::mock_model model(q.size());
 
   // typedef to use within Google Test macros
-  typedef stan::mcmc::dense_e_metric<stan::mcmc::mock_model, rng_t> dense_e;
+  typedef stan::mcmc::dense_e_metric<stan::mcmc::mock_model, stan::rng_t>
+      dense_e;
 
   EXPECT_NO_THROW(dense_e metric(model));
 
diff --git a/src/test/unit/mcmc/hmc/hamiltonians/diag_e_metric_test.cpp b/src/test/unit/mcmc/hmc/hamiltonians/diag_e_metric_test.cpp
index 1d11e26db95..475afbf56d5 100644
--- a/src/test/unit/mcmc/hmc/hamiltonians/diag_e_metric_test.cpp
+++ b/src/test/unit/mcmc/hmc/hamiltonians/diag_e_metric_test.cpp
@@ -1,5 +1,5 @@
 #include <string>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <test/unit/mcmc/hmc/mock_hmc.hpp>
 #include <stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp>
@@ -8,17 +8,15 @@
 #include <test/unit/util.hpp>
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcDiagEMetric, sample_p) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
   q(1) = 1;
 
   stan::mcmc::mock_model model(q.size());
-  stan::mcmc::diag_e_metric<stan::mcmc::mock_model, rng_t> metric(model);
+  stan::mcmc::diag_e_metric<stan::mcmc::mock_model, stan::rng_t> metric(model);
   stan::mcmc::diag_e_point z(q.size());
 
   int n_samples = 1000;
@@ -44,7 +42,6 @@ TEST(McmcDiagEMetric, sample_p) {
 }
 
 TEST(McmcDiagEMetric, gradients) {
-  rng_t base_rng(0);
 
   Eigen::VectorXd q = Eigen::VectorXd::Ones(11);
 
@@ -61,8 +58,8 @@ TEST(McmcDiagEMetric, gradients) {
   std::stringstream debug, info, warn, error, fatal;
   stan::callbacks::stream_logger logger(debug, info, warn, error, fatal);
 
-  stan::mcmc::diag_e_metric<funnel_model_namespace::funnel_model, rng_t> metric(
-      model);
+  stan::mcmc::diag_e_metric<funnel_model_namespace::funnel_model, stan::rng_t>
+      metric(model);
 
   double epsilon = 1e-6;
 
@@ -137,7 +134,6 @@ TEST(McmcDiagEMetric, gradients) {
 TEST(McmcDiagEMetric, streams) {
   stan::test::capture_std_streams();
 
-  rng_t base_rng(0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -146,7 +142,7 @@ TEST(McmcDiagEMetric, streams) {
   stan::mcmc::mock_model model(q.size());
 
   // typedef to use within Google Test macros
-  typedef stan::mcmc::diag_e_metric<stan::mcmc::mock_model, rng_t> diag_e;
+  typedef stan::mcmc::diag_e_metric<stan::mcmc::mock_model, stan::rng_t> diag_e;
 
   EXPECT_NO_THROW(diag_e metric(model));
 
diff --git a/src/test/unit/mcmc/hmc/hamiltonians/softabs_metric_test.cpp b/src/test/unit/mcmc/hmc/hamiltonians/softabs_metric_test.cpp
index c1a477c1a1c..0804105101a 100644
--- a/src/test/unit/mcmc/hmc/hamiltonians/softabs_metric_test.cpp
+++ b/src/test/unit/mcmc/hmc/hamiltonians/softabs_metric_test.cpp
@@ -1,27 +1,24 @@
 #include <stan/io/empty_var_context.hpp>
 #include <stan/mcmc/hmc/hamiltonians/softabs_metric.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/callbacks/stream_logger.hpp>
 #include <test/unit/mcmc/hmc/mock_hmc.hpp>
 #include <test/test-models/good/mcmc/hmc/hamiltonians/funnel.hpp>
 #include <test/unit/util.hpp>
 
-#include <boost/random/additive_combine.hpp>
-
 #include <gtest/gtest.h>
 
 #include <string>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcSoftAbs, sample_p) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
   q(1) = 1;
 
   stan::mcmc::mock_model model(q.size());
-  stan::mcmc::softabs_metric<stan::mcmc::mock_model, rng_t> metric(model);
+  stan::mcmc::softabs_metric<stan::mcmc::mock_model, stan::rng_t> metric(model);
   stan::mcmc::softabs_point z(q.size());
 
   int n_samples = 1000;
@@ -60,7 +57,6 @@ TEST(McmcSoftAbs, sample_p) {
 }
 
 TEST(McmcSoftAbs, gradients) {
-  rng_t base_rng(0);
 
   Eigen::VectorXd q = Eigen::VectorXd::Ones(11);
 
@@ -77,7 +73,7 @@ TEST(McmcSoftAbs, gradients) {
   funnel_model_namespace::funnel_model model(data_var_context, 0,
                                              &model_output);
 
-  stan::mcmc::softabs_metric<funnel_model_namespace::funnel_model, rng_t>
+  stan::mcmc::softabs_metric<funnel_model_namespace::funnel_model, stan::rng_t>
       metric(model);
 
   double epsilon = 1e-6;
@@ -152,7 +148,6 @@ TEST(McmcSoftAbs, gradients) {
 
 TEST(McmcSoftAbs, streams) {
   stan::test::capture_std_streams();
-  rng_t base_rng(0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -160,7 +155,8 @@ TEST(McmcSoftAbs, streams) {
   stan::mcmc::mock_model model(q.size());
 
   // for use in Google Test macros below
-  typedef stan::mcmc::softabs_metric<stan::mcmc::mock_model, rng_t> softabs;
+  typedef stan::mcmc::softabs_metric<stan::mcmc::mock_model, stan::rng_t>
+      softabs;
 
   EXPECT_NO_THROW(softabs metric(model));
 
diff --git a/src/test/unit/mcmc/hmc/hamiltonians/unit_e_metric_test.cpp b/src/test/unit/mcmc/hmc/hamiltonians/unit_e_metric_test.cpp
index 635386f55fc..177c8905aa0 100644
--- a/src/test/unit/mcmc/hmc/hamiltonians/unit_e_metric_test.cpp
+++ b/src/test/unit/mcmc/hmc/hamiltonians/unit_e_metric_test.cpp
@@ -1,5 +1,5 @@
 #include <string>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <test/unit/mcmc/hmc/mock_hmc.hpp>
 #include <stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp>
@@ -8,17 +8,15 @@
 #include <test/unit/util.hpp>
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcUnitEMetric, sample_p) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
   q(1) = 1;
 
   stan::mcmc::mock_model model(q.size());
-  stan::mcmc::unit_e_metric<stan::mcmc::mock_model, rng_t> metric(model);
+  stan::mcmc::unit_e_metric<stan::mcmc::mock_model, stan::rng_t> metric(model);
   stan::mcmc::unit_e_point z(q.size());
 
   int n_samples = 1000;
@@ -44,7 +42,6 @@ TEST(McmcUnitEMetric, sample_p) {
 }
 
 TEST(McmcUnitEMetric, gradients) {
-  rng_t base_rng(0);
 
   Eigen::VectorXd q = Eigen::VectorXd::Ones(11);
 
@@ -61,8 +58,8 @@ TEST(McmcUnitEMetric, gradients) {
   funnel_model_namespace::funnel_model model(data_var_context, 0,
                                              &model_output);
 
-  stan::mcmc::unit_e_metric<funnel_model_namespace::funnel_model, rng_t> metric(
-      model);
+  stan::mcmc::unit_e_metric<funnel_model_namespace::funnel_model, stan::rng_t>
+      metric(model);
 
   double epsilon = 1e-6;
 
@@ -136,7 +133,6 @@ TEST(McmcUnitEMetric, gradients) {
 
 TEST(McmcUnitEMetric, streams) {
   stan::test::capture_std_streams();
-  rng_t base_rng(0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -144,7 +140,7 @@ TEST(McmcUnitEMetric, streams) {
   stan::mcmc::mock_model model(q.size());
 
   // for use in Google Test macros below
-  typedef stan::mcmc::unit_e_metric<stan::mcmc::mock_model, rng_t> unit_e;
+  typedef stan::mcmc::unit_e_metric<stan::mcmc::mock_model, stan::rng_t> unit_e;
 
   EXPECT_NO_THROW(unit_e metric(model));
 
diff --git a/src/test/unit/mcmc/hmc/integrators/expl_leapfrog2_test.cpp b/src/test/unit/mcmc/hmc/integrators/expl_leapfrog2_test.cpp
index 489eddf5d28..dee3dc18cf5 100644
--- a/src/test/unit/mcmc/hmc/integrators/expl_leapfrog2_test.cpp
+++ b/src/test/unit/mcmc/hmc/integrators/expl_leapfrog2_test.cpp
@@ -7,12 +7,9 @@
 #include <stan/io/empty_var_context.hpp>
 #include <stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp>
 #include <stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
-
-typedef boost::ecuyer1988 rng_t;
+#include <stan/services/util/create_rng.hpp>
 
 TEST(McmcHmcIntegratorsExplLeapfrog, energy_conservation) {
-  rng_t base_rng(0);
 
   stan::io::empty_var_context data_var_context;
 
@@ -22,12 +19,12 @@ TEST(McmcHmcIntegratorsExplLeapfrog, energy_conservation) {
 
   gauss_model_namespace::gauss_model model(data_var_context, 0, &model_output);
 
-  stan::mcmc::expl_leapfrog<
-      stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, rng_t> >
+  stan::mcmc::expl_leapfrog<stan::mcmc::unit_e_metric<
+      gauss_model_namespace::gauss_model, stan::rng_t> >
       integrator;
 
-  stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, rng_t> metric(
-      model);
+  stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, stan::rng_t>
+      metric(model);
 
   stan::mcmc::unit_e_point z(1);
   z.q(0) = 1;
@@ -61,7 +58,6 @@ TEST(McmcHmcIntegratorsExplLeapfrog, energy_conservation) {
 }
 
 TEST(McmcHmcIntegratorsExplLeapfrog, symplecticness) {
-  rng_t base_rng(0);
 
   stan::io::empty_var_context data_var_context;
 
@@ -71,12 +67,12 @@ TEST(McmcHmcIntegratorsExplLeapfrog, symplecticness) {
 
   gauss_model_namespace::gauss_model model(data_var_context, 0, &model_output);
 
-  stan::mcmc::expl_leapfrog<
-      stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, rng_t> >
+  stan::mcmc::expl_leapfrog<stan::mcmc::unit_e_metric<
+      gauss_model_namespace::gauss_model, stan::rng_t> >
       integrator;
 
-  stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, rng_t> metric(
-      model);
+  stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, stan::rng_t>
+      metric(model);
 
   // Create a circle of points
   const int n_points = 1000;
diff --git a/src/test/unit/mcmc/hmc/integrators/expl_leapfrog_test.cpp b/src/test/unit/mcmc/hmc/integrators/expl_leapfrog_test.cpp
index 6dd6eaf419e..24a8b388e5b 100644
--- a/src/test/unit/mcmc/hmc/integrators/expl_leapfrog_test.cpp
+++ b/src/test/unit/mcmc/hmc/integrators/expl_leapfrog_test.cpp
@@ -9,13 +9,12 @@
 
 #include <stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp>
 #include <stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
+#include <stan/services/util/create_rng.hpp>
 #include <test/unit/util.hpp>
 
 // namespace
 //************************************************************
 
-typedef boost::ecuyer1988 rng_t;
 
 class McmcHmcIntegratorsExplLeapfrogF : public testing::Test {
  public:
@@ -45,11 +44,11 @@ class McmcHmcIntegratorsExplLeapfrogF : public testing::Test {
 
   // integrator under test
   stan::mcmc::expl_leapfrog<
-      stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t> >
+      stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t> >
       unit_e_integrator;
 
   stan::mcmc::expl_leapfrog<
-      stan::mcmc::diag_e_metric<command_model_namespace::command_model, rng_t> >
+      stan::mcmc::diag_e_metric<command_model_namespace::command_model, stan::rng_t> >
       diag_e_integrator;
 
   // model
@@ -72,7 +71,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, begin_update_p) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -104,7 +103,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, update_q) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -136,7 +135,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, end_update_p) {
   EXPECT_NEAR(z.g(0), 1.67264975797776, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -168,7 +167,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_1) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -200,7 +199,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_2) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -232,7 +231,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_3) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -264,7 +263,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_4) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -296,7 +295,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_5) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -327,7 +326,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_6) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -359,7 +358,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_7) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -391,7 +390,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_8) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -424,7 +423,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_9) {
   EXPECT_NEAR(z.g(0), 1.27097196280777, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::diag_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::diag_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -447,7 +446,7 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, streams) {
   stan::test::capture_std_streams();
 
   typedef stan::mcmc::expl_leapfrog<
-      stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t> >
+      stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t> >
       integrator;
 
   EXPECT_NO_THROW(integrator i);
diff --git a/src/test/unit/mcmc/hmc/integrators/impl_leapfrog2_test.cpp b/src/test/unit/mcmc/hmc/integrators/impl_leapfrog2_test.cpp
index 47606be8eb4..088de2066e1 100644
--- a/src/test/unit/mcmc/hmc/integrators/impl_leapfrog2_test.cpp
+++ b/src/test/unit/mcmc/hmc/integrators/impl_leapfrog2_test.cpp
@@ -1,4 +1,5 @@
 #include <stan/io/empty_var_context.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/mcmc/hmc/hamiltonians/softabs_metric.hpp>
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp>
@@ -6,16 +7,11 @@
 
 #include <test/test-models/good/mcmc/hmc/integrators/gauss.hpp>
 
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
-
 #include <gtest/gtest.h>
 
 #include <sstream>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcHmcIntegratorsImplLeapfrog, unit_e_energy_conservation) {
-  rng_t base_rng(0);
 
   stan::io::empty_var_context data_var_context;
 
@@ -25,12 +21,12 @@ TEST(McmcHmcIntegratorsImplLeapfrog, unit_e_energy_conservation) {
 
   gauss_model_namespace::gauss_model model(data_var_context, 0, &model_output);
 
-  stan::mcmc::impl_leapfrog<
-      stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, rng_t> >
+  stan::mcmc::impl_leapfrog<stan::mcmc::unit_e_metric<
+      gauss_model_namespace::gauss_model, stan::rng_t> >
       integrator;
 
-  stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, rng_t> metric(
-      model);
+  stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, stan::rng_t>
+      metric(model);
 
   stan::mcmc::unit_e_point z(1);
   z.q(0) = 1;
@@ -64,7 +60,6 @@ TEST(McmcHmcIntegratorsImplLeapfrog, unit_e_energy_conservation) {
 }
 
 TEST(McmcHmcIntegratorsImplLeapfrog, unit_e_symplecticness) {
-  rng_t base_rng(0);
 
   stan::io::empty_var_context data_var_context;
 
@@ -74,12 +69,12 @@ TEST(McmcHmcIntegratorsImplLeapfrog, unit_e_symplecticness) {
 
   gauss_model_namespace::gauss_model model(data_var_context, 0, &model_output);
 
-  stan::mcmc::impl_leapfrog<
-      stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, rng_t> >
+  stan::mcmc::impl_leapfrog<stan::mcmc::unit_e_metric<
+      gauss_model_namespace::gauss_model, stan::rng_t> >
       integrator;
 
-  stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, rng_t> metric(
-      model);
+  stan::mcmc::unit_e_metric<gauss_model_namespace::gauss_model, stan::rng_t>
+      metric(model);
 
   // Create a circle of points
   const int n_points = 1000;
@@ -151,7 +146,6 @@ TEST(McmcHmcIntegratorsImplLeapfrog, unit_e_symplecticness) {
 }
 
 TEST(McmcHmcIntegratorsImplLeapfrog, softabs_energy_conservation) {
-  rng_t base_rng(0);
 
   stan::io::empty_var_context data_var_context;
 
@@ -161,12 +155,12 @@ TEST(McmcHmcIntegratorsImplLeapfrog, softabs_energy_conservation) {
 
   gauss_model_namespace::gauss_model model(data_var_context, 0, &model_output);
 
-  stan::mcmc::impl_leapfrog<
-      stan::mcmc::softabs_metric<gauss_model_namespace::gauss_model, rng_t> >
+  stan::mcmc::impl_leapfrog<stan::mcmc::softabs_metric<
+      gauss_model_namespace::gauss_model, stan::rng_t> >
       integrator;
 
-  stan::mcmc::softabs_metric<gauss_model_namespace::gauss_model, rng_t> metric(
-      model);
+  stan::mcmc::softabs_metric<gauss_model_namespace::gauss_model, stan::rng_t>
+      metric(model);
 
   stan::mcmc::softabs_point z(1);
   z.q(0) = 1;
@@ -200,7 +194,6 @@ TEST(McmcHmcIntegratorsImplLeapfrog, softabs_energy_conservation) {
 }
 
 TEST(McmcHmcIntegratorsImplLeapfrog, softabs_symplecticness) {
-  rng_t base_rng(0);
 
   stan::io::empty_var_context data_var_context;
 
@@ -210,12 +203,12 @@ TEST(McmcHmcIntegratorsImplLeapfrog, softabs_symplecticness) {
 
   gauss_model_namespace::gauss_model model(data_var_context, 0, &model_output);
 
-  stan::mcmc::impl_leapfrog<
-      stan::mcmc::softabs_metric<gauss_model_namespace::gauss_model, rng_t> >
+  stan::mcmc::impl_leapfrog<stan::mcmc::softabs_metric<
+      gauss_model_namespace::gauss_model, stan::rng_t> >
       integrator;
 
-  stan::mcmc::softabs_metric<gauss_model_namespace::gauss_model, rng_t> metric(
-      model);
+  stan::mcmc::softabs_metric<gauss_model_namespace::gauss_model, stan::rng_t>
+      metric(model);
 
   // Create a circle of points
   const int n_points = 1000;
diff --git a/src/test/unit/mcmc/hmc/integrators/impl_leapfrog_test.cpp b/src/test/unit/mcmc/hmc/integrators/impl_leapfrog_test.cpp
index 8aba1b21dff..ce5a6325287 100644
--- a/src/test/unit/mcmc/hmc/integrators/impl_leapfrog_test.cpp
+++ b/src/test/unit/mcmc/hmc/integrators/impl_leapfrog_test.cpp
@@ -1,5 +1,6 @@
 #include <stan/io/dump.hpp>
 #include <stan/callbacks/stream_logger.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/mcmc/hmc/hamiltonians/softabs_metric.hpp>
 #include <stan/mcmc/hmc/integrators/impl_leapfrog.hpp>
 #include <stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp>
@@ -8,7 +9,6 @@
 #include <test/test-models/good/mcmc/hmc/integrators/command.hpp>
 #include <test/unit/util.hpp>
 
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
 
 #include <gtest/gtest.h>
 
@@ -17,8 +17,6 @@
 // namespace
 //************************************************************
 
-typedef boost::ecuyer1988 rng_t;
-
 class McmcHmcIntegratorsImplLeapfrogF : public testing::Test {
  public:
   McmcHmcIntegratorsImplLeapfrogF()
@@ -47,15 +45,15 @@ class McmcHmcIntegratorsImplLeapfrogF : public testing::Test {
 
   // integrator under test
   stan::mcmc::impl_leapfrog<
-      stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t> >
+      stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t> >
       unit_e_integrator;
 
   stan::mcmc::impl_leapfrog<
-      stan::mcmc::diag_e_metric<command_model_namespace::command_model, rng_t> >
+      stan::mcmc::diag_e_metric<command_model_namespace::command_model, stan::rng_t> >
       diag_e_integrator;
 
   stan::mcmc::impl_leapfrog<stan::mcmc::softabs_metric<
-      command_model_namespace::command_model, rng_t> >
+      command_model_namespace::command_model, stan::rng_t> >
       softabs_integrator;
 
   // model
@@ -75,7 +73,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, begin_update_p) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -107,7 +105,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, softabs_begin_update_p) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::softabs_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::softabs_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -141,7 +139,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, update_q) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -173,7 +171,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, softabs_update_q) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::softabs_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::softabs_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -207,7 +205,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, end_update_p) {
   EXPECT_NEAR(z.g(0), 1.67264975797776, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -239,7 +237,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, softabs_end_update_p) {
   EXPECT_NEAR(z.g(0), 1.67264975797776, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::softabs_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::softabs_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -273,7 +271,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, evolve_1) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -305,7 +303,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, evolve_2) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -337,7 +335,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, evolve_3) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -369,7 +367,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, evolve_4) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -401,7 +399,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, evolve_5) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -432,7 +430,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, evolve_6) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -464,7 +462,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, evolve_7) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -496,7 +494,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, evolve_8) {
   EXPECT_NEAR(z.g(0), 1.99987371079118, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -529,7 +527,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, evolve_9) {
   EXPECT_NEAR(z.g(0), 1.27097196280777, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::diag_e_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::diag_e_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -561,7 +559,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, softabs_evolve) {
   EXPECT_NEAR(z.g(0), 1.27097196280777, 1e-15);
 
   // setup hamiltonian
-  stan::mcmc::softabs_metric<command_model_namespace::command_model, rng_t>
+  stan::mcmc::softabs_metric<command_model_namespace::command_model, stan::rng_t>
       hamiltonian(*model);
 
   // setup epsilon
@@ -586,7 +584,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, streams) {
   stan::test::capture_std_streams();
 
   typedef stan::mcmc::impl_leapfrog<
-      stan::mcmc::unit_e_metric<command_model_namespace::command_model, rng_t> >
+      stan::mcmc::unit_e_metric<command_model_namespace::command_model, stan::rng_t> >
       integrator;
 
   EXPECT_NO_THROW(integrator i);
@@ -600,7 +598,7 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, softabs_streams) {
   stan::test::capture_std_streams();
 
   typedef stan::mcmc::impl_leapfrog<stan::mcmc::softabs_metric<
-      command_model_namespace::command_model, rng_t> >
+      command_model_namespace::command_model, stan::rng_t> >
       integrator;
 
   EXPECT_NO_THROW(integrator i);
diff --git a/src/test/unit/mcmc/hmc/nuts/base_nuts_test.cpp b/src/test/unit/mcmc/hmc/nuts/base_nuts_test.cpp
index 50527c03cda..87b2297b5cf 100644
--- a/src/test/unit/mcmc/hmc/nuts/base_nuts_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts/base_nuts_test.cpp
@@ -3,20 +3,18 @@
 #include <stan/mcmc/hmc/nuts/base_nuts.hpp>
 #include <stan/mcmc/hmc/integrators/expl_leapfrog.hpp>
 #include <vector>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 namespace stan {
 namespace mcmc {
 
-class mock_nuts
-    : public base_nuts<mock_model, mock_hamiltonian, mock_integrator, rng_t> {
+class mock_nuts : public base_nuts<mock_model, mock_hamiltonian,
+                                   mock_integrator, stan::rng_t> {
  public:
-  mock_nuts(const mock_model& m, rng_t& rng)
-      : base_nuts<mock_model, mock_hamiltonian, mock_integrator, rng_t>(m,
-                                                                        rng) {}
+  mock_nuts(const mock_model& m, stan::rng_t& rng)
+      : base_nuts<mock_model, mock_hamiltonian, mock_integrator, stan::rng_t>(
+          m, rng) {}
 
   bool compute_criterion(Eigen::VectorXd& p_sharp_minus,
                          Eigen::VectorXd& p_sharp_plus, Eigen::VectorXd& rho) {
@@ -24,13 +22,13 @@ class mock_nuts
   }
 };
 
-class rho_inspector_mock_nuts
-    : public base_nuts<mock_model, mock_hamiltonian, mock_integrator, rng_t> {
+class rho_inspector_mock_nuts : public base_nuts<mock_model, mock_hamiltonian,
+                                                 mock_integrator, stan::rng_t> {
  public:
   std::vector<double> rho_values;
-  rho_inspector_mock_nuts(const mock_model& m, rng_t& rng)
-      : base_nuts<mock_model, mock_hamiltonian, mock_integrator, rng_t>(m,
-                                                                        rng) {}
+  rho_inspector_mock_nuts(const mock_model& m, stan::rng_t& rng)
+      : base_nuts<mock_model, mock_hamiltonian, mock_integrator, stan::rng_t>(
+          m, rng) {}
 
   bool compute_criterion(Eigen::VectorXd& p_sharp_minus,
                          Eigen::VectorXd& p_sharp_plus, Eigen::VectorXd& rho) {
@@ -40,13 +38,14 @@ class rho_inspector_mock_nuts
 };
 
 class edge_inspector_mock_nuts
-    : public base_nuts<mock_model, mock_hamiltonian, mock_integrator, rng_t> {
+    : public base_nuts<mock_model, mock_hamiltonian, mock_integrator,
+                       stan::rng_t> {
  public:
   std::vector<double> p_sharp_minus_values;
   std::vector<double> p_sharp_plus_values;
-  edge_inspector_mock_nuts(const mock_model& m, rng_t& rng)
-      : base_nuts<mock_model, mock_hamiltonian, mock_integrator, rng_t>(m,
-                                                                        rng) {}
+  edge_inspector_mock_nuts(const mock_model& m, stan::rng_t& rng)
+      : base_nuts<mock_model, mock_hamiltonian, mock_integrator, stan::rng_t>(
+          m, rng) {}
 
   bool compute_criterion(Eigen::VectorXd& p_sharp_minus,
                          Eigen::VectorXd& p_sharp_plus, Eigen::VectorXd& rho) {
@@ -92,18 +91,18 @@ class divergent_hamiltonian : public base_hamiltonian<M, ps_point, BaseRNG> {
 };
 
 class divergent_nuts : public base_nuts<mock_model, divergent_hamiltonian,
-                                        expl_leapfrog, rng_t> {
+                                        expl_leapfrog, stan::rng_t> {
  public:
-  divergent_nuts(const mock_model& m, rng_t& rng)
-      : base_nuts<mock_model, divergent_hamiltonian, expl_leapfrog, rng_t>(
-          m, rng) {}
+  divergent_nuts(const mock_model& m, stan::rng_t& rng)
+      : base_nuts<mock_model, divergent_hamiltonian, expl_leapfrog,
+                  stan::rng_t>(m, rng) {}
 };
 
 }  // namespace mcmc
 }  // namespace stan
 
 TEST(McmcNutsBaseNuts, set_max_depth_test) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -123,7 +122,7 @@ TEST(McmcNutsBaseNuts, set_max_depth_test) {
 }
 
 TEST(McmcNutsBaseNuts, set_max_delta_test) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -138,7 +137,7 @@ TEST(McmcNutsBaseNuts, set_max_delta_test) {
 }
 
 TEST(McmcNutsBaseNuts, build_tree_test) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
   double init_momentum = 1.5;
@@ -199,7 +198,7 @@ TEST(McmcNutsBaseNuts, build_tree_test) {
 }
 
 TEST(McmcNutsBaseNuts, rho_aggregation_test) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
   double init_momentum = 1.5;
@@ -267,7 +266,7 @@ TEST(McmcNutsBaseNuts, rho_aggregation_test) {
 }
 
 TEST(McmcNutsBaseNuts, divergence_test) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
   double init_momentum = 1.5;
@@ -334,7 +333,7 @@ TEST(McmcNutsBaseNuts, divergence_test) {
 }
 
 TEST(McmcNutsBaseNuts, transition) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
   double init_momentum = 1.5;
@@ -374,7 +373,7 @@ TEST(McmcNutsBaseNuts, transition) {
 }
 
 TEST(McmcNutsBaseNuts, transition_egde_momenta) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
   double init_momentum = 1.5;
diff --git a/src/test/unit/mcmc/hmc/nuts/derived_nuts_test.cpp b/src/test/unit/mcmc/hmc/nuts/derived_nuts_test.cpp
index 1bf6e3dbb1c..a6d0669ed6b 100644
--- a/src/test/unit/mcmc/hmc/nuts/derived_nuts_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts/derived_nuts_test.cpp
@@ -5,17 +5,14 @@
 #include <stan/mcmc/hmc/hamiltonians/unit_e_point.hpp>
 #include <stan/mcmc/hmc/hamiltonians/diag_e_point.hpp>
 #include <stan/mcmc/hmc/hamiltonians/dense_e_point.hpp>
+#include <stan/services/util/create_rng.hpp>
 
 #include <test/unit/mcmc/hmc/mock_hmc.hpp>
 
-#include <boost/random/additive_combine.hpp>
-
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcNutsDerivedNuts, compute_criterion_unit_e) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
 
@@ -26,8 +23,8 @@ TEST(McmcNutsDerivedNuts, compute_criterion_unit_e) {
   Eigen::VectorXd rho(model_size);
 
   stan::mcmc::mock_model model(model_size);
-  stan::mcmc::unit_e_nuts<stan::mcmc::mock_model, rng_t> sampler(model,
-                                                                 base_rng);
+  stan::mcmc::unit_e_nuts<stan::mcmc::mock_model, stan::rng_t> sampler(
+      model, base_rng);
 
   start.q(0) = 1;
   start.p(0) = 1;
@@ -55,7 +52,7 @@ TEST(McmcNutsDerivedNuts, compute_criterion_unit_e) {
 }
 
 TEST(McmcNutsDerivedNuts, compute_criterion_diag_e) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
 
@@ -66,8 +63,8 @@ TEST(McmcNutsDerivedNuts, compute_criterion_diag_e) {
   Eigen::VectorXd rho(model_size);
 
   stan::mcmc::mock_model model(model_size);
-  stan::mcmc::diag_e_nuts<stan::mcmc::mock_model, rng_t> sampler(model,
-                                                                 base_rng);
+  stan::mcmc::diag_e_nuts<stan::mcmc::mock_model, stan::rng_t> sampler(
+      model, base_rng);
 
   start.q(0) = 1;
   start.p(0) = 1;
@@ -95,7 +92,7 @@ TEST(McmcNutsDerivedNuts, compute_criterion_diag_e) {
 }
 
 TEST(McmcNutsDerivedNuts, compute_criterion_dense_e) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
 
@@ -106,8 +103,8 @@ TEST(McmcNutsDerivedNuts, compute_criterion_dense_e) {
   Eigen::VectorXd rho(model_size);
 
   stan::mcmc::mock_model model(model_size);
-  stan::mcmc::dense_e_nuts<stan::mcmc::mock_model, rng_t> sampler(model,
-                                                                  base_rng);
+  stan::mcmc::dense_e_nuts<stan::mcmc::mock_model, stan::rng_t> sampler(
+      model, base_rng);
 
   start.q(0) = 1;
   start.p(0) = 1;
@@ -135,7 +132,7 @@ TEST(McmcNutsDerivedNuts, compute_criterion_dense_e) {
 }
 
 TEST(McmcNutsDerivedNuts, compute_criterion_softabs) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
 
@@ -146,8 +143,8 @@ TEST(McmcNutsDerivedNuts, compute_criterion_softabs) {
   Eigen::VectorXd rho(model_size);
 
   stan::mcmc::mock_model model(model_size);
-  stan::mcmc::softabs_nuts<stan::mcmc::mock_model, rng_t> sampler(model,
-                                                                  base_rng);
+  stan::mcmc::softabs_nuts<stan::mcmc::mock_model, stan::rng_t> sampler(
+      model, base_rng);
 
   start.q(0) = 1;
   start.p(0) = 1;
diff --git a/src/test/unit/mcmc/hmc/nuts/instantiation_test.cpp b/src/test/unit/mcmc/hmc/nuts/instantiation_test.cpp
index 241d2e38f69..bd0070944fd 100644
--- a/src/test/unit/mcmc/hmc/nuts/instantiation_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts/instantiation_test.cpp
@@ -6,16 +6,14 @@
 #include <stan/mcmc/hmc/nuts/adapt_unit_e_nuts.hpp>
 #include <stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp>
 #include <stan/mcmc/hmc/nuts/adapt_dense_e_nuts.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <fstream>
 
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcNuts, instantiaton_test) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   std::stringstream output;
   stan::callbacks::stream_writer writer(output);
@@ -25,21 +23,24 @@ TEST(McmcNuts, instantiaton_test) {
   stan::io::empty_var_context data_var_context;
   gauss3D_model_namespace::gauss3D_model model(data_var_context);
 
-  stan::mcmc::unit_e_nuts<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::unit_e_nuts<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       unit_e_sampler(model, base_rng);
 
-  stan::mcmc::diag_e_nuts<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::diag_e_nuts<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       diag_e_sampler(model, base_rng);
 
-  stan::mcmc::dense_e_nuts<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::dense_e_nuts<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       dense_e_sampler(model, base_rng);
 
-  stan::mcmc::adapt_unit_e_nuts<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::adapt_unit_e_nuts<gauss3D_model_namespace::gauss3D_model,
+                                stan::rng_t>
       adapt_unit_e_sampler(model, base_rng);
 
-  stan::mcmc::adapt_diag_e_nuts<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::adapt_diag_e_nuts<gauss3D_model_namespace::gauss3D_model,
+                                stan::rng_t>
       adapt_diag_e_sampler(model, base_rng);
 
-  stan::mcmc::adapt_dense_e_nuts<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::adapt_dense_e_nuts<gauss3D_model_namespace::gauss3D_model,
+                                 stan::rng_t>
       adapt_dense_e_sampler(model, base_rng);
 }
diff --git a/src/test/unit/mcmc/hmc/nuts/softabs_nuts_test.cpp b/src/test/unit/mcmc/hmc/nuts/softabs_nuts_test.cpp
index 7d4506b788b..f537eea022c 100644
--- a/src/test/unit/mcmc/hmc/nuts/softabs_nuts_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts/softabs_nuts_test.cpp
@@ -1,16 +1,14 @@
 #include <test/test-models/good/mcmc/hmc/common/gauss3D.hpp>
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/mcmc/hmc/nuts/softabs_nuts.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <fstream>
 
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcSoftAbsNuts, build_tree_test) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::softabs_point z_init(3);
   z_init.q(0) = 1;
@@ -26,7 +24,7 @@ TEST(McmcSoftAbsNuts, build_tree_test) {
   stan::io::empty_var_context data_var_context;
   gauss3D_model_namespace::gauss3D_model model(data_var_context);
 
-  stan::mcmc::softabs_nuts<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::softabs_nuts<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -97,7 +95,7 @@ TEST(McmcSoftAbsNuts, build_tree_test) {
 }
 
 TEST(McmcSoftAbsNuts, tree_boundary_test) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::softabs_point z_init(3);
   z_init.q(0) = 1;
@@ -116,7 +114,7 @@ TEST(McmcSoftAbsNuts, tree_boundary_test) {
   model_t model(data_var_context);
 
   // Compute expected tree boundaries
-  typedef stan::mcmc::softabs_metric<model_t, rng_t> metric_t;
+  typedef stan::mcmc::softabs_metric<model_t, stan::rng_t> metric_t;
   metric_t metric(model);
 
   stan::mcmc::impl_leapfrog<metric_t> softabs_integrator;
@@ -169,7 +167,7 @@ TEST(McmcSoftAbsNuts, tree_boundary_test) {
   Eigen::VectorXd p_sharp_backward_4 = metric.dtau_dp(z_test);
 
   // Check expected tree boundaries to those dynamically generated by NUTS
-  stan::mcmc::softabs_nuts<model_t, rng_t> sampler(model, base_rng);
+  stan::mcmc::softabs_nuts<model_t, stan::rng_t> sampler(model, base_rng);
 
   sampler.set_nominal_stepsize(epsilon);
   sampler.set_stepsize_jitter(0);
@@ -311,7 +309,7 @@ TEST(McmcSoftAbsNuts, tree_boundary_test) {
 }
 
 TEST(McmcSoftAbsNuts, transition_test) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::softabs_point z_init(3);
   z_init.q(0) = 1;
@@ -327,7 +325,7 @@ TEST(McmcSoftAbsNuts, transition_test) {
   stan::io::empty_var_context data_var_context;
   gauss3D_model_namespace::gauss3D_model model(data_var_context);
 
-  stan::mcmc::softabs_nuts<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::softabs_nuts<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
diff --git a/src/test/unit/mcmc/hmc/nuts/unit_e_nuts_test.cpp b/src/test/unit/mcmc/hmc/nuts/unit_e_nuts_test.cpp
index bd77c25b9a1..aeccff91bf7 100644
--- a/src/test/unit/mcmc/hmc/nuts/unit_e_nuts_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts/unit_e_nuts_test.cpp
@@ -1,16 +1,14 @@
 #include <test/test-models/good/mcmc/hmc/common/gauss3D.hpp>
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/mcmc/hmc/nuts/unit_e_nuts.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <fstream>
 
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcUnitENuts, build_tree_test) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::unit_e_point z_init(3);
   z_init.q(0) = 1;
@@ -26,7 +24,7 @@ TEST(McmcUnitENuts, build_tree_test) {
   stan::io::empty_var_context data_var_context;
   gauss3D_model_namespace::gauss3D_model model(data_var_context);
 
-  stan::mcmc::unit_e_nuts<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::unit_e_nuts<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -97,7 +95,7 @@ TEST(McmcUnitENuts, build_tree_test) {
 }
 
 TEST(McmcUnitENuts, tree_boundary_test) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::unit_e_point z_init(3);
   z_init.q(0) = 1;
@@ -116,7 +114,7 @@ TEST(McmcUnitENuts, tree_boundary_test) {
   model_t model(data_var_context);
 
   // Compute expected tree boundaries
-  typedef stan::mcmc::unit_e_metric<model_t, rng_t> metric_t;
+  typedef stan::mcmc::unit_e_metric<model_t, stan::rng_t> metric_t;
   metric_t metric(model);
 
   stan::mcmc::expl_leapfrog<metric_t> unit_e_integrator;
@@ -169,7 +167,7 @@ TEST(McmcUnitENuts, tree_boundary_test) {
   Eigen::VectorXd p_sharp_backward_4 = metric.dtau_dp(z_test);
 
   // Check expected tree boundaries to those dynamically generated by NUTS
-  stan::mcmc::unit_e_nuts<model_t, rng_t> sampler(model, base_rng);
+  stan::mcmc::unit_e_nuts<model_t, stan::rng_t> sampler(model, base_rng);
 
   sampler.set_nominal_stepsize(epsilon);
   sampler.set_stepsize_jitter(0);
@@ -311,7 +309,7 @@ TEST(McmcUnitENuts, tree_boundary_test) {
 }
 
 TEST(McmcUnitENuts, transition_test) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::unit_e_point z_init(3);
   z_init.q(0) = 1;
@@ -327,7 +325,7 @@ TEST(McmcUnitENuts, transition_test) {
   stan::io::empty_var_context data_var_context;
   gauss3D_model_namespace::gauss3D_model model(data_var_context);
 
-  stan::mcmc::unit_e_nuts<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::unit_e_nuts<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
diff --git a/src/test/unit/mcmc/hmc/nuts_classic/base_nuts_classic_test.cpp b/src/test/unit/mcmc/hmc/nuts_classic/base_nuts_classic_test.cpp
index 7fbd6304cf1..c3862c2c889 100644
--- a/src/test/unit/mcmc/hmc/nuts_classic/base_nuts_classic_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts_classic/base_nuts_classic_test.cpp
@@ -2,20 +2,19 @@
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/mcmc/hmc/nuts_classic/base_nuts_classic.hpp>
 #include <stan/mcmc/hmc/integrators/expl_leapfrog.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 namespace stan {
 namespace mcmc {
 
-class mock_nuts_classic : public base_nuts_classic<mock_model, mock_hamiltonian,
-                                                   mock_integrator, rng_t> {
+class mock_nuts_classic
+    : public base_nuts_classic<mock_model, mock_hamiltonian, mock_integrator,
+                               stan::rng_t> {
  public:
-  mock_nuts_classic(const mock_model& m, rng_t& rng)
-      : base_nuts_classic<mock_model, mock_hamiltonian, mock_integrator, rng_t>(
-          m, rng) {}
+  mock_nuts_classic(const mock_model& m, stan::rng_t& rng)
+      : base_nuts_classic<mock_model, mock_hamiltonian, mock_integrator,
+                          stan::rng_t>(m, rng) {}
 
  private:
   bool compute_criterion(ps_point& start, ps_point& finish,
@@ -61,11 +60,11 @@ class divergent_hamiltonian : public base_hamiltonian<M, ps_point, BaseRNG> {
 
 class divergent_nuts_classic
     : public base_nuts_classic<mock_model, divergent_hamiltonian, expl_leapfrog,
-                               rng_t> {
+                               stan::rng_t> {
  public:
-  divergent_nuts_classic(const mock_model& m, rng_t& rng)
+  divergent_nuts_classic(const mock_model& m, stan::rng_t& rng)
       : base_nuts_classic<mock_model, divergent_hamiltonian, expl_leapfrog,
-                          rng_t>(m, rng) {}
+                          stan::rng_t>(m, rng) {}
 
  private:
   bool compute_criterion(ps_point& start, ps_point& finish,
@@ -78,7 +77,7 @@ class divergent_nuts_classic
 }  // namespace stan
 
 TEST(McmcNutsBaseNutsClassic, set_max_depth) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -98,7 +97,7 @@ TEST(McmcNutsBaseNutsClassic, set_max_depth) {
 }
 
 TEST(McmcNutsBaseNuts, set_max_delta) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -113,7 +112,7 @@ TEST(McmcNutsBaseNuts, set_max_delta) {
 }
 
 TEST(McmcNutsBaseNutsClassic, build_tree) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
   double init_momentum = 1.5;
@@ -167,7 +166,7 @@ TEST(McmcNutsBaseNutsClassic, build_tree) {
 }
 
 TEST(McmcNutsBaseNutsClassic, slice_criterion) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
   double init_momentum = 1.5;
diff --git a/src/test/unit/mcmc/hmc/nuts_classic/derived_nuts_classic_test.cpp b/src/test/unit/mcmc/hmc/nuts_classic/derived_nuts_classic_test.cpp
index e9f4637ab4e..593c1e56189 100644
--- a/src/test/unit/mcmc/hmc/nuts_classic/derived_nuts_classic_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts_classic/derived_nuts_classic_test.cpp
@@ -7,12 +7,10 @@
 #include <stan/mcmc/hmc/nuts_classic/diag_e_nuts_classic.hpp>
 #include <stan/mcmc/hmc/nuts_classic/dense_e_nuts_classic.hpp>
 #include <test/unit/mcmc/hmc/mock_hmc.hpp>
-#include <boost/random/additive_combine.hpp>
-
-typedef boost::ecuyer1988 rng_t;
+#include <stan/services/util/create_rng.hpp>
 
 TEST(McmcDerivedNutsClassic, compute_criterion_unit_e) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
 
@@ -21,7 +19,7 @@ TEST(McmcDerivedNutsClassic, compute_criterion_unit_e) {
   Eigen::VectorXd rho(model_size);
 
   stan::mcmc::mock_model model(model_size);
-  stan::mcmc::unit_e_nuts_classic<stan::mcmc::mock_model, rng_t> sampler(
+  stan::mcmc::unit_e_nuts_classic<stan::mcmc::mock_model, stan::rng_t> sampler(
       model, base_rng);
 
   start.q(0) = 1;
@@ -46,7 +44,7 @@ TEST(McmcDerivedNutsClassic, compute_criterion_unit_e) {
 }
 
 TEST(McmcDerivedNutsClassic, compute_criterion_diag_e) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
 
@@ -55,7 +53,7 @@ TEST(McmcDerivedNutsClassic, compute_criterion_diag_e) {
   Eigen::VectorXd rho(model_size);
 
   stan::mcmc::mock_model model(model_size);
-  stan::mcmc::diag_e_nuts_classic<stan::mcmc::mock_model, rng_t> sampler(
+  stan::mcmc::diag_e_nuts_classic<stan::mcmc::mock_model, stan::rng_t> sampler(
       model, base_rng);
 
   start.q(0) = 1;
@@ -80,7 +78,7 @@ TEST(McmcDerivedNutsClassic, compute_criterion_diag_e) {
 }
 
 TEST(McmcDerivedNutsClassic, compute_criterion_dense_e) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
 
@@ -89,7 +87,7 @@ TEST(McmcDerivedNutsClassic, compute_criterion_dense_e) {
   Eigen::VectorXd rho(model_size);
 
   stan::mcmc::mock_model model(model_size);
-  stan::mcmc::dense_e_nuts_classic<stan::mcmc::mock_model, rng_t> sampler(
+  stan::mcmc::dense_e_nuts_classic<stan::mcmc::mock_model, stan::rng_t> sampler(
       model, base_rng);
 
   start.q(0) = 1;
diff --git a/src/test/unit/mcmc/hmc/nuts_classic/instantiation_test.cpp b/src/test/unit/mcmc/hmc/nuts_classic/instantiation_test.cpp
index 70f8f38977d..654f75079c1 100644
--- a/src/test/unit/mcmc/hmc/nuts_classic/instantiation_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts_classic/instantiation_test.cpp
@@ -6,16 +6,14 @@
 #include <stan/mcmc/hmc/nuts_classic/adapt_unit_e_nuts_classic.hpp>
 #include <stan/mcmc/hmc/nuts_classic/adapt_diag_e_nuts_classic.hpp>
 #include <stan/mcmc/hmc/nuts_classic/adapt_dense_e_nuts_classic.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <fstream>
 
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcNutsClassic, instantiaton_test) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   std::stringstream output;
   stan::callbacks::stream_writer writer(output);
@@ -25,25 +23,27 @@ TEST(McmcNutsClassic, instantiaton_test) {
   stan::io::empty_var_context data_var_context;
   gauss3D_model_namespace::gauss3D_model model(data_var_context);
 
-  stan::mcmc::unit_e_nuts_classic<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::unit_e_nuts_classic<gauss3D_model_namespace::gauss3D_model,
+                                  stan::rng_t>
       unit_e_sampler(model, base_rng);
 
-  stan::mcmc::diag_e_nuts_classic<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::diag_e_nuts_classic<gauss3D_model_namespace::gauss3D_model,
+                                  stan::rng_t>
       diag_e_sampler(model, base_rng);
 
   stan::mcmc::dense_e_nuts_classic<gauss3D_model_namespace::gauss3D_model,
-                                   rng_t>
+                                   stan::rng_t>
       dense_e_sampler(model, base_rng);
 
   stan::mcmc::adapt_unit_e_nuts_classic<gauss3D_model_namespace::gauss3D_model,
-                                        rng_t>
+                                        stan::rng_t>
       adapt_unit_e_sampler(model, base_rng);
 
   stan::mcmc::adapt_diag_e_nuts_classic<gauss3D_model_namespace::gauss3D_model,
-                                        rng_t>
+                                        stan::rng_t>
       adapt_diag_e_sampler(model, base_rng);
 
   stan::mcmc::adapt_dense_e_nuts_classic<gauss3D_model_namespace::gauss3D_model,
-                                         rng_t>
+                                         stan::rng_t>
       adapt_dense_e_sampler(model, base_rng);
 }
diff --git a/src/test/unit/mcmc/hmc/static/base_static_hmc_test.cpp b/src/test/unit/mcmc/hmc/static/base_static_hmc_test.cpp
index 632426ff270..4d571f6a927 100644
--- a/src/test/unit/mcmc/hmc/static/base_static_hmc_test.cpp
+++ b/src/test/unit/mcmc/hmc/static/base_static_hmc_test.cpp
@@ -2,28 +2,26 @@
 #include <stan/mcmc/hmc/static/base_static_hmc.hpp>
 #include <stan/callbacks/writer.hpp>
 #include <stan/callbacks/stream_writer.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 namespace stan {
 namespace mcmc {
 
 // Mock Static HMC
 class mock_static_hmc : public base_static_hmc<mock_model, mock_hamiltonian,
-                                               mock_integrator, rng_t> {
+                                               mock_integrator, stan::rng_t> {
  public:
-  mock_static_hmc(const mock_model& m, rng_t& rng)
-      : base_static_hmc<mock_model, mock_hamiltonian, mock_integrator, rng_t>(
-          m, rng) {}
+  mock_static_hmc(const mock_model& m, stan::rng_t& rng)
+      : base_static_hmc<mock_model, mock_hamiltonian, mock_integrator,
+                        stan::rng_t>(m, rng) {}
 };
 
 }  // namespace mcmc
 }  // namespace stan
 
 TEST(McmcStaticBaseStaticHMC, set_nominal_stepsize) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   std::vector<double> q(5, 1.0);
   std::vector<int> r(2, 2);
@@ -43,7 +41,7 @@ TEST(McmcStaticBaseStaticHMC, set_nominal_stepsize) {
 }
 
 TEST(McmcStaticBaseStaticHMC, set_T) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   std::vector<double> q(5, 1.0);
   std::vector<int> r(2, 2);
@@ -63,7 +61,7 @@ TEST(McmcStaticBaseStaticHMC, set_T) {
 }
 
 TEST(McmcStaticBaseStaticHMC, set_nominal_stepsize_and_T) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   std::vector<double> q(5, 1.0);
   std::vector<int> r(2, 2);
@@ -90,7 +88,7 @@ TEST(McmcStaticBaseStaticHMC, set_nominal_stepsize_and_T) {
 }
 
 TEST(McmcStaticBaseStaticHMC, set_nominal_stepsize_and_L) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   std::vector<double> q(5, 1.0);
   std::vector<int> r(2, 2);
diff --git a/src/test/unit/mcmc/hmc/static_uniform/base_static_uniform_test.cpp b/src/test/unit/mcmc/hmc/static_uniform/base_static_uniform_test.cpp
index 4b8892eb767..a9f81458d2c 100644
--- a/src/test/unit/mcmc/hmc/static_uniform/base_static_uniform_test.cpp
+++ b/src/test/unit/mcmc/hmc/static_uniform/base_static_uniform_test.cpp
@@ -5,27 +5,25 @@
 #include <stan/callbacks/stream_writer.hpp>
 
 #include <test/unit/mcmc/hmc/mock_hmc.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 namespace stan {
 namespace mcmc {
 // Mock Static HMC
 class mock_static_uniform
     : public base_static_uniform<mock_model, mock_hamiltonian, mock_integrator,
-                                 rng_t> {
+                                 stan::rng_t> {
  public:
-  mock_static_uniform(const mock_model& m, rng_t& rng)
+  mock_static_uniform(const mock_model& m, stan::rng_t& rng)
       : base_static_uniform<mock_model, mock_hamiltonian, mock_integrator,
-                            rng_t>(m, rng) {}
+                            stan::rng_t>(m, rng) {}
 };
 }  // namespace mcmc
 }  // namespace stan
 
 TEST(McmcBaseStaticUniform, set_nominal_stepsize) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   std::vector<double> q(5, 1.0);
   std::vector<int> r(2, 2);
@@ -45,7 +43,7 @@ TEST(McmcBaseStaticUniform, set_nominal_stepsize) {
 }
 
 TEST(McmcBaseStaticUniform, set_T) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   std::vector<double> q(5, 1.0);
   std::vector<int> r(2, 2);
@@ -65,7 +63,7 @@ TEST(McmcBaseStaticUniform, set_T) {
 }
 
 TEST(McmcBaseStaticUniform, set_nominal_stepsize_and_T) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   std::vector<double> q(5, 1.0);
   std::vector<int> r(2, 2);
@@ -92,7 +90,7 @@ TEST(McmcBaseStaticUniform, set_nominal_stepsize_and_T) {
 }
 
 TEST(McmcBaseStaticUniform, set_nominal_stepsize_and_L) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   std::vector<double> q(5, 1.0);
   std::vector<int> r(2, 2);
diff --git a/src/test/unit/mcmc/hmc/static_uniform/derived_static_uniform_test.cpp b/src/test/unit/mcmc/hmc/static_uniform/derived_static_uniform_test.cpp
index 89f16df1edc..f1c9057d168 100644
--- a/src/test/unit/mcmc/hmc/static_uniform/derived_static_uniform_test.cpp
+++ b/src/test/unit/mcmc/hmc/static_uniform/derived_static_uniform_test.cpp
@@ -8,17 +8,14 @@
 #include <stan/mcmc/hmc/static_uniform/adapt_softabs_static_uniform.hpp>
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/io/empty_var_context.hpp>
+#include <stan/services/util/create_rng.hpp>
 
 #include <test/test-models/good/mcmc/hmc/common/gauss.hpp>
 
-#include <boost/random/additive_combine.hpp>
-
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcStaticUniform, unit_e_transition) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::unit_e_point z_init(1);
   z_init.q(0) = 1;
@@ -30,7 +27,8 @@ TEST(McmcStaticUniform, unit_e_transition) {
   stan::io::empty_var_context data_var_context;
   gauss_model_namespace::gauss_model model(data_var_context);
 
-  stan::mcmc::unit_e_static_uniform<gauss_model_namespace::gauss_model, rng_t>
+  stan::mcmc::unit_e_static_uniform<gauss_model_namespace::gauss_model,
+                                    stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -54,7 +52,7 @@ TEST(McmcStaticUniform, unit_e_transition) {
 }
 
 TEST(McmcStaticUniform, diag_e_transition) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::diag_e_point z_init(1);
   z_init.q(0) = 1;
@@ -66,7 +64,8 @@ TEST(McmcStaticUniform, diag_e_transition) {
   stan::io::empty_var_context data_var_context;
   gauss_model_namespace::gauss_model model(data_var_context);
 
-  stan::mcmc::diag_e_static_uniform<gauss_model_namespace::gauss_model, rng_t>
+  stan::mcmc::diag_e_static_uniform<gauss_model_namespace::gauss_model,
+                                    stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -90,7 +89,7 @@ TEST(McmcStaticUniform, diag_e_transition) {
 }
 
 TEST(McmcStaticUniform, dense_e_transition) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::dense_e_point z_init(1);
   z_init.q(0) = 1;
@@ -102,7 +101,8 @@ TEST(McmcStaticUniform, dense_e_transition) {
   stan::io::empty_var_context data_var_context;
   gauss_model_namespace::gauss_model model(data_var_context);
 
-  stan::mcmc::dense_e_static_uniform<gauss_model_namespace::gauss_model, rng_t>
+  stan::mcmc::dense_e_static_uniform<gauss_model_namespace::gauss_model,
+                                     stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -126,7 +126,7 @@ TEST(McmcStaticUniform, dense_e_transition) {
 }
 
 TEST(McmcStaticUniform, softabs_transition) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::softabs_point z_init(1);
   z_init.q(0) = 1;
@@ -138,7 +138,8 @@ TEST(McmcStaticUniform, softabs_transition) {
   stan::io::empty_var_context data_var_context;
   gauss_model_namespace::gauss_model model(data_var_context);
 
-  stan::mcmc::softabs_static_uniform<gauss_model_namespace::gauss_model, rng_t>
+  stan::mcmc::softabs_static_uniform<gauss_model_namespace::gauss_model,
+                                     stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -162,7 +163,7 @@ TEST(McmcStaticUniform, softabs_transition) {
 }
 
 TEST(McmcStaticUniform, adapt_unit_e_transition) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::unit_e_point z_init(1);
   z_init.q(0) = 1;
@@ -175,7 +176,7 @@ TEST(McmcStaticUniform, adapt_unit_e_transition) {
   gauss_model_namespace::gauss_model model(data_var_context);
 
   stan::mcmc::adapt_unit_e_static_uniform<gauss_model_namespace::gauss_model,
-                                          rng_t>
+                                          stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -199,7 +200,7 @@ TEST(McmcStaticUniform, adapt_unit_e_transition) {
 }
 
 TEST(McmcStaticUniform, adapt_diag_e_transition) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::diag_e_point z_init(1);
   z_init.q(0) = 1;
@@ -212,7 +213,7 @@ TEST(McmcStaticUniform, adapt_diag_e_transition) {
   gauss_model_namespace::gauss_model model(data_var_context);
 
   stan::mcmc::adapt_diag_e_static_uniform<gauss_model_namespace::gauss_model,
-                                          rng_t>
+                                          stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -236,7 +237,7 @@ TEST(McmcStaticUniform, adapt_diag_e_transition) {
 }
 
 TEST(McmcStaticUniform, adapt_dense_e_transition) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::dense_e_point z_init(1);
   z_init.q(0) = 1;
@@ -249,7 +250,7 @@ TEST(McmcStaticUniform, adapt_dense_e_transition) {
   gauss_model_namespace::gauss_model model(data_var_context);
 
   stan::mcmc::adapt_dense_e_static_uniform<gauss_model_namespace::gauss_model,
-                                           rng_t>
+                                           stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -273,7 +274,7 @@ TEST(McmcStaticUniform, adapt_dense_e_transition) {
 }
 
 TEST(McmcStaticUniform, adapt_softabs_e_transition) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::softabs_point z_init(1);
   z_init.q(0) = 1;
@@ -286,7 +287,7 @@ TEST(McmcStaticUniform, adapt_softabs_e_transition) {
   gauss_model_namespace::gauss_model model(data_var_context);
 
   stan::mcmc::adapt_softabs_static_uniform<gauss_model_namespace::gauss_model,
-                                           rng_t>
+                                           stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
diff --git a/src/test/unit/mcmc/hmc/xhmc/base_xhmc_test.cpp b/src/test/unit/mcmc/hmc/xhmc/base_xhmc_test.cpp
index f27eb9f37ca..e2bb9200a57 100644
--- a/src/test/unit/mcmc/hmc/xhmc/base_xhmc_test.cpp
+++ b/src/test/unit/mcmc/hmc/xhmc/base_xhmc_test.cpp
@@ -2,20 +2,18 @@
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/mcmc/hmc/xhmc/base_xhmc.hpp>
 #include <stan/mcmc/hmc/integrators/expl_leapfrog.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 namespace stan {
 namespace mcmc {
 
-class mock_xhmc
-    : public base_xhmc<mock_model, mock_hamiltonian, mock_integrator, rng_t> {
+class mock_xhmc : public base_xhmc<mock_model, mock_hamiltonian,
+                                   mock_integrator, stan::rng_t> {
  public:
-  mock_xhmc(const mock_model& m, rng_t& rng)
-      : base_xhmc<mock_model, mock_hamiltonian, mock_integrator, rng_t>(m,
-                                                                        rng) {}
+  mock_xhmc(const mock_model& m, stan::rng_t& rng)
+      : base_xhmc<mock_model, mock_hamiltonian, mock_integrator, stan::rng_t>(
+          m, rng) {}
 };
 
 // Mock Hamiltonian
@@ -54,18 +52,18 @@ class divergent_hamiltonian : public base_hamiltonian<M, ps_point, BaseRNG> {
 };
 
 class divergent_xhmc : public base_xhmc<mock_model, divergent_hamiltonian,
-                                        expl_leapfrog, rng_t> {
+                                        expl_leapfrog, stan::rng_t> {
  public:
-  divergent_xhmc(const mock_model& m, rng_t& rng)
-      : base_xhmc<mock_model, divergent_hamiltonian, expl_leapfrog, rng_t>(
-          m, rng) {}
+  divergent_xhmc(const mock_model& m, stan::rng_t& rng)
+      : base_xhmc<mock_model, divergent_hamiltonian, expl_leapfrog,
+                  stan::rng_t>(m, rng) {}
 };
 
 }  // namespace mcmc
 }  // namespace stan
 
 TEST(McmcXHMCBaseXHMC, set_max_depth) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -83,7 +81,7 @@ TEST(McmcXHMCBaseXHMC, set_max_depth) {
 }
 
 TEST(McmcXHMCBaseXHMC, set_max_deltaH) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -98,7 +96,7 @@ TEST(McmcXHMCBaseXHMC, set_max_deltaH) {
 }
 
 TEST(McmcXHMCBaseXHMC, set_x_delta) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   Eigen::VectorXd q(2);
   q(0) = 5;
@@ -113,7 +111,7 @@ TEST(McmcXHMCBaseXHMC, set_x_delta) {
 }
 
 TEST(McmcXHMCBaseXHMC, build_tree) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
   double init_momentum = 1.5;
@@ -164,7 +162,7 @@ TEST(McmcXHMCBaseXHMC, build_tree) {
 }
 
 TEST(McmcXHMCBaseXHMC, divergence_test) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
   double init_momentum = 1.5;
@@ -223,7 +221,7 @@ TEST(McmcXHMCBaseXHMC, divergence_test) {
 }
 
 TEST(McmcXHMCBaseXHMC, transition) {
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   int model_size = 1;
   double init_momentum = 1.5;
diff --git a/src/test/unit/mcmc/hmc/xhmc/softabs_xhmc_test.cpp b/src/test/unit/mcmc/hmc/xhmc/softabs_xhmc_test.cpp
index 600785eaf76..9aaf38dfaee 100644
--- a/src/test/unit/mcmc/hmc/xhmc/softabs_xhmc_test.cpp
+++ b/src/test/unit/mcmc/hmc/xhmc/softabs_xhmc_test.cpp
@@ -1,16 +1,14 @@
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/mcmc/hmc/xhmc/softabs_xhmc.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <test/test-models/good/mcmc/hmc/common/gauss3D.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <fstream>
 
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcUnitEXHMC, build_tree) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::softabs_point z_init(3);
   z_init.q(0) = 1;
@@ -26,7 +24,7 @@ TEST(McmcUnitEXHMC, build_tree) {
   stan::io::empty_var_context data_var_context;
   gauss3D_model_namespace::gauss3D_model model(data_var_context);
 
-  stan::mcmc::softabs_xhmc<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::softabs_xhmc<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -81,7 +79,7 @@ TEST(McmcUnitEXHMC, build_tree) {
 }
 
 TEST(McmcUnitEXHMC, transition) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::softabs_point z_init(3);
   z_init.q(0) = 1;
@@ -97,7 +95,7 @@ TEST(McmcUnitEXHMC, transition) {
   stan::io::empty_var_context data_var_context;
   gauss3D_model_namespace::gauss3D_model model(data_var_context);
 
-  stan::mcmc::softabs_xhmc<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::softabs_xhmc<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -121,4 +119,4 @@ TEST(McmcUnitEXHMC, transition) {
   EXPECT_EQ("", warn.str());
   EXPECT_EQ("", error.str());
   EXPECT_EQ("", fatal.str());
-}
\ No newline at end of file
+}
diff --git a/src/test/unit/mcmc/hmc/xhmc/unit_e_xhmc_test.cpp b/src/test/unit/mcmc/hmc/xhmc/unit_e_xhmc_test.cpp
index a9c53459754..8580905f3ff 100644
--- a/src/test/unit/mcmc/hmc/xhmc/unit_e_xhmc_test.cpp
+++ b/src/test/unit/mcmc/hmc/xhmc/unit_e_xhmc_test.cpp
@@ -1,16 +1,14 @@
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/mcmc/hmc/xhmc/unit_e_xhmc.hpp>
-#include <boost/random/additive_combine.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <test/test-models/good/mcmc/hmc/common/gauss3D.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <fstream>
 
 #include <gtest/gtest.h>
 
-typedef boost::ecuyer1988 rng_t;
-
 TEST(McmcUnitEXHMC, build_tree) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::unit_e_point z_init(3);
   z_init.q(0) = 1;
@@ -26,7 +24,7 @@ TEST(McmcUnitEXHMC, build_tree) {
   stan::io::empty_var_context data_var_context;
   gauss3D_model_namespace::gauss3D_model model(data_var_context);
 
-  stan::mcmc::unit_e_xhmc<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::unit_e_xhmc<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
@@ -81,7 +79,7 @@ TEST(McmcUnitEXHMC, build_tree) {
 }
 
 TEST(McmcUnitEXHMC, transition) {
-  rng_t base_rng(4839294);
+  stan::rng_t base_rng = stan::services::util::create_rng(4839294, 0);
 
   stan::mcmc::unit_e_point z_init(3);
   z_init.q(0) = 1;
@@ -97,7 +95,7 @@ TEST(McmcUnitEXHMC, transition) {
   stan::io::empty_var_context data_var_context;
   gauss3D_model_namespace::gauss3D_model model(data_var_context);
 
-  stan::mcmc::unit_e_xhmc<gauss3D_model_namespace::gauss3D_model, rng_t>
+  stan::mcmc::unit_e_xhmc<gauss3D_model_namespace::gauss3D_model, stan::rng_t>
       sampler(model, base_rng);
 
   sampler.z() = z_init;
diff --git a/src/test/unit/model/array_functions_roundtrip_test.cpp b/src/test/unit/model/array_functions_roundtrip_test.cpp
index 0f37719c42e..9c5cffa4a87 100644
--- a/src/test/unit/model/array_functions_roundtrip_test.cpp
+++ b/src/test/unit/model/array_functions_roundtrip_test.cpp
@@ -1,4 +1,5 @@
 #include <stan/model/log_prob_grad.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <stan/io/empty_var_context.hpp>
 #include <stan/io/json/json_data.hpp>
 #include <test/test-models/good/model/parameters.hpp>
@@ -20,7 +21,7 @@ class ModelArrayFunctionsRoundtripTest : public testing::Test {
 
   stan::io::empty_var_context context;
   stan_model model;
-  boost::ecuyer1988 rng;
+  stan::rng_t rng;
   std::unique_ptr<stan::io::var_context> inits;
   std::stringstream out;
 
diff --git a/src/test/unit/services/util/initialize_test.cpp b/src/test/unit/services/util/initialize_test.cpp
index 768929da5c1..a0acf6db685 100644
--- a/src/test/unit/services/util/initialize_test.cpp
+++ b/src/test/unit/services/util/initialize_test.cpp
@@ -1,4 +1,5 @@
 #include <stan/services/util/initialize.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/util.hpp>
 #include <stan/callbacks/stream_writer.hpp>
@@ -24,7 +25,7 @@ class ServicesUtilInitialize : public testing::Test {
   stan::callbacks::stream_writer message;
   stan::test::unit::instrumented_logger logger;
   stan::test::unit::instrumented_writer init;
-  boost::ecuyer1988 rng;
+  stan::rng_t rng;
 };
 
 TEST_F(ServicesUtilInitialize, radius_zero__print_false) {
diff --git a/src/test/unit/services/util/run_adaptive_sampler_test.cpp b/src/test/unit/services/util/run_adaptive_sampler_test.cpp
index 476d4cbad74..9c7d892707b 100644
--- a/src/test/unit/services/util/run_adaptive_sampler_test.cpp
+++ b/src/test/unit/services/util/run_adaptive_sampler_test.cpp
@@ -1,4 +1,5 @@
 #include <stan/services/util/run_adaptive_sampler.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 #include <test/test-models/good/services/test_lp.hpp>
 #include <stan/callbacks/structured_writer.hpp>
@@ -27,12 +28,12 @@ class ServicesUtil : public testing::Test {
   stan::io::empty_var_context context;
   stan_model model;
   std::vector<double> cont_vector;
-  boost::ecuyer1988 rng;
+  stan::rng_t rng;
   stan::test::unit::instrumented_interrupt interrupt;
   stan::test::unit::instrumented_writer sample_writer, diagnostic_writer;
   stan::callbacks::structured_writer dummy_metric_writer;
   stan::test::unit::instrumented_logger logger;
-  stan::mcmc::adapt_unit_e_nuts<stan_model, boost::ecuyer1988> sampler;
+  stan::mcmc::adapt_unit_e_nuts<stan_model, stan::rng_t> sampler;
   int num_warmup, num_samples, num_thin, refresh;
   bool save_warmup;
 };
diff --git a/src/test/unit/services/util/run_sampler_test.cpp b/src/test/unit/services/util/run_sampler_test.cpp
index 1b51781569c..4b63d61d686 100644
--- a/src/test/unit/services/util/run_sampler_test.cpp
+++ b/src/test/unit/services/util/run_sampler_test.cpp
@@ -1,4 +1,5 @@
 #include <stan/services/util/run_sampler.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 #include <test/test-models/good/services/test_lp.hpp>
 #include <stan/io/empty_var_context.hpp>
@@ -73,7 +74,7 @@ class ServicesUtil : public testing::Test {
   stan::io::empty_var_context context;
   stan_model model;
   std::vector<double> cont_vector;
-  boost::ecuyer1988 rng;
+  stan::rng_t rng;
   stan::test::unit::instrumented_interrupt interrupt;
   stan::test::unit::instrumented_writer sample_writer, diagnostic_writer;
   stan::test::unit::instrumented_logger logger;
diff --git a/src/test/unit/variational/advi_messages_test.cpp b/src/test/unit/variational/advi_messages_test.cpp
index 486a8240c7e..5441891dfd5 100644
--- a/src/test/unit/variational/advi_messages_test.cpp
+++ b/src/test/unit/variational/advi_messages_test.cpp
@@ -6,9 +6,8 @@
 #include <vector>
 #include <string>
 #include <iostream>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
+#include <stan/services/util/create_rng.hpp>
 
-typedef boost::ecuyer1988 rng_t;
 
 class advi_test : public ::testing::Test {
  public:
@@ -39,10 +38,10 @@ class advi_test : public ::testing::Test {
     diagnostic_stream_.str("");
 
     advi_meanfield_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_meanfield, rng_t>(
+        stan_model, stan::variational::normal_meanfield, stan::rng_t>(
         *model_, cont_params_, base_rng_, 1, 100, 1, 1);
     advi_fullrank_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_fullrank, rng_t>(
+        stan_model, stan::variational::normal_fullrank, stan::rng_t>(
         *model_, cont_params_, base_rng_, 1, 100, 1, 1);
   }
 
@@ -56,8 +55,8 @@ class advi_test : public ::testing::Test {
   std::string err_msg2;
 
   stan::variational::advi<stan_model, stan::variational::normal_meanfield,
-                          rng_t> *advi_meanfield_;
-  stan::variational::advi<stan_model, stan::variational::normal_fullrank, rng_t>
+                          stan::rng_t> *advi_meanfield_;
+  stan::variational::advi<stan_model, stan::variational::normal_fullrank, stan::rng_t>
       *advi_fullrank_;
   std::stringstream model_stream_;
   std::stringstream log_stream_;
@@ -69,7 +68,7 @@ class advi_test : public ::testing::Test {
 
  private:
   stan_model *model_;
-  rng_t base_rng_;
+  stan::rng_t base_rng_;
   Eigen::VectorXd cont_params_;
 };
 
diff --git a/src/test/unit/variational/advi_multivar_no_constraint_test.cpp b/src/test/unit/variational/advi_multivar_no_constraint_test.cpp
index 4310938aec1..9c5fc960073 100644
--- a/src/test/unit/variational/advi_multivar_no_constraint_test.cpp
+++ b/src/test/unit/variational/advi_multivar_no_constraint_test.cpp
@@ -5,9 +5,8 @@
 #include <test/unit/util.hpp>
 #include <vector>
 #include <string>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
+#include <stan/services/util/create_rng.hpp>
 
-typedef boost::ecuyer1988 rng_t;
 typedef multivariate_no_constraint_model_namespace::
     multivariate_no_constraint_model Model;
 
@@ -21,7 +20,7 @@ TEST(advi_test, multivar_no_constraint_fullrank) {
   Model my_model(dummy_context);
 
   // RNG
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   // Other params
   int n_monte_carlo_grad = 10;
@@ -36,7 +35,7 @@ TEST(advi_test, multivar_no_constraint_fullrank) {
   cont_params(1) = 0.75;
 
   // ADVI
-  stan::variational::advi<Model, stan::variational::normal_fullrank, rng_t>
+  stan::variational::advi<Model, stan::variational::normal_fullrank, stan::rng_t>
       test_advi(my_model, cont_params, base_rng, n_monte_carlo_grad,
                 n_grad_samples, 100, 1);
 
@@ -139,7 +138,7 @@ TEST(advi_test, multivar_no_constraint_meanfield) {
   Model my_model(dummy_context);
 
   // RNG
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   // Other params
   int n_monte_carlo_grad = 10;
@@ -153,7 +152,7 @@ TEST(advi_test, multivar_no_constraint_meanfield) {
   cont_params(1) = 0.75;
 
   // ADVI
-  stan::variational::advi<Model, stan::variational::normal_meanfield, rng_t>
+  stan::variational::advi<Model, stan::variational::normal_meanfield, stan::rng_t>
       test_advi(my_model, cont_params, base_rng, n_monte_carlo_grad,
                 1e4,  // absurdly high!
                 100, 1);
diff --git a/src/test/unit/variational/advi_multivar_with_constraint_test.cpp b/src/test/unit/variational/advi_multivar_with_constraint_test.cpp
index 8730b3b0672..ae6bcf6d882 100644
--- a/src/test/unit/variational/advi_multivar_with_constraint_test.cpp
+++ b/src/test/unit/variational/advi_multivar_with_constraint_test.cpp
@@ -5,9 +5,8 @@
 #include <test/unit/util.hpp>
 #include <vector>
 #include <string>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
+#include <stan/services/util/create_rng.hpp>
 
-typedef boost::ecuyer1988 rng_t;
 typedef multivariate_with_constraint_model_namespace::
     multivariate_with_constraint_model Model;
 
@@ -21,7 +20,7 @@ TEST(advi_test, multivar_with_constraint_fullrank) {
   Model my_model(dummy_context);
 
   // RNG
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   // Other params
   int n_monte_carlo_grad = 10;
@@ -36,7 +35,8 @@ TEST(advi_test, multivar_with_constraint_fullrank) {
   cont_params(1) = 0.75;
 
   // ADVI
-  stan::variational::advi<Model, stan::variational::normal_fullrank, rng_t>
+  stan::variational::advi<Model, stan::variational::normal_fullrank,
+                          stan::rng_t>
       test_advi(my_model, cont_params, base_rng, n_monte_carlo_grad,
                 n_monte_carlo_elbo, 100, 1);
 
@@ -78,7 +78,7 @@ TEST(advi_test, multivar_with_constraint_meanfield) {
   Model my_model(dummy_context);
 
   // RNG
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   // Other params
   int n_monte_carlo_grad = 10;
@@ -93,7 +93,8 @@ TEST(advi_test, multivar_with_constraint_meanfield) {
   cont_params(1) = 0.75;
 
   // ADVI
-  stan::variational::advi<Model, stan::variational::normal_meanfield, rng_t>
+  stan::variational::advi<Model, stan::variational::normal_meanfield,
+                          stan::rng_t>
       test_advi(my_model, cont_params, base_rng, n_monte_carlo_grad,
                 n_monte_carlo_elbo, 100, 1);
 
diff --git a/src/test/unit/variational/advi_univar_no_constraint_test.cpp b/src/test/unit/variational/advi_univar_no_constraint_test.cpp
index f1394948ff6..ef8d5fd1dce 100644
--- a/src/test/unit/variational/advi_univar_no_constraint_test.cpp
+++ b/src/test/unit/variational/advi_univar_no_constraint_test.cpp
@@ -2,13 +2,12 @@
 #include <stan/variational/advi.hpp>
 #include <stan/callbacks/stream_logger.hpp>
 #include <test/unit/util.hpp>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 #include <cmath>
 #include <string>
 #include <vector>
 
-typedef boost::ecuyer1988 rng_t;
 typedef univariate_no_constraint_model_namespace::univariate_no_constraint_model
     Model;
 
@@ -22,7 +21,7 @@ TEST(advi_test, univar_no_constraint_fullrank) {
   Model my_model(dummy_context);
 
   // RNG
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   // Other params
   int n_monte_carlo_grad = 10;
@@ -35,7 +34,8 @@ TEST(advi_test, univar_no_constraint_fullrank) {
   cont_params(0) = -0.75;
 
   // ADVI
-  stan::variational::advi<Model, stan::variational::normal_fullrank, rng_t>
+  stan::variational::advi<Model, stan::variational::normal_fullrank,
+                          stan::rng_t>
       test_advi(my_model, cont_params, base_rng, n_monte_carlo_grad,
                 1e4,  // absurdly high!
                 100, 1);
@@ -146,7 +146,7 @@ TEST(advi_test, univar_no_constraint_meanfield) {
   Model my_model(dummy_context);
 
   // RNG
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   // Other params
   int n_monte_carlo_grad = 10;
@@ -159,7 +159,8 @@ TEST(advi_test, univar_no_constraint_meanfield) {
   cont_params(0) = -0.75;
 
   // ADVI
-  stan::variational::advi<Model, stan::variational::normal_meanfield, rng_t>
+  stan::variational::advi<Model, stan::variational::normal_meanfield,
+                          stan::rng_t>
       test_advi(my_model, cont_params, base_rng, n_monte_carlo_grad,
                 1e4,  // absurdly high!
                 100, 1);
diff --git a/src/test/unit/variational/advi_univar_with_constraint_test.cpp b/src/test/unit/variational/advi_univar_with_constraint_test.cpp
index c12ac4f4b27..d28844a01fa 100644
--- a/src/test/unit/variational/advi_univar_with_constraint_test.cpp
+++ b/src/test/unit/variational/advi_univar_with_constraint_test.cpp
@@ -2,13 +2,12 @@
 #include <stan/variational/advi.hpp>
 #include <stan/callbacks/stream_logger.hpp>
 #include <test/unit/util.hpp>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 #include <cmath>
 #include <string>
 #include <vector>
 
-typedef boost::ecuyer1988 rng_t;
 typedef univariate_with_constraint_model_namespace::
     univariate_with_constraint_model Model;
 
@@ -22,7 +21,7 @@ TEST(advi_test, univar_with_constraint_fullrank) {
   Model my_model(dummy_context);
 
   // RNG
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   // Other params
   int n_monte_carlo_grad = 10;
@@ -35,7 +34,8 @@ TEST(advi_test, univar_with_constraint_fullrank) {
   cont_params(0) = -0.75;
 
   // ADVI
-  stan::variational::advi<Model, stan::variational::normal_fullrank, rng_t>
+  stan::variational::advi<Model, stan::variational::normal_fullrank,
+                          stan::rng_t>
       test_advi(my_model, cont_params, base_rng, n_monte_carlo_grad,
                 5e5,  // absurdly high!
                 100, 1);
@@ -151,7 +151,7 @@ TEST(advi_test, univar_with_constraint_meanfield) {
   Model my_model(dummy_context);
 
   // RNG
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   // Other params
   int n_monte_carlo_grad = 10;
@@ -164,7 +164,8 @@ TEST(advi_test, univar_with_constraint_meanfield) {
   cont_params(0) = -0.75;
 
   // ADVI
-  stan::variational::advi<Model, stan::variational::normal_meanfield, rng_t>
+  stan::variational::advi<Model, stan::variational::normal_meanfield,
+                          stan::rng_t>
       test_advi(my_model, cont_params, base_rng, n_monte_carlo_grad,
                 5e5,  // absurdly high!
                 100, 1);
diff --git a/src/test/unit/variational/eta_adapt_big_test.cpp b/src/test/unit/variational/eta_adapt_big_test.cpp
index 080c666b328..b07808b7165 100644
--- a/src/test/unit/variational/eta_adapt_big_test.cpp
+++ b/src/test/unit/variational/eta_adapt_big_test.cpp
@@ -1,14 +1,12 @@
 #include <test/test-models/good/variational/eta_should_be_big.hpp>
 #include <stan/variational/advi.hpp>
 #include <stan/callbacks/stream_logger.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/util.hpp>
 #include <vector>
 #include <string>
 #include <iostream>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
-
-typedef boost::ecuyer1988 rng_t;
 
 class eta_adapt_big_test : public ::testing::Test {
  public:
@@ -28,11 +26,11 @@ class eta_adapt_big_test : public ::testing::Test {
     log_stream_.str("");
 
     advi_meanfield_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_meanfield, rng_t>(
+        stan_model, stan::variational::normal_meanfield, stan::rng_t>(
         *model_, cont_params_, base_rng_, 1, 100, 100, 1);
 
     advi_fullrank_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_fullrank, rng_t>(
+        stan_model, stan::variational::normal_fullrank, stan::rng_t>(
         *model_, cont_params_, base_rng_, 1, 100, 100, 1);
   }
 
@@ -43,15 +41,15 @@ class eta_adapt_big_test : public ::testing::Test {
   }
 
   stan::variational::advi<stan_model, stan::variational::normal_meanfield,
-                          rng_t> *advi_meanfield_;
-  stan::variational::advi<stan_model, stan::variational::normal_fullrank, rng_t>
+                          stan::rng_t> *advi_meanfield_;
+  stan::variational::advi<stan_model, stan::variational::normal_fullrank, stan::rng_t>
       *advi_fullrank_;
   std::stringstream model_stream_;
   std::stringstream log_stream_;
   stan::callbacks::stream_logger logger;
 
   stan_model *model_;
-  rng_t base_rng_;
+  stan::rng_t base_rng_;
   Eigen::VectorXd cont_params_;
 };
 
diff --git a/src/test/unit/variational/eta_adapt_fail_test.cpp b/src/test/unit/variational/eta_adapt_fail_test.cpp
index 1b98f930f6b..36d770aa028 100644
--- a/src/test/unit/variational/eta_adapt_fail_test.cpp
+++ b/src/test/unit/variational/eta_adapt_fail_test.cpp
@@ -1,14 +1,12 @@
 #include <test/test-models/good/variational/eta_should_fail.hpp>
 #include <stan/variational/advi.hpp>
 #include <stan/callbacks/stream_logger.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/util.hpp>
 #include <vector>
 #include <string>
 #include <iostream>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
-
-typedef boost::ecuyer1988 rng_t;
 
 class eta_should_fail_test : public ::testing::Test {
  public:
@@ -28,11 +26,11 @@ class eta_should_fail_test : public ::testing::Test {
     log_stream_.str("");
 
     advi_meanfield_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_meanfield, rng_t>(
+        stan_model, stan::variational::normal_meanfield, stan::rng_t>(
         *model_, cont_params_, base_rng_, 1, 1, 100, 1);
 
     advi_fullrank_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_fullrank, rng_t>(
+        stan_model, stan::variational::normal_fullrank, stan::rng_t>(
         *model_, cont_params_, base_rng_, 1, 1, 100, 1);
   }
 
@@ -43,15 +41,15 @@ class eta_should_fail_test : public ::testing::Test {
   }
 
   stan::variational::advi<stan_model, stan::variational::normal_meanfield,
-                          rng_t> *advi_meanfield_;
-  stan::variational::advi<stan_model, stan::variational::normal_fullrank, rng_t>
+                          stan::rng_t> *advi_meanfield_;
+  stan::variational::advi<stan_model, stan::variational::normal_fullrank, stan::rng_t>
       *advi_fullrank_;
   std::stringstream model_stream_;
   std::stringstream log_stream_;
   stan::callbacks::stream_logger logger;
 
   stan_model *model_;
-  rng_t base_rng_;
+  stan::rng_t base_rng_;
   Eigen::VectorXd cont_params_;
 };
 
diff --git a/src/test/unit/variational/eta_adapt_mock_models_test.cpp b/src/test/unit/variational/eta_adapt_mock_models_test.cpp
index 26d9ecb4b57..e4f7b3a627d 100644
--- a/src/test/unit/variational/eta_adapt_mock_models_test.cpp
+++ b/src/test/unit/variational/eta_adapt_mock_models_test.cpp
@@ -7,9 +7,6 @@
 #include <gtest/gtest.h>
 #include <test/unit/util.hpp>
 #include <stan/variational/advi.hpp>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
-
-typedef boost::ecuyer1988 rng_t;
 
 // Mock Model
 class mock_model : public stan::model::prob_grad {
diff --git a/src/test/unit/variational/eta_adapt_small_test.cpp b/src/test/unit/variational/eta_adapt_small_test.cpp
index 0b2fd28f79e..4328e6ff6a0 100644
--- a/src/test/unit/variational/eta_adapt_small_test.cpp
+++ b/src/test/unit/variational/eta_adapt_small_test.cpp
@@ -6,11 +6,9 @@
 #include <vector>
 #include <string>
 #include <iostream>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
+#include <stan/services/util/create_rng.hpp>
 #include <boost/version.hpp>
 
-typedef boost::ecuyer1988 rng_t;
-
 class eta_adapt_small_test : public ::testing::Test {
  public:
   eta_adapt_small_test()
@@ -29,11 +27,11 @@ class eta_adapt_small_test : public ::testing::Test {
     log_stream_.str("");
 
     advi_meanfield_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_meanfield, rng_t>(
+        stan_model, stan::variational::normal_meanfield, stan::rng_t>(
         *model_, cont_params_, base_rng_, 1, 100, 100, 1);
 
     advi_fullrank_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_fullrank, rng_t>(
+        stan_model, stan::variational::normal_fullrank, stan::rng_t>(
         *model_, cont_params_, base_rng_, 1, 100, 100, 1);
   }
 
@@ -44,15 +42,15 @@ class eta_adapt_small_test : public ::testing::Test {
   }
 
   stan::variational::advi<stan_model, stan::variational::normal_meanfield,
-                          rng_t> *advi_meanfield_;
-  stan::variational::advi<stan_model, stan::variational::normal_fullrank, rng_t>
+                          stan::rng_t> *advi_meanfield_;
+  stan::variational::advi<stan_model, stan::variational::normal_fullrank, stan::rng_t>
       *advi_fullrank_;
   std::stringstream model_stream_;
   std::stringstream log_stream_;
   stan::callbacks::stream_logger logger;
 
   stan_model *model_;
-  rng_t base_rng_;
+  stan::rng_t base_rng_;
   Eigen::VectorXd cont_params_;
 };
 
diff --git a/src/test/unit/variational/gradient_warn_test.cpp b/src/test/unit/variational/gradient_warn_test.cpp
index 2dab70c96dc..c8ac490b88a 100644
--- a/src/test/unit/variational/gradient_warn_test.cpp
+++ b/src/test/unit/variational/gradient_warn_test.cpp
@@ -2,14 +2,14 @@
 #include <stan/variational/advi.hpp>
 #include <stan/callbacks/stream_writer.hpp>
 #include <stan/callbacks/stream_logger.hpp>
+#include <stan/services/util/create_rng.hpp>
+
 #include <gtest/gtest.h>
 #include <test/unit/util.hpp>
 #include <vector>
 #include <string>
 #include <iostream>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
 
-typedef boost::ecuyer1988 rng_t;
 
 class advi_test : public ::testing::Test {
  public:
@@ -35,10 +35,10 @@ class advi_test : public ::testing::Test {
     diagnostic_stream_.str("");
 
     advi_meanfield_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_meanfield, rng_t>(
+        stan_model, stan::variational::normal_meanfield, stan::rng_t>(
         *model_, cont_params_, base_rng_, 1, 100, 100, 1);
     advi_fullrank_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_fullrank, rng_t>(
+        stan_model, stan::variational::normal_fullrank, stan::rng_t>(
         *model_, cont_params_, base_rng_, 1, 100, 100, 1);
   }
 
@@ -49,8 +49,8 @@ class advi_test : public ::testing::Test {
   }
 
   stan::variational::advi<stan_model, stan::variational::normal_meanfield,
-                          rng_t> *advi_meanfield_;
-  stan::variational::advi<stan_model, stan::variational::normal_fullrank, rng_t>
+                          stan::rng_t> *advi_meanfield_;
+  stan::variational::advi<stan_model, stan::variational::normal_fullrank, stan::rng_t>
       *advi_fullrank_;
   std::stringstream model_stream_;
   std::stringstream message_stream_;
@@ -63,7 +63,7 @@ class advi_test : public ::testing::Test {
 
  private:
   stan_model *model_;
-  rng_t base_rng_;
+  stan::rng_t base_rng_;
   Eigen::VectorXd cont_params_;
 };
 
diff --git a/src/test/unit/variational/hier_logistic_cp_test.cpp b/src/test/unit/variational/hier_logistic_cp_test.cpp
index 0905aeedc80..1c1fea788bf 100644
--- a/src/test/unit/variational/hier_logistic_cp_test.cpp
+++ b/src/test/unit/variational/hier_logistic_cp_test.cpp
@@ -2,14 +2,13 @@
 #include <stan/variational/advi.hpp>
 #include <stan/callbacks/writer.hpp>
 #include <stan/callbacks/stream_logger.hpp>
+#include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/util.hpp>
 #include <vector>
 #include <string>
 #include <iostream>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
 
-typedef boost::ecuyer1988 rng_t;
 typedef hier_logistic_cp_model_namespace::hier_logistic_cp_model Model_cp;
 
 TEST(advi_test, hier_logistic_cp_constraint_meanfield) {
@@ -27,13 +26,14 @@ TEST(advi_test, hier_logistic_cp_constraint_meanfield) {
   Model_cp my_model(data_var_context);
 
   // RNG
-  rng_t base_rng(0);
+  stan::rng_t base_rng = stan::services::util::create_rng(0, 0);
 
   // Dummy input
   Eigen::VectorXd cont_params = Eigen::VectorXd::Zero(my_model.num_params_r());
 
   // ADVI
-  stan::variational::advi<Model_cp, stan::variational::normal_meanfield, rng_t>
+  stan::variational::advi<Model_cp, stan::variational::normal_meanfield,
+                          stan::rng_t>
       test_advi(my_model, cont_params, base_rng, 10, 100, 100, 1);
 
   stan::callbacks::writer writer;
diff --git a/src/test/unit/variational/hier_logistic_test.cpp b/src/test/unit/variational/hier_logistic_test.cpp
index df0092aa120..c7e8ada23d0 100644
--- a/src/test/unit/variational/hier_logistic_test.cpp
+++ b/src/test/unit/variational/hier_logistic_test.cpp
@@ -5,11 +5,9 @@
 #include <vector>
 #include <string>
 #include <iostream>
-#include <boost/random/additive_combine.hpp>  // L'Ecuyer RNG
+#include <stan/services/util/create_rng.hpp>
 #include <stan/callbacks/stream_logger.hpp>
 
-typedef boost::ecuyer1988 rng_t;
-
 class advi_test : public ::testing::Test {
  public:
   advi_test()
@@ -38,10 +36,10 @@ class advi_test : public ::testing::Test {
     diagnostic_stream_.str("");
 
     advi_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_meanfield, rng_t>(
+        stan_model, stan::variational::normal_meanfield, stan::rng_t>(
         *model_, cont_params_, base_rng_, 10, 100, 100, 1);
     advi_fullrank_ = new stan::variational::advi<
-        stan_model, stan::variational::normal_fullrank, rng_t>(
+        stan_model, stan::variational::normal_fullrank, stan::rng_t>(
         *model_, cont_params_, base_rng_, 10, 100, 100, 1);
   }
 
@@ -53,8 +51,8 @@ class advi_test : public ::testing::Test {
   }
 
   stan::variational::advi<stan_model, stan::variational::normal_meanfield,
-                          rng_t> *advi_;
-  stan::variational::advi<stan_model, stan::variational::normal_fullrank, rng_t>
+                          stan::rng_t> *advi_;
+  stan::variational::advi<stan_model, stan::variational::normal_fullrank, stan::rng_t>
       *advi_fullrank_;
   std::stringstream model_stream_;
   std::stringstream message_stream_;
@@ -67,7 +65,7 @@ class advi_test : public ::testing::Test {
  private:
   stan_model *model_;
   stan_model *model_null_stream_;
-  rng_t base_rng_;
+  stan::rng_t base_rng_;
   Eigen::VectorXd cont_params_;
 };