diff --git a/stan/math/prim/core/init_threadpool_tbb.hpp b/stan/math/prim/core/init_threadpool_tbb.hpp
index 26c8915f387..cf1eec8e5e3 100644
--- a/stan/math/prim/core/init_threadpool_tbb.hpp
+++ b/stan/math/prim/core/init_threadpool_tbb.hpp
@@ -3,7 +3,7 @@
 
 #include <stan/math/prim/err/invalid_argument.hpp>
 
-#include <boost/lexical_cast.hpp>
+#include <charconv>
 
 #ifndef TBB_INTERFACE_NEW
 #include <tbb/tbb_stddef.h>
@@ -21,6 +21,7 @@
 #endif
 
 #include <cstdlib>
+#include <cstring>
 #include <thread>
 
 namespace stan {
@@ -48,9 +49,10 @@ inline int get_num_threads() {
 #ifdef STAN_THREADS
   const char* env_stan_num_threads = std::getenv("STAN_NUM_THREADS");
   if (env_stan_num_threads != nullptr) {
-    try {
-      const int env_num_threads
-          = boost::lexical_cast<int>(env_stan_num_threads);
+    int env_num_threads = 0;
+    const char* end = env_stan_num_threads + std::strlen(env_stan_num_threads);
+    auto result = std::from_chars(env_stan_num_threads, end, env_num_threads);
+    if (result.ec == std::errc() && result.ptr == end) {
       if (env_num_threads > 0) {
         num_threads = env_num_threads;
       } else if (env_num_threads == -1) {
@@ -61,7 +63,7 @@ inline int get_num_threads() {
                          "The STAN_NUM_THREADS environment variable is '",
                          "' but it must be positive or -1");
       }
-    } catch (const boost::bad_lexical_cast&) {
+    } else {
       invalid_argument("get_num_threads(int)", "STAN_NUM_THREADS",
                        env_stan_num_threads,
                        "The STAN_NUM_THREADS environment variable is '",
diff --git a/stan/math/prim/fun/grad_2F1.hpp b/stan/math/prim/fun/grad_2F1.hpp
index 42d5d2a50dd..8731c6d30b7 100644
--- a/stan/math/prim/fun/grad_2F1.hpp
+++ b/stan/math/prim/fun/grad_2F1.hpp
@@ -13,7 +13,6 @@
 #include <stan/math/prim/fun/sign.hpp>
 #include <stan/math/prim/fun/hypergeometric_2F1.hpp>
 #include <cmath>
-#include <boost/optional.hpp>
 
 namespace stan {
 namespace math {
diff --git a/stan/math/prim/fun/hypergeometric_2F1.hpp b/stan/math/prim/fun/hypergeometric_2F1.hpp
index ae327e033f1..674e6e46ab2 100644
--- a/stan/math/prim/fun/hypergeometric_2F1.hpp
+++ b/stan/math/prim/fun/hypergeometric_2F1.hpp
@@ -17,7 +17,7 @@
 #include <stan/math/prim/fun/sqrt.hpp>
 #include <stan/math/prim/fun/square.hpp>
 #include <stan/math/prim/fun/hypergeometric_pFq.hpp>
-#include <boost/optional.hpp>
+#include <optional>
 
 namespace stan {
 namespace math {
@@ -43,7 +43,7 @@ namespace internal {
  * @return Gauss hypergeometric function
  */
 template <typename Ta1, typename Ta2, typename Tb, typename Tz,
-          typename RtnT = boost::optional<return_type_t<Ta1, Ta2, Tb, Tz>>,
+          typename RtnT = std::optional<return_type_t<Ta1, Ta2, Tb, Tz>>,
           require_all_arithmetic_t<Ta1, Ta2, Tb, Tz>* = nullptr>
 inline RtnT hyper_2F1_special_cases(const Ta1& a1, const Ta2& a2, const Tb& b,
                                     const Tz& z) {
@@ -149,7 +149,7 @@ inline RtnT hyper_2F1_special_cases(const Ta1& a1, const Ta2& a2, const Tb& b,
  */
 template <typename Ta1, typename Ta2, typename Tb, typename Tz,
           typename ScalarT = return_type_t<Ta1, Ta2, Tb, Tz>,
-          typename OptT = boost::optional<ScalarT>,
+          typename OptT = std::optional<ScalarT>,
           require_all_arithmetic_t<Ta1, Ta2, Tb, Tz>* = nullptr>
 inline return_type_t<Ta1, Ta2, Tb, Tz> hypergeometric_2F1(const Ta1& a1,
                                                           const Ta2& a2,
@@ -168,15 +168,15 @@ inline return_type_t<Ta1, Ta2, Tb, Tz> hypergeometric_2F1(const Ta1& a1,
   // Check whether value can be calculated by any special-case rules
   // before estimating infinite sum
   OptT special_case_a1a2 = internal::hyper_2F1_special_cases(a1, a2, b, z);
-  if (special_case_a1a2.is_initialized()) {
-    return special_case_a1a2.get();
+  if (special_case_a1a2.has_value()) {
+    return special_case_a1a2.value();
   }
 
   // Check whether any special case rules apply with 'a' arguments reversed
   // as 2F1(a1, a2, b, z) = 2F1(a2, a1, b, z)
   OptT special_case_a2a1 = internal::hyper_2F1_special_cases(a2, a1, b, z);
-  if (special_case_a2a1.is_initialized()) {
-    return special_case_a2a1.get();
+  if (special_case_a2a1.has_value()) {
+    return special_case_a2a1.value();
   }
 
   Eigen::Matrix<double, 2, 1> a_args(2);
diff --git a/stan/math/prim/prob/bernoulli_logit_glm_rng.hpp b/stan/math/prim/prob/bernoulli_logit_glm_rng.hpp
index 5a65aad3e57..85fb4240627 100644
--- a/stan/math/prim/prob/bernoulli_logit_glm_rng.hpp
+++ b/stan/math/prim/prob/bernoulli_logit_glm_rng.hpp
@@ -8,8 +8,7 @@
 #include <stan/math/prim/fun/as_array_or_scalar.hpp>
 #include <stan/math/prim/fun/inv_logit.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/bernoulli_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 #include <vector>
 
 namespace stan {
@@ -42,8 +41,6 @@ namespace math {
 template <typename T_x, typename T_alpha, typename T_beta, class RNG>
 inline typename VectorBuilder<true, int, T_alpha>::type bernoulli_logit_glm_rng(
     const T_x& x, const T_alpha& alpha, const T_beta& beta, RNG& rng) {
-  using boost::bernoulli_distribution;
-  using boost::variate_generator;
   using T_x_ref = ref_type_t<T_x>;
   using T_alpha_ref = ref_type_t<T_alpha>;
   using T_beta_ref = ref_type_t<T_beta>;
@@ -76,9 +73,8 @@ inline typename VectorBuilder<true, int, T_alpha>::type bernoulli_logit_glm_rng(
 
   for (size_t m = 0; m < M; ++m) {
     double theta_m = alpha_vec[m] + x_beta(m);
-    variate_generator<RNG&, bernoulli_distribution<>> bernoulli_rng(
-        rng, bernoulli_distribution<>(inv_logit(theta_m)));
-    output[m] = bernoulli_rng();
+    std::bernoulli_distribution bernoulli_rng(inv_logit(theta_m));
+    output[m] = bernoulli_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/bernoulli_logit_rng.hpp b/stan/math/prim/prob/bernoulli_logit_rng.hpp
index 662657149a7..8e99a647c24 100644
--- a/stan/math/prim/prob/bernoulli_logit_rng.hpp
+++ b/stan/math/prim/prob/bernoulli_logit_rng.hpp
@@ -6,8 +6,7 @@
 #include <stan/math/prim/fun/inv_logit.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size.hpp>
-#include <boost/random/bernoulli_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -29,8 +28,6 @@ namespace math {
 template <typename T_t, class RNG>
 inline typename VectorBuilder<true, int, T_t>::type bernoulli_logit_rng(
     const T_t& t, RNG& rng) {
-  using boost::bernoulli_distribution;
-  using boost::variate_generator;
   ref_type_t<T_t> t_ref = t;
   check_finite("bernoulli_logit_rng", "Logit transformed probability parameter",
                t_ref);
@@ -40,9 +37,8 @@ inline typename VectorBuilder<true, int, T_t>::type bernoulli_logit_rng(
   VectorBuilder<true, int, T_t> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, bernoulli_distribution<> > bernoulli_rng(
-        rng, bernoulli_distribution<>(inv_logit(t_vec[n])));
-    output[n] = bernoulli_rng();
+    std::bernoulli_distribution bernoulli_rng(inv_logit(t_vec[n]));
+    output[n] = bernoulli_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/bernoulli_rng.hpp b/stan/math/prim/prob/bernoulli_rng.hpp
index 49848f04ed2..dea48804380 100644
--- a/stan/math/prim/prob/bernoulli_rng.hpp
+++ b/stan/math/prim/prob/bernoulli_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size.hpp>
-#include <boost/random/bernoulli_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -28,8 +27,6 @@ namespace math {
 template <typename T_theta, class RNG>
 inline typename VectorBuilder<true, int, T_theta>::type bernoulli_rng(
     const T_theta& theta, RNG& rng) {
-  using boost::bernoulli_distribution;
-  using boost::variate_generator;
   static constexpr const char* function = "bernoulli_rng";
   ref_type_t<T_theta> theta_ref = theta;
   check_bounded(function, "Probability parameter", value_of(theta_ref), 0.0,
@@ -40,9 +37,8 @@ inline typename VectorBuilder<true, int, T_theta>::type bernoulli_rng(
   VectorBuilder<true, int, T_theta> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, bernoulli_distribution<> > bernoulli_rng(
-        rng, bernoulli_distribution<>(theta_vec[n]));
-    output[n] = bernoulli_rng();
+    std::bernoulli_distribution bernoulli_rng(theta_vec[n]);
+    output[n] = bernoulli_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/beta_proportion_rng.hpp b/stan/math/prim/prob/beta_proportion_rng.hpp
index 61bb114f00f..de71fe8c6a6 100644
--- a/stan/math/prim/prob/beta_proportion_rng.hpp
+++ b/stan/math/prim/prob/beta_proportion_rng.hpp
@@ -6,7 +6,6 @@
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
-#include <boost/random/variate_generator.hpp>
 
 namespace stan {
 namespace math {
diff --git a/stan/math/prim/prob/beta_rng.hpp b/stan/math/prim/prob/beta_rng.hpp
index 9fdc839bf4e..899e633cc22 100644
--- a/stan/math/prim/prob/beta_rng.hpp
+++ b/stan/math/prim/prob/beta_rng.hpp
@@ -6,10 +6,8 @@
 #include <stan/math/prim/fun/log_sum_exp.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/gamma_distribution.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -35,9 +33,6 @@ namespace math {
 template <typename T_shape1, typename T_shape2, class RNG>
 inline typename VectorBuilder<true, double, T_shape1, T_shape2>::type beta_rng(
     const T_shape1 &alpha, const T_shape2 &beta, RNG &rng) {
-  using boost::variate_generator;
-  using boost::random::gamma_distribution;
-  using boost::random::uniform_real_distribution;
   using T_alpha_ref = ref_type_t<T_shape1>;
   using T_beta_ref = ref_type_t<T_shape2>;
   static constexpr const char *function = "beta_rng";
@@ -53,29 +48,24 @@ inline typename VectorBuilder<true, double, T_shape1, T_shape2>::type beta_rng(
   size_t N = max_size(alpha, beta);
   VectorBuilder<true, double, T_shape1, T_shape2> output(N);
 
-  variate_generator<RNG &, uniform_real_distribution<>> uniform_rng(
-      rng, uniform_real_distribution<>(0.0, 1.0));
+  std::uniform_real_distribution<> uniform_rng(0.0, 1.0);
   for (size_t n = 0; n < N; ++n) {
     // If alpha and beta are large, trust the usual ratio of gammas
     // method for generating beta random variables. If any parameter
     // is small, work in log space and use Marsaglia and Tsang's trick
     if (alpha_vec[n] > 1.0 && beta_vec[n] > 1.0) {
-      variate_generator<RNG &, gamma_distribution<>> rng_gamma_alpha(
-          rng, gamma_distribution<>(alpha_vec[n], 1.0));
-      variate_generator<RNG &, gamma_distribution<>> rng_gamma_beta(
-          rng, gamma_distribution<>(beta_vec[n], 1.0));
-      double a = rng_gamma_alpha();
-      double b = rng_gamma_beta();
+      std::gamma_distribution<> rng_gamma_alpha(alpha_vec[n], 1.0);
+      std::gamma_distribution<> rng_gamma_beta(beta_vec[n], 1.0);
+      double a = rng_gamma_alpha(rng);
+      double b = rng_gamma_beta(rng);
       output[n] = a / (a + b);
     } else {
-      variate_generator<RNG &, gamma_distribution<>> rng_gamma_alpha(
-          rng, gamma_distribution<>(alpha_vec[n] + 1, 1.0));
-      variate_generator<RNG &, gamma_distribution<>> rng_gamma_beta(
-          rng, gamma_distribution<>(beta_vec[n] + 1, 1.0));
-      double log_a = std::log(uniform_rng()) / alpha_vec[n]
-                     + std::log(rng_gamma_alpha());
-      double log_b
-          = std::log(uniform_rng()) / beta_vec[n] + std::log(rng_gamma_beta());
+      std::gamma_distribution<> rng_gamma_alpha(alpha_vec[n] + 1, 1.0);
+      std::gamma_distribution<> rng_gamma_beta(beta_vec[n] + 1, 1.0);
+      double log_a = std::log(uniform_rng(rng)) / alpha_vec[n]
+                     + std::log(rng_gamma_alpha(rng));
+      double log_b = std::log(uniform_rng(rng)) / beta_vec[n]
+                     + std::log(rng_gamma_beta(rng));
       double log_sum = log_sum_exp(log_a, log_b);
       output[n] = std::exp(log_a - log_sum);
     }
diff --git a/stan/math/prim/prob/binomial_rng.hpp b/stan/math/prim/prob/binomial_rng.hpp
index 3dab00a3649..f29ad9521a0 100644
--- a/stan/math/prim/prob/binomial_rng.hpp
+++ b/stan/math/prim/prob/binomial_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/binomial_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -31,8 +30,6 @@ namespace math {
 template <typename T_N, typename T_theta, class RNG>
 inline typename VectorBuilder<true, int, T_N, T_theta>::type binomial_rng(
     const T_N& N, const T_theta& theta, RNG& rng) {
-  using boost::binomial_distribution;
-  using boost::variate_generator;
   using T_N_ref = ref_type_t<T_N>;
   using T_theta_ref = ref_type_t<T_theta>;
   static constexpr const char* function = "binomial_rng";
@@ -50,10 +47,9 @@ inline typename VectorBuilder<true, int, T_N, T_theta>::type binomial_rng(
   VectorBuilder<true, int, T_N, T_theta> output(M);
 
   for (size_t m = 0; m < M; ++m) {
-    variate_generator<RNG&, binomial_distribution<> > binomial_rng(
-        rng, binomial_distribution<>(N_vec[m], theta_vec[m]));
+    std::binomial_distribution<> binomial_rng(N_vec[m], theta_vec[m]);
 
-    output[m] = binomial_rng();
+    output[m] = binomial_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/categorical_logit_rng.hpp b/stan/math/prim/prob/categorical_logit_rng.hpp
index 155cc80684b..396ef3a8c48 100644
--- a/stan/math/prim/prob/categorical_logit_rng.hpp
+++ b/stan/math/prim/prob/categorical_logit_rng.hpp
@@ -6,8 +6,7 @@
 #include <stan/math/prim/fun/cumulative_sum.hpp>
 #include <stan/math/prim/fun/softmax.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
-#include <boost/random/uniform_01.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -27,16 +26,14 @@ namespace math {
  */
 template <class RNG>
 inline int categorical_logit_rng(const Eigen::VectorXd& beta, RNG& rng) {
-  using boost::uniform_01;
-  using boost::variate_generator;
   static constexpr const char* function = "categorical_logit_rng";
   check_finite(function, "Log odds parameter", beta);
 
-  variate_generator<RNG&, uniform_01<> > uniform01_rng(rng, uniform_01<>());
+  std::uniform_real_distribution<> uniform01_rng(0, 1);
   Eigen::VectorXd theta = softmax(beta);
   Eigen::VectorXd index = cumulative_sum(theta);
 
-  double c = uniform01_rng();
+  double c = uniform01_rng(rng);
   int b = 0;
   while (c > index(b)) {
     b++;
diff --git a/stan/math/prim/prob/categorical_rng.hpp b/stan/math/prim/prob/categorical_rng.hpp
index fac6d198aea..87eb36bf075 100644
--- a/stan/math/prim/prob/categorical_rng.hpp
+++ b/stan/math/prim/prob/categorical_rng.hpp
@@ -4,8 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/cumulative_sum.hpp>
-#include <boost/random/uniform_01.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -13,19 +12,17 @@ namespace math {
 template <class RNG>
 inline int categorical_rng(
     const Eigen::Matrix<double, Eigen::Dynamic, 1>& theta, RNG& rng) {
-  using boost::uniform_01;
-  using boost::variate_generator;
   static constexpr const char* function = "categorical_rng";
   check_simplex(function, "Probabilities parameter", theta);
 
-  variate_generator<RNG&, uniform_01<> > uniform01_rng(rng, uniform_01<>());
+  std::uniform_real_distribution<> uniform01_rng(0, 1);
 
   Eigen::VectorXd index(theta.rows());
   index.setZero();
 
   index = cumulative_sum(theta);
 
-  double c = uniform01_rng();
+  double c = uniform01_rng(rng);
   int b = 0;
   while (c > index(b, 0)) {
     b++;
diff --git a/stan/math/prim/prob/cauchy_rng.hpp b/stan/math/prim/prob/cauchy_rng.hpp
index 21ddc51d505..4ed7a12854c 100644
--- a/stan/math/prim/prob/cauchy_rng.hpp
+++ b/stan/math/prim/prob/cauchy_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/cauchy_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -32,8 +31,6 @@ namespace math {
 template <typename T_loc, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_scale>::type cauchy_rng(
     const T_loc& mu, const T_scale& sigma, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::cauchy_distribution;
   static constexpr const char* function = "cauchy_rng";
   using T_mu_ref = ref_type_t<T_loc>;
   using T_sigma_ref = ref_type_t<T_scale>;
@@ -50,9 +47,8 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type cauchy_rng(
   VectorBuilder<true, double, T_loc, T_scale> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, cauchy_distribution<> > cauchy_rng(
-        rng, cauchy_distribution<>(mu_vec[n], sigma_vec[n]));
-    output[n] = cauchy_rng();
+    std::cauchy_distribution<> cauchy_rng(mu_vec[n], sigma_vec[n]);
+    output[n] = cauchy_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/chi_square_rng.hpp b/stan/math/prim/prob/chi_square_rng.hpp
index 559b67bd88e..1884233401f 100644
--- a/stan/math/prim/prob/chi_square_rng.hpp
+++ b/stan/math/prim/prob/chi_square_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size.hpp>
-#include <boost/random/chi_squared_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -27,8 +26,6 @@ namespace math {
 template <typename T_deg, class RNG>
 inline typename VectorBuilder<true, double, T_deg>::type chi_square_rng(
     const T_deg& nu, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::chi_squared_distribution;
   using T_nu_ref = ref_type_t<T_deg>;
   static constexpr const char* function = "chi_square_rng";
   T_nu_ref nu_ref = nu;
@@ -39,9 +36,8 @@ inline typename VectorBuilder<true, double, T_deg>::type chi_square_rng(
   VectorBuilder<true, double, T_deg> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, chi_squared_distribution<> > chi_square_rng(
-        rng, chi_squared_distribution<>(nu_vec[n]));
-    output[n] = chi_square_rng();
+    std::chi_squared_distribution<> chi_square_rng(nu_vec[n]);
+    output[n] = chi_square_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/dirichlet_rng.hpp b/stan/math/prim/prob/dirichlet_rng.hpp
index bd4cd0c6232..81f6a7d5796 100644
--- a/stan/math/prim/prob/dirichlet_rng.hpp
+++ b/stan/math/prim/prob/dirichlet_rng.hpp
@@ -6,10 +6,8 @@
 #include <stan/math/prim/fun/log.hpp>
 #include <stan/math/prim/fun/log_sum_exp.hpp>
 #include <boost/math/special_functions/gamma.hpp>
-#include <boost/random/gamma_distribution.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -38,23 +36,18 @@ namespace math {
 template <class RNG>
 inline Eigen::VectorXd dirichlet_rng(
     const Eigen::Matrix<double, Eigen::Dynamic, 1>& alpha, RNG& rng) {
-  using boost::gamma_distribution;
-  using boost::variate_generator;
-  using boost::random::uniform_real_distribution;
   using Eigen::VectorXd;
   using std::exp;
   using std::log;
 
   // separate algorithm if any parameter is less than 1
   if (alpha.minCoeff() < 1) {
-    variate_generator<RNG&, uniform_real_distribution<> > uniform_rng(
-        rng, uniform_real_distribution<>(0.0, 1.0));
+    std::uniform_real_distribution<> uniform_rng(0.0, 1.0);
     VectorXd log_y(alpha.size());
     for (int i = 0; i < alpha.size(); ++i) {
-      variate_generator<RNG&, gamma_distribution<> > gamma_rng(
-          rng, gamma_distribution<>(alpha(i) + 1, 1));
-      double log_u = log(uniform_rng());
-      log_y(i) = log(gamma_rng()) + log_u / alpha(i);
+      std::gamma_distribution<> gamma_rng(alpha(i) + 1, 1);
+      double log_u = log(uniform_rng(rng));
+      log_y(i) = log(gamma_rng(rng)) + log_u / alpha(i);
     }
     double log_sum_y = log_sum_exp(log_y);
     VectorXd theta(alpha.size());
@@ -67,9 +60,8 @@ inline Eigen::VectorXd dirichlet_rng(
   // standard normalized gamma algorithm
   Eigen::VectorXd y(alpha.rows());
   for (int i = 0; i < alpha.rows(); i++) {
-    variate_generator<RNG&, gamma_distribution<> > gamma_rng(
-        rng, gamma_distribution<>(alpha(i, 0), 1e-7));
-    y(i) = gamma_rng();
+    std::gamma_distribution<> gamma_rng(alpha(i, 0), 1e-7);
+    y(i) = gamma_rng(rng);
   }
   return y / y.sum();
 }
diff --git a/stan/math/prim/prob/discrete_range_rng.hpp b/stan/math/prim/prob/discrete_range_rng.hpp
index 2dacf745339..e410afea65b 100644
--- a/stan/math/prim/prob/discrete_range_rng.hpp
+++ b/stan/math/prim/prob/discrete_range_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/uniform_int_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -35,8 +34,6 @@ template <typename T_lower, typename T_upper, class RNG>
 inline typename VectorBuilder<true, int, T_lower, T_upper>::type
 discrete_range_rng(const T_lower& lower, const T_upper& upper, RNG& rng) {
   static constexpr const char* function = "discrete_range_rng";
-  using boost::variate_generator;
-  using boost::random::uniform_int_distribution;
   check_consistent_sizes(function, "Lower bound parameter", lower,
                          "Upper bound parameter", upper);
   check_greater_or_equal(function, "Upper bound parameter", upper, lower);
@@ -47,10 +44,10 @@ discrete_range_rng(const T_lower& lower, const T_upper& upper, RNG& rng) {
   VectorBuilder<true, int, T_lower, T_upper> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, uniform_int_distribution<>> discrete_range_rng(
-        rng, uniform_int_distribution<>(lower_vec[n], upper_vec[n]));
+    std::uniform_int_distribution<> discrete_range_rng(lower_vec[n],
+                                                       upper_vec[n]);
 
-    output[n] = discrete_range_rng();
+    output[n] = discrete_range_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/double_exponential_rng.hpp b/stan/math/prim/prob/double_exponential_rng.hpp
index a5ee59a89ac..bee581c0fdb 100644
--- a/stan/math/prim/prob/double_exponential_rng.hpp
+++ b/stan/math/prim/prob/double_exponential_rng.hpp
@@ -6,9 +6,8 @@
 #include <stan/math/prim/fun/log1m.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -34,8 +33,6 @@ namespace math {
 template <typename T_loc, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_scale>::type
 double_exponential_rng(const T_loc& mu, const T_scale& sigma, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::uniform_real_distribution;
   using T_mu_ref = ref_type_t<T_loc>;
   using T_sigma_ref = ref_type_t<T_scale>;
   static constexpr const char* function = "double_exponential_rng";
@@ -51,10 +48,9 @@ double_exponential_rng(const T_loc& mu, const T_scale& sigma, RNG& rng) {
   size_t N = max_size(mu, sigma);
   VectorBuilder<true, double, T_loc, T_scale> output(N);
 
-  variate_generator<RNG&, uniform_real_distribution<> > z_rng(
-      rng, uniform_real_distribution<>(-1.0, 1.0));
+  std::uniform_real_distribution<> z_rng(-1.0, 1.0);
   for (size_t n = 0; n < N; ++n) {
-    double z = z_rng();
+    double z = z_rng(rng);
     output[n] = mu_vec[n]
                 - ((z > 0) ? 1.0 : -1.0) * sigma_vec[n] * std::log(std::abs(z));
   }
diff --git a/stan/math/prim/prob/exp_mod_normal_rng.hpp b/stan/math/prim/prob/exp_mod_normal_rng.hpp
index 84df8fd1733..6da82dfa54b 100644
--- a/stan/math/prim/prob/exp_mod_normal_rng.hpp
+++ b/stan/math/prim/prob/exp_mod_normal_rng.hpp
@@ -7,8 +7,6 @@
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/prob/normal_rng.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 
 namespace stan {
 namespace math {
diff --git a/stan/math/prim/prob/exponential_rng.hpp b/stan/math/prim/prob/exponential_rng.hpp
index 5d97ca02916..ffde3c8341e 100644
--- a/stan/math/prim/prob/exponential_rng.hpp
+++ b/stan/math/prim/prob/exponential_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size.hpp>
-#include <boost/random/exponential_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -27,8 +26,6 @@ namespace math {
 template <typename T_inv, class RNG>
 inline typename VectorBuilder<true, double, T_inv>::type exponential_rng(
     const T_inv& beta, RNG& rng) {
-  using boost::exponential_distribution;
-  using boost::variate_generator;
   static constexpr const char* function = "exponential_rng";
   using T_beta_ref = ref_type_t<T_inv>;
   T_beta_ref beta_ref = beta;
@@ -39,9 +36,8 @@ inline typename VectorBuilder<true, double, T_inv>::type exponential_rng(
   VectorBuilder<true, double, T_inv> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, exponential_distribution<> > exp_rng(
-        rng, exponential_distribution<>(beta_vec[n]));
-    output[n] = exp_rng();
+    std::exponential_distribution<> exp_rng(beta_vec[n]);
+    output[n] = exp_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/frechet_cdf.hpp b/stan/math/prim/prob/frechet_cdf.hpp
index e57a90f0bb0..b77e5df7a7a 100644
--- a/stan/math/prim/prob/frechet_cdf.hpp
+++ b/stan/math/prim/prob/frechet_cdf.hpp
@@ -14,8 +14,6 @@
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
-#include <boost/random/weibull_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
 #include <cmath>
 
diff --git a/stan/math/prim/prob/frechet_lccdf.hpp b/stan/math/prim/prob/frechet_lccdf.hpp
index 6226d64c92d..b0c43c07178 100644
--- a/stan/math/prim/prob/frechet_lccdf.hpp
+++ b/stan/math/prim/prob/frechet_lccdf.hpp
@@ -14,8 +14,6 @@
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
-#include <boost/random/weibull_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
 #include <cmath>
 
diff --git a/stan/math/prim/prob/frechet_lcdf.hpp b/stan/math/prim/prob/frechet_lcdf.hpp
index 4770deaa202..4d96c4c9bff 100644
--- a/stan/math/prim/prob/frechet_lcdf.hpp
+++ b/stan/math/prim/prob/frechet_lcdf.hpp
@@ -12,8 +12,6 @@
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
-#include <boost/random/weibull_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
 #include <cmath>
 
diff --git a/stan/math/prim/prob/frechet_lpdf.hpp b/stan/math/prim/prob/frechet_lpdf.hpp
index 3dba3f3ffd9..5569f21a51c 100644
--- a/stan/math/prim/prob/frechet_lpdf.hpp
+++ b/stan/math/prim/prob/frechet_lpdf.hpp
@@ -15,8 +15,6 @@
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
-#include <boost/random/weibull_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
 #include <cmath>
 
diff --git a/stan/math/prim/prob/frechet_rng.hpp b/stan/math/prim/prob/frechet_rng.hpp
index 0c15df9bc35..16ef0a61634 100644
--- a/stan/math/prim/prob/frechet_rng.hpp
+++ b/stan/math/prim/prob/frechet_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/weibull_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -31,8 +30,6 @@ namespace math {
 template <typename T_shape, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_shape, T_scale>::type frechet_rng(
     const T_shape& alpha, const T_scale& sigma, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::weibull_distribution;
   using T_alpha_ref = ref_type_t<T_shape>;
   using T_sigma_ref = ref_type_t<T_scale>;
   static constexpr const char* function = "frechet_rng";
@@ -49,9 +46,8 @@ inline typename VectorBuilder<true, double, T_shape, T_scale>::type frechet_rng(
   VectorBuilder<true, double, T_shape, T_scale> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, weibull_distribution<> > weibull_rng(
-        rng, weibull_distribution<>(alpha_vec[n], 1.0 / sigma_vec[n]));
-    output[n] = 1 / weibull_rng();
+    std::weibull_distribution<> weibull_rng(alpha_vec[n], 1.0 / sigma_vec[n]);
+    output[n] = 1 / weibull_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/gamma_cdf.hpp b/stan/math/prim/prob/gamma_cdf.hpp
index 4112a61ac39..2487a69e28a 100644
--- a/stan/math/prim/prob/gamma_cdf.hpp
+++ b/stan/math/prim/prob/gamma_cdf.hpp
@@ -15,8 +15,6 @@
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/tgamma.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
-#include <boost/random/gamma_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <stan/math/prim/functor/partials_propagator.hpp>
 #include <cmath>
 
diff --git a/stan/math/prim/prob/gamma_rng.hpp b/stan/math/prim/prob/gamma_rng.hpp
index c609e34214e..5f50038c59f 100644
--- a/stan/math/prim/prob/gamma_rng.hpp
+++ b/stan/math/prim/prob/gamma_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
-#include <boost/random/gamma_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -32,8 +31,6 @@ namespace math {
 template <typename T_shape, typename T_inv, class RNG>
 inline typename VectorBuilder<true, double, T_shape, T_inv>::type gamma_rng(
     const T_shape& alpha, const T_inv& beta, RNG& rng) {
-  using boost::gamma_distribution;
-  using boost::variate_generator;
   using T_alpha_ref = ref_type_t<T_shape>;
   using T_beta_ref = ref_type_t<T_inv>;
   static constexpr const char* function = "gamma_rng";
@@ -50,11 +47,10 @@ inline typename VectorBuilder<true, double, T_shape, T_inv>::type gamma_rng(
   VectorBuilder<true, double, T_shape, T_inv> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    // Convert rate (inverse scale) argument to scale for boost
-    variate_generator<RNG&, gamma_distribution<> > gamma_rng(
-        rng, gamma_distribution<>(alpha_vec[n],
-                                  1 / static_cast<double>(beta_vec[n])));
-    output[n] = gamma_rng();
+    // Convert rate (inverse scale) argument to scale for std
+    std::gamma_distribution<> gamma_rng(alpha_vec[n],
+                                        1 / static_cast<double>(beta_vec[n]));
+    output[n] = gamma_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/gaussian_dlm_obs_rng.hpp b/stan/math/prim/prob/gaussian_dlm_obs_rng.hpp
index b25d77fc07c..edfc6c6fc77 100644
--- a/stan/math/prim/prob/gaussian_dlm_obs_rng.hpp
+++ b/stan/math/prim/prob/gaussian_dlm_obs_rng.hpp
@@ -4,8 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 #include <vector>
 
 namespace stan {
@@ -29,17 +28,13 @@ template <class RNG>
 inline Eigen::VectorXd multi_normal_semidefinite_rng(
     const Eigen::VectorXd &mu, const Eigen::LDLT<Eigen::MatrixXd> &S_ldlt,
     RNG &rng) {
-  using boost::normal_distribution;
-  using boost::variate_generator;
-
-  variate_generator<RNG &, normal_distribution<>> std_normal_rng(
-      rng, normal_distribution<>(0, 1));
+  std::normal_distribution<> std_normal_rng(0, 1);
 
   Eigen::VectorXd stddev = S_ldlt.vectorD().array().sqrt().matrix();
   size_t M = S_ldlt.vectorD().size();
   Eigen::VectorXd z(M);
   for (int i = 0; i < M; i++) {
-    z(i) = stddev(i) * std_normal_rng();
+    z(i) = stddev(i) * std_normal_rng(rng);
   }
 
   Eigen::VectorXd Y
diff --git a/stan/math/prim/prob/gumbel_rng.hpp b/stan/math/prim/prob/gumbel_rng.hpp
index b6ece64a4ec..87eb21d7d16 100644
--- a/stan/math/prim/prob/gumbel_rng.hpp
+++ b/stan/math/prim/prob/gumbel_rng.hpp
@@ -5,9 +5,8 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
-#include <boost/random/uniform_01.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -33,8 +32,6 @@ namespace math {
 template <typename T_loc, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_scale>::type gumbel_rng(
     const T_loc& mu, const T_scale& beta, RNG& rng) {
-  using boost::uniform_01;
-  using boost::variate_generator;
   using T_mu_ref = ref_type_t<T_loc>;
   using T_beta_ref = ref_type_t<T_scale>;
   static constexpr const char* function = "gumbel_rng";
@@ -50,9 +47,10 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type gumbel_rng(
   size_t N = max_size(mu, beta);
   VectorBuilder<true, double, T_loc, T_scale> output(N);
 
-  variate_generator<RNG&, uniform_01<> > uniform01_rng(rng, uniform_01<>());
+  std::uniform_real_distribution<> uniform01_rng(0, 1);
   for (size_t n = 0; n < N; ++n) {
-    output[n] = mu_vec[n] - beta_vec[n] * std::log(-std::log(uniform01_rng()));
+    output[n]
+        = mu_vec[n] - beta_vec[n] * std::log(-std::log(uniform01_rng(rng)));
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/hmm_hidden_state_prob.hpp b/stan/math/prim/prob/hmm_hidden_state_prob.hpp
index d90731d23ad..083f5829e4a 100644
--- a/stan/math/prim/prob/hmm_hidden_state_prob.hpp
+++ b/stan/math/prim/prob/hmm_hidden_state_prob.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err/hmm_check.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
-#include <boost/random.hpp>
 #include <vector>
 
 namespace stan {
diff --git a/stan/math/prim/prob/hmm_latent_rng.hpp b/stan/math/prim/prob/hmm_latent_rng.hpp
index 24037e37df3..c7c0ad2177a 100644
--- a/stan/math/prim/prob/hmm_latent_rng.hpp
+++ b/stan/math/prim/prob/hmm_latent_rng.hpp
@@ -5,7 +5,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err/hmm_check.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
-#include <boost/random.hpp>
+#include <random>
 #include <vector>
 
 namespace stan {
@@ -69,7 +69,7 @@ inline std::vector<int> hmm_latent_rng(const T_omega& log_omegas,
   std::vector<double> probs(n_states);
   Eigen::Map<Eigen::VectorXd> probs_vec(probs.data(), n_states);
   probs_vec = alphas.col(n_transitions) / alphas.col(n_transitions).sum();
-  boost::random::discrete_distribution<> cat_hidden(probs);
+  std::discrete_distribution<> cat_hidden(probs.begin(), probs.end());
   hidden_states[n_transitions] = cat_hidden(rng) + stan::error_index::value;
 
   for (int n = n_transitions; n-- > 0;) {
@@ -84,7 +84,7 @@ inline std::vector<int> hmm_latent_rng(const T_omega& log_omegas,
 
     // discrete_distribution produces samples in [0, K), so
     // we need to add 1 to generate over [1, K).
-    boost::random::discrete_distribution<> cat_hidden(probs);
+    std::discrete_distribution<> cat_hidden(probs.begin(), probs.end());
     hidden_states[n] = cat_hidden(rng) + stan::error_index::value;
 
     // update backwards state
diff --git a/stan/math/prim/prob/hypergeometric_rng.hpp b/stan/math/prim/prob/hypergeometric_rng.hpp
index b6d7e410f46..96238b6d5fa 100644
--- a/stan/math/prim/prob/hypergeometric_rng.hpp
+++ b/stan/math/prim/prob/hypergeometric_rng.hpp
@@ -12,7 +12,6 @@ namespace math {
 
 template <class RNG>
 inline int hypergeometric_rng(int N, int a, int b, RNG& rng) {
-  using boost::variate_generator;
   using boost::math::hypergeometric_distribution;
   static constexpr const char* function = "hypergeometric_rng";
   check_bounded(function, "Draws parameter", value_of(N), 0, a + b);
diff --git a/stan/math/prim/prob/inv_chi_square_rng.hpp b/stan/math/prim/prob/inv_chi_square_rng.hpp
index a15c23d323e..0166bf918fa 100644
--- a/stan/math/prim/prob/inv_chi_square_rng.hpp
+++ b/stan/math/prim/prob/inv_chi_square_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size.hpp>
-#include <boost/random/chi_squared_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -27,8 +26,6 @@ namespace math {
 template <typename T_deg, class RNG>
 inline typename VectorBuilder<true, double, T_deg>::type inv_chi_square_rng(
     const T_deg& nu, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::chi_squared_distribution;
   using T_nu_ref = ref_type_t<T_deg>;
   static constexpr const char* function = "inv_chi_square_rng";
   T_nu_ref nu_ref = nu;
@@ -39,9 +36,8 @@ inline typename VectorBuilder<true, double, T_deg>::type inv_chi_square_rng(
   VectorBuilder<true, double, T_deg> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, chi_squared_distribution<> > chi_square_rng(
-        rng, chi_squared_distribution<>(nu_vec[n]));
-    output[n] = 1 / chi_square_rng();
+    std::chi_squared_distribution<> chi_square_rng(nu_vec[n]);
+    output[n] = 1 / chi_square_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/inv_gamma_rng.hpp b/stan/math/prim/prob/inv_gamma_rng.hpp
index cd7baeb7bc1..fb0423172bb 100644
--- a/stan/math/prim/prob/inv_gamma_rng.hpp
+++ b/stan/math/prim/prob/inv_gamma_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
-#include <boost/random/gamma_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -32,8 +31,6 @@ namespace math {
 template <typename T_shape, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_shape, T_scale>::type
 inv_gamma_rng(const T_shape& alpha, const T_scale& beta, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::gamma_distribution;
   using T_alpha_ref = ref_type_t<T_shape>;
   using T_beta_ref = ref_type_t<T_scale>;
   static constexpr const char* function = "inv_gamma_rng";
@@ -50,10 +47,9 @@ inv_gamma_rng(const T_shape& alpha, const T_scale& beta, RNG& rng) {
   VectorBuilder<true, double, T_shape, T_scale> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, gamma_distribution<> > gamma_rng(
-        rng, gamma_distribution<>(alpha_vec[n],
-                                  1 / static_cast<double>(beta_vec[n])));
-    output[n] = 1 / gamma_rng();
+    std::gamma_distribution<> gamma_rng(alpha_vec[n],
+                                        1 / static_cast<double>(beta_vec[n]));
+    output[n] = 1 / gamma_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/logistic_rng.hpp b/stan/math/prim/prob/logistic_rng.hpp
index 9b7e94f1484..71667ff1b61 100644
--- a/stan/math/prim/prob/logistic_rng.hpp
+++ b/stan/math/prim/prob/logistic_rng.hpp
@@ -5,9 +5,8 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/exponential_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -33,8 +32,6 @@ namespace math {
 template <typename T_loc, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_scale>::type logistic_rng(
     const T_loc& mu, const T_scale& sigma, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::exponential_distribution;
   using T_mu_ref = ref_type_t<T_loc>;
   using T_sigma_ref = ref_type_t<T_scale>;
   static constexpr const char* function = "logistic_rng";
@@ -50,10 +47,10 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type logistic_rng(
   size_t N = max_size(mu, sigma);
   VectorBuilder<true, double, T_loc, T_scale> output(N);
 
-  variate_generator<RNG&, exponential_distribution<> > exp_rng(
-      rng, exponential_distribution<>(1));
+  std::exponential_distribution<> exp_rng(1);
   for (size_t n = 0; n < N; ++n) {
-    output[n] = mu_vec[n] - sigma_vec[n] * std::log(exp_rng() / exp_rng());
+    output[n]
+        = mu_vec[n] - sigma_vec[n] * std::log(exp_rng(rng) / exp_rng(rng));
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/loglogistic_rng.hpp b/stan/math/prim/prob/loglogistic_rng.hpp
index fa92d93826e..db04934b5a8 100644
--- a/stan/math/prim/prob/loglogistic_rng.hpp
+++ b/stan/math/prim/prob/loglogistic_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
-#include <boost/random/uniform_01.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -32,8 +31,6 @@ namespace math {
 template <typename T_scale, typename T_shape, class RNG>
 inline typename VectorBuilder<true, double, T_scale, T_shape>::type
 loglogistic_rng(const T_scale& alpha, const T_shape& beta, RNG& rng) {
-  using boost::uniform_01;
-  using boost::variate_generator;
   using std::pow;
   using T_alpha_ref = ref_type_t<T_scale>;
   using T_beta_ref = ref_type_t<T_shape>;
@@ -51,8 +48,8 @@ loglogistic_rng(const T_scale& alpha, const T_shape& beta, RNG& rng) {
   VectorBuilder<true, double, T_scale, T_shape> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, uniform_01<> > uniform01_rng(rng, uniform_01<>());
-    const double tmp = uniform01_rng();
+    std::uniform_real_distribution<> uniform01_rng(0, 1);
+    const double tmp = uniform01_rng(rng);
     output[n] = alpha_vec[n] * pow(tmp / (1 - tmp), 1 / beta_vec[n]);
   }
   return output.data();
diff --git a/stan/math/prim/prob/lognormal_rng.hpp b/stan/math/prim/prob/lognormal_rng.hpp
index 64adc0516c5..583d4a83bc1 100644
--- a/stan/math/prim/prob/lognormal_rng.hpp
+++ b/stan/math/prim/prob/lognormal_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/lognormal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -32,8 +31,6 @@ namespace math {
 template <typename T_loc, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_scale>::type lognormal_rng(
     const T_loc& mu, const T_scale& sigma, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::lognormal_distribution;
   using T_mu_ref = ref_type_t<T_loc>;
   using T_sigma_ref = ref_type_t<T_scale>;
   static constexpr const char* function = "lognormal_rng";
@@ -50,9 +47,8 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type lognormal_rng(
   VectorBuilder<true, double, T_loc, T_scale> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, lognormal_distribution<> > lognorm_rng(
-        rng, lognormal_distribution<>(mu_vec[n], sigma_vec[n]));
-    output[n] = lognorm_rng();
+    std::lognormal_distribution<> lognorm_rng(mu_vec[n], sigma_vec[n]);
+    output[n] = lognorm_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/matrix_normal_prec_rng.hpp b/stan/math/prim/prob/matrix_normal_prec_rng.hpp
index e7a1d275fb2..ae478bd6c18 100644
--- a/stan/math/prim/prob/matrix_normal_prec_rng.hpp
+++ b/stan/math/prim/prob/matrix_normal_prec_rng.hpp
@@ -4,8 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -33,8 +32,6 @@ inline Eigen::MatrixXd matrix_normal_prec_rng(const Eigen::MatrixXd &Mu,
                                               const Eigen::MatrixXd &Sigma,
                                               const Eigen::MatrixXd &D,
                                               RNG &rng) {
-  using boost::normal_distribution;
-  using boost::variate_generator;
   static constexpr const char *function = "matrix_normal_prec_rng";
   check_positive(function, "Sigma rows", Sigma.rows());
   check_finite(function, "Sigma", Sigma);
@@ -79,8 +76,7 @@ inline Eigen::MatrixXd matrix_normal_prec_rng(const Eigen::MatrixXd &Mu,
   int m = Sigma.rows();
   int n = D.rows();
 
-  variate_generator<RNG &, normal_distribution<>> std_normal_rng(
-      rng, normal_distribution<>(0, 1));
+  std::normal_distribution<> std_normal_rng(0, 1);
 
   // X = sqrt[DS]^(-1) C sqrt[DD]^(-1)
   // X ~ N[0, DS, DD]
@@ -93,7 +89,7 @@ inline Eigen::MatrixXd matrix_normal_prec_rng(const Eigen::MatrixXd &Mu,
     for (int row = 0; row < m; ++row) {
       double stddev = row_stddev(row) * col_stddev(col);
       // C(row, col) = std_normal_rng();
-      X(row, col) = stddev * std_normal_rng();
+      X(row, col) = stddev * std_normal_rng(rng);
     }
   }
 
diff --git a/stan/math/prim/prob/multi_normal_cholesky_rng.hpp b/stan/math/prim/prob/multi_normal_cholesky_rng.hpp
index fb1e490c683..7a29eb7a32c 100644
--- a/stan/math/prim/prob/multi_normal_cholesky_rng.hpp
+++ b/stan/math/prim/prob/multi_normal_cholesky_rng.hpp
@@ -7,8 +7,7 @@
 #include <stan/math/prim/fun/size_mvt.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/vector_seq_view.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -34,9 +33,6 @@ inline typename StdVectorBuilder<true, Eigen::VectorXd, T_loc>::type
 multi_normal_cholesky_rng(
     const T_loc& mu,
     const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& L, RNG& rng) {
-  using boost::normal_distribution;
-  using boost::variate_generator;
-
   static constexpr const char* function = "multi_normal_cholesky_rng";
   vector_seq_view<T_loc> mu_vec(mu);
   size_t size_mu = mu_vec[0].size();
@@ -61,13 +57,12 @@ multi_normal_cholesky_rng(
 
   StdVectorBuilder<true, Eigen::VectorXd, T_loc> output(N);
 
-  variate_generator<RNG&, normal_distribution<> > std_normal_rng(
-      rng, normal_distribution<>(0, 1));
+  std::normal_distribution<> std_normal_rng(0, 1);
 
   for (size_t n = 0; n < N; ++n) {
     Eigen::VectorXd z(L.cols());
     for (int i = 0; i < L.cols(); i++) {
-      z(i) = std_normal_rng();
+      z(i) = std_normal_rng(rng);
     }
 
     output[n] = as_column_vector_or_scalar(mu_vec[n]) + L_ref * z;
diff --git a/stan/math/prim/prob/multi_normal_prec_rng.hpp b/stan/math/prim/prob/multi_normal_prec_rng.hpp
index b3c7c63f518..f323f3a1674 100644
--- a/stan/math/prim/prob/multi_normal_prec_rng.hpp
+++ b/stan/math/prim/prob/multi_normal_prec_rng.hpp
@@ -7,8 +7,7 @@
 #include <stan/math/prim/fun/as_column_vector_or_scalar.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/vector_seq_view.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -32,8 +31,6 @@ namespace math {
 template <typename T_loc, class RNG>
 inline typename StdVectorBuilder<true, Eigen::VectorXd, T_loc>::type
 multi_normal_prec_rng(const T_loc &mu, const Eigen::MatrixXd &S, RNG &rng) {
-  using boost::normal_distribution;
-  using boost::variate_generator;
   static constexpr const char *function = "multi_normal_prec_rng";
   check_positive(function, "Precision matrix rows", S.rows());
 
@@ -67,13 +64,12 @@ multi_normal_prec_rng(const T_loc &mu, const Eigen::MatrixXd &S, RNG &rng) {
 
   StdVectorBuilder<true, Eigen::VectorXd, T_loc> output(N);
 
-  variate_generator<RNG &, normal_distribution<>> std_normal_rng(
-      rng, normal_distribution<>(0, 1));
+  std::normal_distribution<> std_normal_rng(0, 1);
 
   for (size_t n = 0; n < N; ++n) {
     Eigen::VectorXd z(S.cols());
     for (int i = 0; i < S.cols(); i++) {
-      z(i) = std_normal_rng();
+      z(i) = std_normal_rng(rng);
     }
 
     output[n]
diff --git a/stan/math/prim/prob/multi_normal_rng.hpp b/stan/math/prim/prob/multi_normal_rng.hpp
index b30d5e1ae3a..0d006a3fc58 100644
--- a/stan/math/prim/prob/multi_normal_rng.hpp
+++ b/stan/math/prim/prob/multi_normal_rng.hpp
@@ -7,8 +7,7 @@
 #include <stan/math/prim/fun/size_mvt.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/vector_seq_view.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -34,8 +33,6 @@ inline typename StdVectorBuilder<true, Eigen::VectorXd, T_loc>::type
 multi_normal_rng(const T_loc& mu,
                  const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& S,
                  RNG& rng) {
-  using boost::normal_distribution;
-  using boost::variate_generator;
   static constexpr const char* function = "multi_normal_rng";
   check_positive(function, "Covariance matrix rows", S.rows());
 
@@ -65,13 +62,12 @@ multi_normal_rng(const T_loc& mu,
 
   StdVectorBuilder<true, Eigen::VectorXd, T_loc> output(N);
 
-  variate_generator<RNG&, normal_distribution<> > std_normal_rng(
-      rng, normal_distribution<>(0, 1));
+  std::normal_distribution<> std_normal_rng(0, 1);
 
   for (size_t n = 0; n < N; ++n) {
     Eigen::VectorXd z(S.cols());
     for (int i = 0; i < S.cols(); i++) {
-      z(i) = std_normal_rng();
+      z(i) = std_normal_rng(rng);
     }
 
     output[n] = as_column_vector_or_scalar(mu_vec[n]) + llt_of_S.matrixL() * z;
diff --git a/stan/math/prim/prob/multi_student_t_cholesky_rng.hpp b/stan/math/prim/prob/multi_student_t_cholesky_rng.hpp
index 074f14f337a..f935679d715 100644
--- a/stan/math/prim/prob/multi_student_t_cholesky_rng.hpp
+++ b/stan/math/prim/prob/multi_student_t_cholesky_rng.hpp
@@ -8,9 +8,8 @@
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/vector_seq_view.hpp>
 #include <stan/math/prim/prob/inv_gamma_rng.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -41,10 +40,6 @@ template <typename T_loc, class RNG>
 inline typename StdVectorBuilder<true, Eigen::VectorXd, T_loc>::type
 multi_student_t_cholesky_rng(double nu, const T_loc& mu,
                              const Eigen::MatrixXd& L, RNG& rng) {
-  using boost::normal_distribution;
-  using boost::variate_generator;
-  using boost::random::gamma_distribution;
-
   static constexpr const char* function = "multi_student_t_cholesky_rng";
   check_not_nan(function, "Degrees of freedom parameter", nu);
   check_positive(function, "Degrees of freedom parameter", nu);
@@ -73,14 +68,13 @@ multi_student_t_cholesky_rng(double nu, const T_loc& mu,
 
   StdVectorBuilder<true, Eigen::VectorXd, T_loc> output(N);
 
-  variate_generator<RNG&, normal_distribution<> > std_normal_rng(
-      rng, normal_distribution<>(0, 1));
+  std::normal_distribution<> std_normal_rng(0, 1);
 
   double w = inv_gamma_rng(nu / 2, nu / 2, rng);
   for (size_t n = 0; n < N; ++n) {
     Eigen::VectorXd z(L.cols());
     for (int i = 0; i < L.cols(); i++) {
-      z(i) = std_normal_rng();
+      z(i) = std_normal_rng(rng);
     }
     z *= std::sqrt(w);
     output[n] = as_column_vector_or_scalar(mu_vec[n]) + L * z;
diff --git a/stan/math/prim/prob/multi_student_t_rng.hpp b/stan/math/prim/prob/multi_student_t_rng.hpp
index cffb45338eb..ef0def652dc 100644
--- a/stan/math/prim/prob/multi_student_t_rng.hpp
+++ b/stan/math/prim/prob/multi_student_t_rng.hpp
@@ -8,9 +8,8 @@
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/vector_seq_view.hpp>
 #include <stan/math/prim/prob/inv_gamma_rng.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -42,10 +41,6 @@ inline typename StdVectorBuilder<true, Eigen::VectorXd, T_loc>::type
 multi_student_t_rng(
     double nu, const T_loc& mu,
     const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& S, RNG& rng) {
-  using boost::normal_distribution;
-  using boost::variate_generator;
-  using boost::random::gamma_distribution;
-
   static constexpr const char* function = "multi_student_t_rng";
   check_not_nan(function, "Degrees of freedom parameter", nu);
   check_positive(function, "Degrees of freedom parameter", nu);
@@ -78,14 +73,13 @@ multi_student_t_rng(
 
   StdVectorBuilder<true, Eigen::VectorXd, T_loc> output(N);
 
-  variate_generator<RNG&, normal_distribution<> > std_normal_rng(
-      rng, normal_distribution<>(0, 1));
+  std::normal_distribution<> std_normal_rng(0, 1);
 
   double w = inv_gamma_rng(nu / 2, nu / 2, rng);
   for (size_t n = 0; n < N; ++n) {
     Eigen::VectorXd z(S.cols());
     for (int i = 0; i < S.cols(); i++) {
-      z(i) = std_normal_rng();
+      z(i) = std_normal_rng(rng);
     }
     z *= std::sqrt(w);
     output[n] = as_column_vector_or_scalar(mu_vec[n]) + llt_of_S.matrixL() * z;
diff --git a/stan/math/prim/prob/neg_binomial_2_log_rng.hpp b/stan/math/prim/prob/neg_binomial_2_log_rng.hpp
index aefce01d132..7806f261d93 100644
--- a/stan/math/prim/prob/neg_binomial_2_log_rng.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_log_rng.hpp
@@ -6,10 +6,8 @@
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/gamma_distribution.hpp>
-#include <boost/random/poisson_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -35,9 +33,6 @@ namespace math {
 template <typename T_loc, typename T_inv, class RNG>
 inline typename VectorBuilder<true, int, T_loc, T_inv>::type
 neg_binomial_2_log_rng(const T_loc& eta, const T_inv& phi, RNG& rng) {
-  using boost::gamma_distribution;
-  using boost::variate_generator;
-  using boost::random::poisson_distribution;
   using T_eta_ref = ref_type_t<T_loc>;
   using T_phi_ref = ref_type_t<T_inv>;
   static constexpr const char* function = "neg_binomial_2_log_rng";
@@ -63,8 +58,8 @@ neg_binomial_2_log_rng(const T_loc& eta, const T_inv& phi, RNG& rng) {
                           "divided by the precision parameter",
                           exp_eta_div_phi);
 
-    double rng_from_gamma = variate_generator<RNG&, gamma_distribution<> >(
-        rng, gamma_distribution<>(phi_vec[n], exp_eta_div_phi))();
+    double rng_from_gamma
+        = std::gamma_distribution<>(phi_vec[n], exp_eta_div_phi)(rng);
 
     // same as the constraints for poisson_rng
     check_less(function, "Random number that came from gamma distribution",
@@ -75,8 +70,7 @@ neg_binomial_2_log_rng(const T_loc& eta, const T_inv& phi, RNG& rng) {
                       "Random number that came from gamma distribution",
                       rng_from_gamma);
 
-    output[n] = variate_generator<RNG&, poisson_distribution<> >(
-        rng, poisson_distribution<>(rng_from_gamma))();
+    output[n] = std::poisson_distribution<>(rng_from_gamma)(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/neg_binomial_2_rng.hpp b/stan/math/prim/prob/neg_binomial_2_rng.hpp
index 67f6b69687d..3d2f1d6a164 100644
--- a/stan/math/prim/prob/neg_binomial_2_rng.hpp
+++ b/stan/math/prim/prob/neg_binomial_2_rng.hpp
@@ -6,9 +6,7 @@
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/gamma_distribution.hpp>
-#include <boost/random/poisson_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -34,9 +32,6 @@ namespace math {
 template <typename T_loc, typename T_prec, class RNG>
 inline typename VectorBuilder<true, int, T_loc, T_prec>::type
 neg_binomial_2_rng(const T_loc& mu, const T_prec& phi, RNG& rng) {
-  using boost::gamma_distribution;
-  using boost::variate_generator;
-  using boost::random::poisson_distribution;
   using T_mu_ref = ref_type_t<T_loc>;
   using T_phi_ref = ref_type_t<T_prec>;
   static constexpr const char* function = "neg_binomial_2_rng";
@@ -61,8 +56,8 @@ neg_binomial_2_rng(const T_loc& mu, const T_prec& phi, RNG& rng) {
                           "precision parameter",
                           mu_div_phi);
 
-    double rng_from_gamma = variate_generator<RNG&, gamma_distribution<> >(
-        rng, gamma_distribution<>(phi_vec[n], mu_div_phi))();
+    double rng_from_gamma
+        = std::gamma_distribution<>(phi_vec[n], mu_div_phi)(rng);
 
     // same as the constraints for poisson_rng
     check_less(function, "Random number that came from gamma distribution",
@@ -73,8 +68,7 @@ neg_binomial_2_rng(const T_loc& mu, const T_prec& phi, RNG& rng) {
                       "Random number that came from gamma distribution",
                       rng_from_gamma);
 
-    output[n] = variate_generator<RNG&, poisson_distribution<> >(
-        rng, poisson_distribution<>(rng_from_gamma))();
+    output[n] = std::poisson_distribution<>(rng_from_gamma)(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/neg_binomial_rng.hpp b/stan/math/prim/prob/neg_binomial_rng.hpp
index 0e01f6693c8..7fcb21ea1c3 100644
--- a/stan/math/prim/prob/neg_binomial_rng.hpp
+++ b/stan/math/prim/prob/neg_binomial_rng.hpp
@@ -6,9 +6,7 @@
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/gamma_distribution.hpp>
-#include <boost/random/poisson_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -34,9 +32,6 @@ namespace math {
 template <typename T_shape, typename T_inv, class RNG>
 inline typename VectorBuilder<true, int, T_shape, T_inv>::type neg_binomial_rng(
     const T_shape& alpha, const T_inv& beta, RNG& rng) {
-  using boost::gamma_distribution;
-  using boost::variate_generator;
-  using boost::random::poisson_distribution;
   using T_alpha_ref = ref_type_t<T_shape>;
   using T_beta_ref = ref_type_t<T_inv>;
   static constexpr const char* function = "neg_binomial_rng";
@@ -53,8 +48,8 @@ inline typename VectorBuilder<true, int, T_shape, T_inv>::type neg_binomial_rng(
   VectorBuilder<true, int, T_shape, T_inv> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    double rng_from_gamma = variate_generator<RNG&, gamma_distribution<> >(
-        rng, gamma_distribution<>(alpha_vec[n], 1.0 / beta_vec[n]))();
+    double rng_from_gamma
+        = std::gamma_distribution<>(alpha_vec[n], 1.0 / beta_vec[n])(rng);
 
     // same as the constraints for poisson_rng
     check_less(function, "Random number that came from gamma distribution",
@@ -65,8 +60,7 @@ inline typename VectorBuilder<true, int, T_shape, T_inv>::type neg_binomial_rng(
                       "Random number that came from gamma distribution",
                       rng_from_gamma);
 
-    output[n] = variate_generator<RNG&, poisson_distribution<> >(
-        rng, poisson_distribution<>(rng_from_gamma))();
+    output[n] = std::poisson_distribution<>(rng_from_gamma)(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/normal_rng.hpp b/stan/math/prim/prob/normal_rng.hpp
index 73c328eb9e4..e106e20e3a6 100644
--- a/stan/math/prim/prob/normal_rng.hpp
+++ b/stan/math/prim/prob/normal_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -32,8 +31,6 @@ namespace math {
 template <typename T_loc, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_scale>::type normal_rng(
     const T_loc& mu, const T_scale& sigma, RNG& rng) {
-  using boost::normal_distribution;
-  using boost::variate_generator;
   using T_mu_ref = ref_type_t<T_loc>;
   using T_sigma_ref = ref_type_t<T_scale>;
   static constexpr const char* function = "normal_rng";
@@ -50,9 +47,8 @@ inline typename VectorBuilder<true, double, T_loc, T_scale>::type normal_rng(
   VectorBuilder<true, double, T_loc, T_scale> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, normal_distribution<> > norm_rng(
-        rng, normal_distribution<>(mu_vec[n], sigma_vec[n]));
-    output[n] = norm_rng();
+    std::normal_distribution<> norm_rng(mu_vec[n], sigma_vec[n]);
+    output[n] = norm_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/ordered_logistic_rng.hpp b/stan/math/prim/prob/ordered_logistic_rng.hpp
index 59921364a26..ad1376fda83 100644
--- a/stan/math/prim/prob/ordered_logistic_rng.hpp
+++ b/stan/math/prim/prob/ordered_logistic_rng.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/inv_logit.hpp>
 #include <stan/math/prim/prob/categorical_rng.hpp>
-#include <boost/random/variate_generator.hpp>
 
 namespace stan {
 namespace math {
@@ -13,7 +12,6 @@ namespace math {
 template <class RNG>
 inline int ordered_logistic_rng(
     double eta, const Eigen::Matrix<double, Eigen::Dynamic, 1>& c, RNG& rng) {
-  using boost::variate_generator;
   static constexpr const char* function = "ordered_logistic";
   check_finite(function, "Location parameter", eta);
   check_greater(function, "Size of cut points parameter", c.size(), 0);
diff --git a/stan/math/prim/prob/pareto_rng.hpp b/stan/math/prim/prob/pareto_rng.hpp
index aac348d7c0e..f3ffe795d98 100644
--- a/stan/math/prim/prob/pareto_rng.hpp
+++ b/stan/math/prim/prob/pareto_rng.hpp
@@ -6,9 +6,8 @@
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
-#include <boost/random/exponential_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -34,8 +33,6 @@ namespace math {
 template <typename T_shape, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_shape, T_scale>::type pareto_rng(
     const T_scale& y_min, const T_shape& alpha, RNG& rng) {
-  using boost::exponential_distribution;
-  using boost::variate_generator;
   static constexpr const char* function = "pareto_rng";
   check_consistent_sizes(function, "Scale Parameter", y_min, "Shape parameter",
                          alpha);
@@ -50,9 +47,8 @@ inline typename VectorBuilder<true, double, T_shape, T_scale>::type pareto_rng(
   VectorBuilder<true, double, T_scale, T_shape> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, exponential_distribution<> > exp_rng(
-        rng, exponential_distribution<>(alpha_vec[n]));
-    output[n] = y_min_vec[n] * std::exp(exp_rng());
+    std::exponential_distribution<> exp_rng(alpha_vec[n]);
+    output[n] = y_min_vec[n] * std::exp(exp_rng(rng));
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/pareto_type_2_rng.hpp b/stan/math/prim/prob/pareto_type_2_rng.hpp
index e2bcd5f0db8..874d72810c3 100644
--- a/stan/math/prim/prob/pareto_type_2_rng.hpp
+++ b/stan/math/prim/prob/pareto_type_2_rng.hpp
@@ -8,9 +8,8 @@
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/prob/exponential_rng.hpp>
 #include <stan/math/prim/prob/normal_rng.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -41,8 +40,6 @@ template <typename T_loc, typename T_scale, typename T_shape, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_scale, T_shape>::type
 pareto_type_2_rng(const T_loc& mu, const T_scale& lambda, const T_shape& alpha,
                   RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::uniform_real_distribution;
   static constexpr const char* function = "pareto_type_2_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          lambda, "Shape Parameter", alpha);
@@ -59,10 +56,9 @@ pareto_type_2_rng(const T_loc& mu, const T_scale& lambda, const T_shape& alpha,
   size_t N = max_size(mu, lambda, alpha);
   VectorBuilder<true, double, T_loc, T_scale, T_shape> output(N);
 
-  variate_generator<RNG&, uniform_real_distribution<> > uniform_rng(
-      rng, uniform_real_distribution<>(0.0, 1.0));
+  std::uniform_real_distribution<> uniform_rng(0.0, 1.0);
   for (size_t n = 0; n < N; ++n) {
-    output[n] = (std::pow(1.0 - uniform_rng(), -1.0 / alpha_vec[n]) - 1.0)
+    output[n] = (std::pow(1.0 - uniform_rng(rng), -1.0 / alpha_vec[n]) - 1.0)
                     * lambda_vec[n]
                 + mu_vec[n];
   }
diff --git a/stan/math/prim/prob/poisson_binomial_rng.hpp b/stan/math/prim/prob/poisson_binomial_rng.hpp
index e05bedb07a4..15ac8bb16e6 100644
--- a/stan/math/prim/prob/poisson_binomial_rng.hpp
+++ b/stan/math/prim/prob/poisson_binomial_rng.hpp
@@ -4,8 +4,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/size.hpp>
-#include <boost/random/bernoulli_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -30,9 +29,8 @@ inline int poisson_binomial_rng(const T_theta& theta, RNG& rng) {
 
   int y = 0;
   for (size_t i = 0; i < theta.size(); ++i) {
-    boost::variate_generator<RNG&, boost::bernoulli_distribution<> >
-        bernoulli_rng(rng, boost::bernoulli_distribution<>(theta_ref[i]));
-    y += bernoulli_rng();
+    std::bernoulli_distribution bernoulli_rng(theta_ref[i]);
+    y += bernoulli_rng(rng);
   }
 
   return y;
diff --git a/stan/math/prim/prob/poisson_log_rng.hpp b/stan/math/prim/prob/poisson_log_rng.hpp
index 262bccb25b2..fa2b5643d95 100644
--- a/stan/math/prim/prob/poisson_log_rng.hpp
+++ b/stan/math/prim/prob/poisson_log_rng.hpp
@@ -7,9 +7,8 @@
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
-#include <boost/random/poisson_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -30,8 +29,6 @@ namespace math {
 template <typename T_rate, class RNG>
 inline typename VectorBuilder<true, int, T_rate>::type poisson_log_rng(
     const T_rate& alpha, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::poisson_distribution;
   static constexpr const char* function = "poisson_log_rng";
   static constexpr double POISSON_MAX_LOG_RATE = 30 * LOG_TWO;
   const auto& alpha_ref = to_ref(alpha);
@@ -43,9 +40,8 @@ inline typename VectorBuilder<true, int, T_rate>::type poisson_log_rng(
   VectorBuilder<true, int, T_rate> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, poisson_distribution<> > poisson_rng(
-        rng, poisson_distribution<>(std::exp(alpha_vec[n])));
-    output[n] = poisson_rng();
+    std::poisson_distribution<> poisson_rng(std::exp(alpha_vec[n]));
+    output[n] = poisson_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/poisson_rng.hpp b/stan/math/prim/prob/poisson_rng.hpp
index 83052699ee2..8dbe124664b 100644
--- a/stan/math/prim/prob/poisson_rng.hpp
+++ b/stan/math/prim/prob/poisson_rng.hpp
@@ -7,8 +7,7 @@
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
-#include <boost/random/poisson_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -29,8 +28,6 @@ namespace math {
 template <typename T_rate, class RNG>
 inline typename VectorBuilder<true, int, T_rate>::type poisson_rng(
     const T_rate& lambda, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::poisson_distribution;
   static constexpr const char* function = "poisson_rng";
   const auto& lambda_ref = to_ref(lambda);
   check_positive(function, "Rate parameter", lambda_ref);
@@ -41,9 +38,8 @@ inline typename VectorBuilder<true, int, T_rate>::type poisson_rng(
   VectorBuilder<true, int, T_rate> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, poisson_distribution<> > poisson_rng(
-        rng, poisson_distribution<>(lambda_vec[n]));
-    output[n] = poisson_rng();
+    std::poisson_distribution<> poisson_rng(lambda_vec[n]);
+    output[n] = poisson_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/rayleigh_rng.hpp b/stan/math/prim/prob/rayleigh_rng.hpp
index 55b382be7fc..0547e8d6579 100644
--- a/stan/math/prim/prob/rayleigh_rng.hpp
+++ b/stan/math/prim/prob/rayleigh_rng.hpp
@@ -6,9 +6,8 @@
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -29,8 +28,6 @@ namespace math {
 template <typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_scale>::type rayleigh_rng(
     const T_scale& sigma, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::uniform_real_distribution;
   static constexpr const char* function = "rayleigh_rng";
   const auto& sigma_ref = to_ref(sigma);
   check_positive_finite(function, "Scale parameter", sigma_ref);
@@ -39,10 +36,9 @@ inline typename VectorBuilder<true, double, T_scale>::type rayleigh_rng(
   size_t N = stan::math::size(sigma);
   VectorBuilder<true, double, T_scale> output(N);
 
-  variate_generator<RNG&, uniform_real_distribution<> > uniform_rng(
-      rng, uniform_real_distribution<>(0.0, 1.0));
+  std::uniform_real_distribution<> uniform_rng(0.0, 1.0);
   for (size_t n = 0; n < N; ++n) {
-    output[n] = sigma_vec[n] * std::sqrt(-2.0 * std::log(uniform_rng()));
+    output[n] = sigma_vec[n] * std::sqrt(-2.0 * std::log(uniform_rng(rng)));
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp b/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp
index 819255a94a1..40f38b0e011 100644
--- a/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp
+++ b/stan/math/prim/prob/scaled_inv_chi_square_rng.hpp
@@ -6,8 +6,7 @@
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
-#include <boost/random/chi_squared_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -34,8 +33,6 @@ namespace math {
 template <typename T_deg, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_deg, T_scale>::type
 scaled_inv_chi_square_rng(const T_deg& nu, const T_scale& s, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::chi_squared_distribution;
   static constexpr const char* function = "scaled_inv_chi_square_rng";
   check_consistent_sizes(function, "Location parameter", nu, "Scale Parameter",
                          s);
@@ -50,9 +47,8 @@ scaled_inv_chi_square_rng(const T_deg& nu, const T_scale& s, RNG& rng) {
   VectorBuilder<true, double, T_deg, T_scale> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, chi_squared_distribution<> > chi_square_rng(
-        rng, chi_squared_distribution<>(nu_vec[n]));
-    output[n] = nu_vec[n] * s_vec[n] * s_vec[n] / chi_square_rng();
+    std::chi_squared_distribution<> chi_square_rng(nu_vec[n]);
+    output[n] = nu_vec[n] * s_vec[n] * s_vec[n] / chi_square_rng(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/skew_double_exponential_rng.hpp b/stan/math/prim/prob/skew_double_exponential_rng.hpp
index 7b238d71189..6e90c69bc5b 100644
--- a/stan/math/prim/prob/skew_double_exponential_rng.hpp
+++ b/stan/math/prim/prob/skew_double_exponential_rng.hpp
@@ -5,9 +5,8 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/log1m.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -37,8 +36,6 @@ template <typename T_loc, typename T_scale, typename T_skewness, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_scale, T_skewness>::type
 skew_double_exponential_rng(const T_loc& mu, const T_scale& sigma,
                             const T_skewness& tau, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::uniform_real_distribution;
   using T_mu_ref = ref_type_t<T_loc>;
   using T_sigma_ref = ref_type_t<T_scale>;
   using T_tau_ref = ref_type_t<T_skewness>;
@@ -58,10 +55,9 @@ skew_double_exponential_rng(const T_loc& mu, const T_scale& sigma,
   size_t N = max_size(mu, sigma, tau);
   VectorBuilder<true, double, T_loc, T_scale, T_skewness> output(N);
 
-  variate_generator<RNG&, uniform_real_distribution<> > z_rng(
-      rng, uniform_real_distribution<>(0.0, 1.0));
+  std::uniform_real_distribution<> z_rng(0.0, 1.0);
   for (size_t n = 0; n < N; ++n) {
-    double z = z_rng();
+    double z = z_rng(rng);
     if (z < tau_vec[n]) {
       output[n]
           = log(z / tau_vec[n]) * sigma_vec[n] / (2.0 * (1.0 - tau_vec[n]))
diff --git a/stan/math/prim/prob/skew_normal_rng.hpp b/stan/math/prim/prob/skew_normal_rng.hpp
index d4cc8590de7..3155e551513 100644
--- a/stan/math/prim/prob/skew_normal_rng.hpp
+++ b/stan/math/prim/prob/skew_normal_rng.hpp
@@ -7,8 +7,7 @@
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/prob/normal_rng.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -39,8 +38,6 @@ template <typename T_loc, typename T_scale, typename T_shape, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_scale, T_shape>::type
 skew_normal_rng(const T_loc& mu, const T_scale& sigma, const T_shape& alpha,
                 RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::normal_distribution;
   static constexpr const char* function = "skew_normal_rng";
   check_consistent_sizes(function, "Location parameter", mu, "Scale Parameter",
                          sigma, "Shape Parameter", alpha);
@@ -57,11 +54,10 @@ skew_normal_rng(const T_loc& mu, const T_scale& sigma, const T_shape& alpha,
   size_t N = max_size(mu, sigma, alpha);
   VectorBuilder<true, double, T_loc, T_scale, T_shape> output(N);
 
-  variate_generator<RNG&, normal_distribution<> > norm_rng(
-      rng, normal_distribution<>(0, 1));
+  std::normal_distribution<> norm_rng(0, 1);
   for (size_t n = 0; n < N; ++n) {
-    double r1 = norm_rng();
-    double r2 = norm_rng();
+    double r1 = norm_rng(rng);
+    double r2 = norm_rng(rng);
 
     if (r2 > alpha_vec[n] * r1) {
       r1 = -r1;
diff --git a/stan/math/prim/prob/std_normal_rng.hpp b/stan/math/prim/prob/std_normal_rng.hpp
index 9a5b5d0bbb3..3cf0d3257ea 100644
--- a/stan/math/prim/prob/std_normal_rng.hpp
+++ b/stan/math/prim/prob/std_normal_rng.hpp
@@ -2,8 +2,7 @@
 #define STAN_MATH_PRIM_PROB_STD_NORMAL_RNG_HPP
 
 #include <stan/math/prim/meta.hpp>
-#include <boost/random/normal_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -18,13 +17,9 @@ namespace math {
  */
 template <class RNG>
 inline double std_normal_rng(RNG& rng) {
-  using boost::normal_distribution;
-  using boost::variate_generator;
+  std::normal_distribution<> norm_rng(0, 1);
 
-  variate_generator<RNG&, normal_distribution<>> norm_rng(
-      rng, normal_distribution<>(0, 1));
-
-  return norm_rng();
+  return norm_rng(rng);
 }
 
 }  // namespace math
diff --git a/stan/math/prim/prob/student_t_rng.hpp b/stan/math/prim/prob/student_t_rng.hpp
index ce7d271ef30..275ff170bed 100644
--- a/stan/math/prim/prob/student_t_rng.hpp
+++ b/stan/math/prim/prob/student_t_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/student_t_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -40,8 +39,6 @@ student_t_rng(const T_deg& nu, const T_loc& mu, const T_scale& sigma,
   using T_nu_ref = ref_type_t<T_deg>;
   using T_mu_ref = ref_type_t<T_loc>;
   using T_sigma_ref = ref_type_t<T_scale>;
-  using boost::variate_generator;
-  using boost::random::student_t_distribution;
   static constexpr const char* function = "student_t_rng";
   check_consistent_sizes(function, "Degrees of freedom parameter", nu,
                          "Location parameter", mu, "Scale Parameter", sigma);
@@ -59,9 +56,8 @@ student_t_rng(const T_deg& nu, const T_loc& mu, const T_scale& sigma,
   VectorBuilder<true, double, T_deg, T_loc, T_scale> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, student_t_distribution<> > rng_unit_student_t(
-        rng, student_t_distribution<>(nu_vec[n]));
-    output[n] = mu_vec[n] + sigma_vec[n] * rng_unit_student_t();
+    std::student_t_distribution<> rng_unit_student_t(nu_vec[n]);
+    output[n] = mu_vec[n] + sigma_vec[n] * rng_unit_student_t(rng);
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/uniform_rng.hpp b/stan/math/prim/prob/uniform_rng.hpp
index ad18cf1add3..e1daa12424b 100644
--- a/stan/math/prim/prob/uniform_rng.hpp
+++ b/stan/math/prim/prob/uniform_rng.hpp
@@ -7,8 +7,7 @@
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -36,8 +35,6 @@ inline typename VectorBuilder<true, double, T_alpha, T_beta>::type uniform_rng(
     const T_alpha& alpha, const T_beta& beta, RNG& rng) {
   using T_alpha_ref = ref_type_t<T_alpha>;
   using T_beta_ref = ref_type_t<T_beta>;
-  using boost::variate_generator;
-  using boost::random::uniform_real_distribution;
   static constexpr const char* function = "uniform_rng";
   check_consistent_sizes(function, "Lower bound parameter", alpha,
                          "Upper bound parameter", beta);
@@ -52,10 +49,9 @@ inline typename VectorBuilder<true, double, T_alpha, T_beta>::type uniform_rng(
   size_t N = max_size(alpha, beta);
   VectorBuilder<true, double, T_alpha, T_beta> output(N);
 
-  variate_generator<RNG&, uniform_real_distribution<> > uniform_rng(
-      rng, uniform_real_distribution<>(0.0, 1.0));
+  std::uniform_real_distribution<> uniform_rng(0.0, 1.0);
   for (size_t n = 0; n < N; ++n) {
-    output[n] = (beta_vec[n] - alpha_vec[n]) * uniform_rng() + alpha_vec[n];
+    output[n] = (beta_vec[n] - alpha_vec[n]) * uniform_rng(rng) + alpha_vec[n];
   }
 
   return output.data();
diff --git a/stan/math/prim/prob/von_mises_rng.hpp b/stan/math/prim/prob/von_mises_rng.hpp
index c68a6227f74..f08ae8c773b 100644
--- a/stan/math/prim/prob/von_mises_rng.hpp
+++ b/stan/math/prim/prob/von_mises_rng.hpp
@@ -6,9 +6,8 @@
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
 #include <cmath>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -45,8 +44,6 @@ namespace math {
 template <typename T_loc, typename T_conc, class RNG>
 inline typename VectorBuilder<true, double, T_loc, T_conc>::type von_mises_rng(
     const T_loc& mu, const T_conc& kappa, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::uniform_real_distribution;
   using T_mu_ref = ref_type_t<T_loc>;
   using T_kappa_ref = ref_type_t<T_conc>;
   static constexpr const char* function = "von_mises_rng";
@@ -64,14 +61,13 @@ inline typename VectorBuilder<true, double, T_loc, T_conc>::type von_mises_rng(
   size_t N = max_size(mu, kappa_ref);
   VectorBuilder<true, double, T_loc, T_conc> output(N);
 
-  variate_generator<RNG&, uniform_real_distribution<> > uniform_rng(
-      rng, uniform_real_distribution<>(0.0, 1.0));
+  std::uniform_real_distribution<> uniform_rng(0.0, 1.0);
 
   for (size_t n = 0; n < N; ++n) {
     // for kappa sufficiently close to zero, it reduces to a
     // circular uniform distribution centered at mu
     if (kappa_vec[n] < 1.4e-8) {
-      output[n] = (uniform_rng() - 0.5) * TWO_PI
+      output[n] = (uniform_rng(rng) - 0.5) * TWO_PI
                   + std::fmod(std::fmod(mu_vec[n], TWO_PI) + TWO_PI, TWO_PI);
       continue;
     }
@@ -83,10 +79,10 @@ inline typename VectorBuilder<true, double, T_loc, T_conc>::type von_mises_rng(
     bool done = false;
     double W;
     while (!done) {
-      double Z = std::cos(pi() * uniform_rng());
+      double Z = std::cos(pi() * uniform_rng(rng));
       W = (1 + s * Z) / (s + Z);
       double Y = kappa_vec[n] * (s - W);
-      double U2 = uniform_rng();
+      double U2 = uniform_rng(rng);
       done = Y * (2 - Y) - U2 > 0;
 
       if (!done) {
@@ -94,7 +90,7 @@ inline typename VectorBuilder<true, double, T_loc, T_conc>::type von_mises_rng(
       }
     }
 
-    double U3 = uniform_rng() - 0.5;
+    double U3 = uniform_rng(rng) - 0.5;
     double sign = ((U3 >= 0) - (U3 <= 0));
 
     //  it's really an fmod() with a positivity constraint
diff --git a/stan/math/prim/prob/weibull_rng.hpp b/stan/math/prim/prob/weibull_rng.hpp
index 0626add9fd8..a216e9298b4 100644
--- a/stan/math/prim/prob/weibull_rng.hpp
+++ b/stan/math/prim/prob/weibull_rng.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/max_size.hpp>
 #include <stan/math/prim/fun/scalar_seq_view.hpp>
-#include <boost/random/weibull_distribution.hpp>
-#include <boost/random/variate_generator.hpp>
+#include <random>
 
 namespace stan {
 namespace math {
@@ -32,8 +31,6 @@ namespace math {
 template <typename T_shape, typename T_scale, class RNG>
 inline typename VectorBuilder<true, double, T_shape, T_scale>::type weibull_rng(
     const T_shape& alpha, const T_scale& sigma, RNG& rng) {
-  using boost::variate_generator;
-  using boost::random::weibull_distribution;
   using T_alpha_ref = ref_type_t<T_shape>;
   using T_sigma_ref = ref_type_t<T_scale>;
   static constexpr const char* function = "weibull_rng";
@@ -50,9 +47,8 @@ inline typename VectorBuilder<true, double, T_shape, T_scale>::type weibull_rng(
   VectorBuilder<true, double, T_shape, T_scale> output(N);
 
   for (size_t n = 0; n < N; ++n) {
-    variate_generator<RNG&, weibull_distribution<> > weibull_rng(
-        rng, weibull_distribution<>(alpha_vec[n], sigma_vec[n]));
-    output[n] = weibull_rng();
+    std::weibull_distribution<> weibull_rng(alpha_vec[n], sigma_vec[n]);
+    output[n] = weibull_rng(rng);
   }
 
   return output.data();
diff --git a/test/unit/math/fwd/prob/categorical_test.cpp b/test/unit/math/fwd/prob/categorical_test.cpp
index f7700bc9389..a44f6ce933c 100644
--- a/test/unit/math/fwd/prob/categorical_test.cpp
+++ b/test/unit/math/fwd/prob/categorical_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
diff --git a/test/unit/math/fwd/prob/dirichlet_test.cpp b/test/unit/math/fwd/prob/dirichlet_test.cpp
index fc844a4f2ec..82b4b4a9721 100644
--- a/test/unit/math/fwd/prob/dirichlet_test.cpp
+++ b/test/unit/math/fwd/prob/dirichlet_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributions, fvar_double) {
diff --git a/test/unit/math/fwd/prob/inv_wishart_test.cpp b/test/unit/math/fwd/prob/inv_wishart_test.cpp
index f996bf1a1bd..fe1431efa3a 100644
--- a/test/unit/math/fwd/prob/inv_wishart_test.cpp
+++ b/test/unit/math/fwd/prob/inv_wishart_test.cpp
@@ -1,7 +1,6 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
 
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 
diff --git a/test/unit/math/fwd/prob/lkj_corr_test.cpp b/test/unit/math/fwd/prob/lkj_corr_test.cpp
index 820dd5fdc74..4a92f74a7db 100644
--- a/test/unit/math/fwd/prob/lkj_corr_test.cpp
+++ b/test/unit/math/fwd/prob/lkj_corr_test.cpp
@@ -1,11 +1,11 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsLkjCorr, fvar_double) {
   using stan::math::fvar;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int K = 4;
   Eigen::Matrix<fvar<double>, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -15,7 +15,7 @@ TEST(ProbDistributionsLkjCorr, fvar_double) {
   fvar<double> eta = stan::math::uniform_rng(0, 2, rng);
   fvar<double> f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_, stan::math::lkj_corr_lpdf(Sigma, eta).val_);
-  EXPECT_FLOAT_EQ(2.5177896, stan::math::lkj_corr_lpdf(Sigma, eta).d_);
+  EXPECT_FLOAT_EQ(-2.9161840, stan::math::lkj_corr_lpdf(Sigma, eta).d_);
   eta = 1.0;
   f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_, stan::math::lkj_corr_lpdf(Sigma, eta).val_);
@@ -24,7 +24,7 @@ TEST(ProbDistributionsLkjCorr, fvar_double) {
 
 TEST(ProbDistributionsLkjCorrCholesky, fvar_double) {
   using stan::math::fvar;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int K = 4;
   Eigen::Matrix<fvar<double>, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -34,7 +34,8 @@ TEST(ProbDistributionsLkjCorrCholesky, fvar_double) {
   fvar<double> eta = stan::math::uniform_rng(0, 2, rng);
   fvar<double> f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_, stan::math::lkj_corr_cholesky_lpdf(Sigma, eta).val_);
-  EXPECT_FLOAT_EQ(6.7766843, stan::math::lkj_corr_cholesky_lpdf(Sigma, eta).d_);
+  EXPECT_FLOAT_EQ(-1.3742759,
+                  stan::math::lkj_corr_cholesky_lpdf(Sigma, eta).d_);
   eta = 1.0;
   f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_, stan::math::lkj_corr_cholesky_lpdf(Sigma, eta).val_);
@@ -43,7 +44,7 @@ TEST(ProbDistributionsLkjCorrCholesky, fvar_double) {
 
 TEST(ProbDistributionsLkjCorr, fvar_fvar_double) {
   using stan::math::fvar;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int K = 4;
   Eigen::Matrix<fvar<fvar<double> >, Eigen::Dynamic, Eigen::Dynamic> Sigma(K,
                                                                            K);
@@ -54,7 +55,7 @@ TEST(ProbDistributionsLkjCorr, fvar_fvar_double) {
   fvar<fvar<double> > eta = stan::math::uniform_rng(0, 2, rng);
   fvar<fvar<double> > f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_.val_, stan::math::lkj_corr_lpdf(Sigma, eta).val_.val_);
-  EXPECT_FLOAT_EQ(2.5177896, stan::math::lkj_corr_lpdf(Sigma, eta).d_.val_);
+  EXPECT_FLOAT_EQ(-2.9161840, stan::math::lkj_corr_lpdf(Sigma, eta).d_.val_);
   eta = 1.0;
   f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_.val_, stan::math::lkj_corr_lpdf(Sigma, eta).val_.val_);
@@ -63,7 +64,7 @@ TEST(ProbDistributionsLkjCorr, fvar_fvar_double) {
 
 TEST(ProbDistributionsLkjCorrCholesky, fvar_fvar_double) {
   using stan::math::fvar;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int K = 4;
   Eigen::Matrix<fvar<fvar<double> >, Eigen::Dynamic, Eigen::Dynamic> Sigma(K,
                                                                            K);
@@ -75,7 +76,7 @@ TEST(ProbDistributionsLkjCorrCholesky, fvar_fvar_double) {
   fvar<fvar<double> > f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_.val_,
                   stan::math::lkj_corr_cholesky_lpdf(Sigma, eta).val_.val_);
-  EXPECT_FLOAT_EQ(6.7766843,
+  EXPECT_FLOAT_EQ(-1.3742759,
                   stan::math::lkj_corr_cholesky_lpdf(Sigma, eta).d_.val_);
   eta = 1.0;
   f = stan::math::do_lkj_constant(eta, K);
diff --git a/test/unit/math/fwd/prob/multi_normal_cholesky_test.cpp b/test/unit/math/fwd/prob/multi_normal_cholesky_test.cpp
index 13a0825c06c..b6692ad788e 100644
--- a/test/unit/math/fwd/prob/multi_normal_cholesky_test.cpp
+++ b/test/unit/math/fwd/prob/multi_normal_cholesky_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsMultiNormalCholesky, fvar_double) {
diff --git a/test/unit/math/fwd/prob/multi_normal_test.cpp b/test/unit/math/fwd/prob/multi_normal_test.cpp
index 1f247efcf0e..ca71c30055f 100644
--- a/test/unit/math/fwd/prob/multi_normal_test.cpp
+++ b/test/unit/math/fwd/prob/multi_normal_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsMultiNormal, fvar_double) {
diff --git a/test/unit/math/fwd/prob/multi_student_t_test.cpp b/test/unit/math/fwd/prob/multi_student_t_test.cpp
index f57b8c0609a..45df24ebc67 100644
--- a/test/unit/math/fwd/prob/multi_student_t_test.cpp
+++ b/test/unit/math/fwd/prob/multi_student_t_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsMultiStudentT, fvar_double) {
diff --git a/test/unit/math/fwd/prob/multinomial_test.cpp b/test/unit/math/fwd/prob/multinomial_test.cpp
index ab55eee58a8..67ee06a63a9 100644
--- a/test/unit/math/fwd/prob/multinomial_test.cpp
+++ b/test/unit/math/fwd/prob/multinomial_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
diff --git a/test/unit/math/fwd/prob/wishart_test.cpp b/test/unit/math/fwd/prob/wishart_test.cpp
index ae9bd5ed936..b379834e233 100644
--- a/test/unit/math/fwd/prob/wishart_test.cpp
+++ b/test/unit/math/fwd/prob/wishart_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/fwd.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 #include <boost/math/distributions.hpp>
 
diff --git a/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp b/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp
index 11c92e17770..90f3f9b058b 100644
--- a/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp
+++ b/test/unit/math/laplace/laplace_bernoulli_logit_rng_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math.hpp>
 #include <stan/math/mix.hpp>
 
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 
 #include <gtest/gtest.h>
@@ -74,7 +74,7 @@ TEST(laplace_bernoulli_logit_rng, two_dim_diag) {
   std::vector<double> d0;
   std::vector<int> di0;
   std::vector<Eigen::VectorXd> x_dummy;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   rng.seed(1954);
   Eigen::MatrixXd theta_pred = laplace_latent_bernoulli_logit_rng(
       y, y_index, mean, 1, diagonal_kernel_functor{},
diff --git a/test/unit/math/laplace/laplace_latent_solve_test.cpp b/test/unit/math/laplace/laplace_latent_solve_test.cpp
index 8bb4795dffb..2bb6c9a2570 100644
--- a/test/unit/math/laplace/laplace_latent_solve_test.cpp
+++ b/test/unit/math/laplace/laplace_latent_solve_test.cpp
@@ -2,8 +2,6 @@
 #include <stan/math/mix.hpp>
 #include <test/unit/math/laplace/laplace_utility.hpp>
 
-#include <boost/random/mersenne_twister.hpp>
-
 #include <gtest/gtest.h>
 #include <stdexcept>
 #include <vector>
diff --git a/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp b/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp
index 3e248b8874e..d1179d491d9 100644
--- a/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp
+++ b/test/unit/math/laplace/laplace_neg_binomial_2_log_rng_test.cpp
@@ -2,7 +2,7 @@
 #include <stan/math/mix.hpp>
 #include <test/unit/math/laplace/laplace_utility.hpp>
 
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 
 #include <gtest/gtest.h>
@@ -96,7 +96,7 @@ TEST(laplace_latent_neg_binomial_2_log_rng, count_two_dim_diag) {
       = algebra_solver(stationary_point_nb(), theta_0, phi, d0, di0);
   Eigen::MatrixXd K_laplace = laplace_covariance_nb(theta_root, phi, eta);
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   rng.seed(1954);
   Eigen::MatrixXd theta_benchmark
       = stan::math::multi_normal_rng(theta_root, K_laplace, rng);
diff --git a/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp b/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp
index 6c5a13670f0..d7f0ea5e9f8 100644
--- a/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp
+++ b/test/unit/math/laplace/laplace_poisson_log_rng_test.cpp
@@ -2,7 +2,7 @@
 #include <stan/math/mix.hpp>
 #include <test/unit/math/laplace/laplace_utility.hpp>
 
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 
 #include <gtest/gtest.h>
@@ -20,7 +20,7 @@ TEST_F(laplace_count_two_dim_diag_test, poisson_log_likelihood) {
   Eigen::MatrixXd K_laplace
       = stan::math::test::laplace_covariance(theta_root, phi);
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   rng.seed(1954);
   Eigen::MatrixXd theta_pred = laplace_latent_poisson_log_rng(
       y, y_index, 0, 1, stan::math::test::diagonal_kernel_functor{},
diff --git a/test/unit/math/laplace/laplace_utility.hpp b/test/unit/math/laplace/laplace_utility.hpp
index b7ec38c1619..46f829310db 100644
--- a/test/unit/math/laplace/laplace_utility.hpp
+++ b/test/unit/math/laplace/laplace_utility.hpp
@@ -4,6 +4,7 @@
 #include <test/unit/math/laplace/aki_disease_data/x1.hpp>
 #include <boost/algorithm/string.hpp>
 #include <iostream>
+#include <random>
 #include <gtest/gtest.h>
 
 namespace stan {
@@ -411,7 +412,7 @@ class laplace_count_two_dim_diag_test : public ::testing::Test {
   std::vector<int> di0;
   Eigen::MatrixXd K_laplace;
   Eigen::MatrixXd theta_benchmark;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   double tol;
   int n_sim;
 };
diff --git a/test/unit/math/mix/prob/categorical_test.cpp b/test/unit/math/mix/prob/categorical_test.cpp
index fa54fc94d75..bb1d2b815eb 100644
--- a/test/unit/math/mix/prob/categorical_test.cpp
+++ b/test/unit/math/mix/prob/categorical_test.cpp
@@ -1,7 +1,6 @@
 #include <stan/math/mix.hpp>
 #include <test/unit/math/test_ad.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
diff --git a/test/unit/math/mix/prob/lkj_corr_test.cpp b/test/unit/math/mix/prob/lkj_corr_test.cpp
index bb96388ade9..827f2f1b20d 100644
--- a/test/unit/math/mix/prob/lkj_corr_test.cpp
+++ b/test/unit/math/mix/prob/lkj_corr_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/test_ad.hpp>
 #include <gtest/gtest.h>
 #include <test/unit/math/rev/prob/lkj_corr_cholesky_test_functors.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <test/unit/math/mix/prob/higher_order_utils.hpp>
 #include <vector>
@@ -10,7 +10,7 @@
 TEST_F(AgradRev, ProbDistributionsLkjCorr_fvar_var) {
   using stan::math::fvar;
   using stan::math::var;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int K = 4;
   Eigen::Matrix<fvar<var>, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -21,7 +21,7 @@ TEST_F(AgradRev, ProbDistributionsLkjCorr_fvar_var) {
   fvar<var> f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_.val(),
                   stan::math::lkj_corr_lpdf(Sigma, eta).val_.val());
-  EXPECT_FLOAT_EQ(2.5177896, stan::math::lkj_corr_lpdf(Sigma, eta).d_.val());
+  EXPECT_FLOAT_EQ(-2.9161840, stan::math::lkj_corr_lpdf(Sigma, eta).d_.val());
   eta = 1.0;
   f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_.val(),
@@ -32,7 +32,7 @@ TEST_F(AgradRev, ProbDistributionsLkjCorr_fvar_var) {
 TEST_F(AgradRev, ProbDistributionsLkjCorrCholesky_fvar_var) {
   using stan::math::fvar;
   using stan::math::var;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int K = 4;
   Eigen::Matrix<fvar<var>, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -43,7 +43,7 @@ TEST_F(AgradRev, ProbDistributionsLkjCorrCholesky_fvar_var) {
   fvar<var> f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_.val(),
                   stan::math::lkj_corr_cholesky_lpdf(Sigma, eta).val_.val());
-  EXPECT_FLOAT_EQ(6.7766843,
+  EXPECT_FLOAT_EQ(-1.3742759,
                   stan::math::lkj_corr_cholesky_lpdf(Sigma, eta).d_.val());
   eta = 1.0;
   f = stan::math::do_lkj_constant(eta, K);
@@ -55,7 +55,7 @@ TEST_F(AgradRev, ProbDistributionsLkjCorrCholesky_fvar_var) {
 TEST_F(AgradRev, ProbDistributionsLkjCorr_fvar_fvar_var) {
   using stan::math::fvar;
   using stan::math::var;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int K = 4;
   Eigen::Matrix<fvar<fvar<var> >, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -66,7 +66,7 @@ TEST_F(AgradRev, ProbDistributionsLkjCorr_fvar_fvar_var) {
   fvar<fvar<var> > f = stan::math::do_lkj_constant(eta, K);
   EXPECT_FLOAT_EQ(f.val_.val_.val(),
                   stan::math::lkj_corr_lpdf(Sigma, eta).val_.val_.val());
-  EXPECT_FLOAT_EQ(2.5177896,
+  EXPECT_FLOAT_EQ(-2.9161840,
                   stan::math::lkj_corr_lpdf(Sigma, eta).d_.val_.val());
   eta = 1.0;
   f = stan::math::do_lkj_constant(eta, K);
@@ -79,7 +79,7 @@ TEST_F(AgradRev, ProbDistributionsLkjCorr_fvar_fvar_var) {
 TEST_F(AgradRev, ProbDistributionsLkjCorrCholesky_fvar_fvar_var) {
   using stan::math::fvar;
   using stan::math::var;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int K = 4;
   Eigen::Matrix<fvar<fvar<var> >, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
@@ -91,7 +91,7 @@ TEST_F(AgradRev, ProbDistributionsLkjCorrCholesky_fvar_fvar_var) {
   EXPECT_FLOAT_EQ(
       f.val_.val_.val(),
       stan::math::lkj_corr_cholesky_lpdf(Sigma, eta).val_.val_.val());
-  EXPECT_FLOAT_EQ(6.7766843,
+  EXPECT_FLOAT_EQ(-1.3742759,
                   stan::math::lkj_corr_cholesky_lpdf(Sigma, eta).d_.val_.val());
   eta = 1.0;
   f = stan::math::do_lkj_constant(eta, K);
diff --git a/test/unit/math/mix/prob/multinomial_test.cpp b/test/unit/math/mix/prob/multinomial_test.cpp
index 96f0bfbd59e..33317116e94 100644
--- a/test/unit/math/mix/prob/multinomial_test.cpp
+++ b/test/unit/math/mix/prob/multinomial_test.cpp
@@ -1,7 +1,6 @@
 #include <stan/math/mix.hpp>
 #include <test/unit/math/test_ad.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
diff --git a/test/unit/math/opencl/multiply_transpose_test.cpp b/test/unit/math/opencl/multiply_transpose_test.cpp
index 5f193aaacc1..69114beaf6d 100644
--- a/test/unit/math/opencl/multiply_transpose_test.cpp
+++ b/test/unit/math/opencl/multiply_transpose_test.cpp
@@ -4,10 +4,10 @@
 #include <stan/math/opencl/multiply_transpose.hpp>
 #include <stan/math/opencl/copy.hpp>
 #include <test/unit/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <algorithm>
-boost::random::mt19937 rng;
+std::mt19937 rng;
 
 TEST(MathMatrixOpenCL, multiply_transpose_exception_fail_zero) {
   stan::math::row_vector_d rv(0);
diff --git a/test/unit/math/opencl/prim/mdivide_left_tri_low_test.cpp b/test/unit/math/opencl/prim/mdivide_left_tri_low_test.cpp
index c0151631ed1..2120f9da5fc 100644
--- a/test/unit/math/opencl/prim/mdivide_left_tri_low_test.cpp
+++ b/test/unit/math/opencl/prim/mdivide_left_tri_low_test.cpp
@@ -3,7 +3,7 @@
 #include <stan/math/opencl/prim.hpp>
 #include <test/unit/math/expect_near_rel.hpp>
 #include <test/unit/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <algorithm>
 
@@ -29,7 +29,7 @@ TEST(MathMatrixCL, mdivide_left_tri_low_cl_exception) {
 }
 
 inline void mdivide_left_tri_low_Ab_test(int size) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   auto m1 = stan::math::matrix_d(size, size);
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < i; j++) {
@@ -57,7 +57,7 @@ inline void mdivide_left_tri_low_Ab_test(int size) {
 }
 
 inline void mdivide_left_tri_low_A_test(int size) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   auto m1 = stan::math::matrix_d(size, size);
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < i; j++) {
diff --git a/test/unit/math/opencl/prim/mdivide_right_tri_low_test.cpp b/test/unit/math/opencl/prim/mdivide_right_tri_low_test.cpp
index a339cfbb107..3674b7df2b8 100644
--- a/test/unit/math/opencl/prim/mdivide_right_tri_low_test.cpp
+++ b/test/unit/math/opencl/prim/mdivide_right_tri_low_test.cpp
@@ -3,7 +3,7 @@
 #include <stan/math/opencl/prim.hpp>
 #include <test/unit/math/expect_near_rel.hpp>
 #include <test/unit/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <algorithm>
 
@@ -29,7 +29,7 @@ TEST(MathMatrixCL, mdivide_right_tri_low_cl_exception) {
 }
 
 inline void mdivide_right_tri_low_Ab_test(int size) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   auto m1 = stan::math::matrix_d(size, size);
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < i; j++) {
diff --git a/test/unit/math/opencl/tri_inverse_test.cpp b/test/unit/math/opencl/tri_inverse_test.cpp
index 98c36794892..8ccce713a9d 100644
--- a/test/unit/math/opencl/tri_inverse_test.cpp
+++ b/test/unit/math/opencl/tri_inverse_test.cpp
@@ -2,7 +2,7 @@
 #include <stan/math/prim.hpp>
 #include <stan/math/opencl/prim.hpp>
 #include <test/unit/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <algorithm>
 
@@ -25,7 +25,7 @@ TEST(MathMatrixCL, inverse_cl_exception) {
 }
 
 inline void lower_inverse_test(int size) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   auto m1 = stan::math::matrix_d(size, size);
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < i; j++) {
@@ -56,7 +56,7 @@ inline void lower_inverse_test(int size) {
 }
 
 inline void upper_inverse_test(int size) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   auto m1 = stan::math::matrix_d(size, size);
   for (int i = 0; i < size; i++) {
     for (int j = 0; j < i; j++) {
diff --git a/test/unit/math/prim/fun/chol2inv_test.cpp b/test/unit/math/prim/fun/chol2inv_test.cpp
index 9942e366ac1..0db3b6a85ab 100644
--- a/test/unit/math/prim/fun/chol2inv_test.cpp
+++ b/test/unit/math/prim/fun/chol2inv_test.cpp
@@ -1,5 +1,5 @@
 #include <stan/math/prim.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 
 TEST(MathMatrixPrimMat, chol2inv_exception) {
@@ -25,7 +25,7 @@ TEST(MathMatrixPrimMat, chol2inv) {
   using stan::math::matrix_d;
   using stan::math::wishart_rng;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   matrix_d I(3, 3);
   I.setZero();
   I.diagonal().setOnes();
diff --git a/test/unit/math/prim/fun/promote_elements_test.cpp b/test/unit/math/prim/fun/promote_elements_test.cpp
index 20656134faf..c12ec59751b 100644
--- a/test/unit/math/prim/fun/promote_elements_test.cpp
+++ b/test/unit/math/prim/fun/promote_elements_test.cpp
@@ -1,5 +1,4 @@
 #include <stan/math/prim.hpp>
-#include <boost/typeof/typeof.hpp>
 #include <gtest/gtest.h>
 #include <type_traits>
 #include <vector>
@@ -10,7 +9,7 @@ TEST(MathFunctionsScalPromote_Elements, int2double) {
   using std::vector;
   int from;
   promote_elements<double, int> p;
-  typedef BOOST_TYPEOF(p.promote(from)) result_t;
+  typedef std::decay_t<decltype(p.promote(from))> result_t;
   bool same = std::is_same<double, result_t>::value;
   EXPECT_TRUE(same);
 }
@@ -21,7 +20,7 @@ TEST(MathFunctionsScalPromote_Elements, double2double) {
   using std::vector;
   double from;
   promote_elements<double, double> p;
-  typedef BOOST_TYPEOF(p.promote(from)) result_t;
+  typedef std::decay_t<decltype(p.promote(from))> result_t;
   bool same = std::is_same<double, result_t>::value;
   EXPECT_TRUE(same);
 }
@@ -35,7 +34,7 @@ TEST(MathFunctionsArrPromote_Elements, intVec2doubleVec) {
   from.push_back(2);
   from.push_back(3);
   promote_elements<vector<double>, vector<int> > p;
-  typedef BOOST_TYPEOF(p.promote(from)) result_t;
+  typedef std::decay_t<decltype(p.promote(from))> result_t;
   bool same = std::is_same<vector<double>, result_t>::value;
   EXPECT_TRUE(same);
 }
@@ -49,7 +48,7 @@ TEST(MathFunctionsArrPromote_Elements, doubleVec2doubleVec) {
   from.push_back(2);
   from.push_back(3);
   promote_elements<vector<double>, vector<double> > p;
-  typedef BOOST_TYPEOF(p.promote(from)) result_t;
+  typedef std::decay_t<decltype(p.promote(from))> result_t;
   bool same = std::is_same<vector<double>, result_t>::value;
   EXPECT_TRUE(same);
 }
@@ -61,7 +60,7 @@ TEST(MathFunctionsMatPromote_Elements, doubleMat2doubleMat) {
   stan::math::matrix_d m1(2, 3);
   m1 << 1, 2, 3, 4, 5, 6;
   promote_elements<Matrix<double, 2, 3>, Matrix<double, 2, 3> > p;
-  typedef BOOST_TYPEOF(p.promote(m1)) result_t;
+  typedef std::decay_t<decltype(p.promote(m1))> result_t;
   bool same = std::is_same<Matrix<double, 2, 3>, result_t>::value;
   EXPECT_TRUE(same);
 }
diff --git a/test/unit/math/prim/prob/bernoulli_logit_glm_rng_test.cpp b/test/unit/math/prim/prob/bernoulli_logit_glm_rng_test.cpp
index 35e800b3a63..87b32442ea0 100644
--- a/test/unit/math/prim/prob/bernoulli_logit_glm_rng_test.cpp
+++ b/test/unit/math/prim/prob/bernoulli_logit_glm_rng_test.cpp
@@ -1,7 +1,7 @@
 #include <gtest/gtest.h>
 #include <stan/math/prim.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <limits>
 #include <vector>
@@ -9,7 +9,7 @@
 TEST(ProbDistributionsBernoulliLogitGlm, vectorized) {
   using stan::math::bernoulli_logit_glm_rng;
   //  Test scalar/vector combinations.
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   Eigen::MatrixXd x(2, 3);
   x << 3.5, -1.5, 0.0, 2.0, 1.0, 3.0;
@@ -70,7 +70,7 @@ TEST(ProbDistributionsBernoulliLogitGlm, vectorized) {
 TEST(ProbDistributionsBernoulliLogitGlm, errorCheck) {
   using stan::math::bernoulli_logit_glm_rng;
   // Check errors for nonfinite and wrong sizes.
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   int N = 3;
   int M = 2;
@@ -107,7 +107,7 @@ TEST(ProbDistributionsBernoulliLogitGlm, errorCheck) {
 TEST(ProbDistributionsBernoulliLogitGlm, marginalChiSquareGoodnessFitTest) {
   using stan::math::bernoulli_logit_glm_rng;
   // Check distribution of result.
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Eigen::MatrixXd x(2, 2);
   x << 3.5, -1.5, 2.0, -1.2;
   std::vector<double> alpha{2.0, 1.0};
diff --git a/test/unit/math/prim/prob/bernoulli_logit_test.cpp b/test/unit/math/prim/prob/bernoulli_logit_test.cpp
index e608decd807..72ba4e87dec 100644
--- a/test/unit/math/prim/prob/bernoulli_logit_test.cpp
+++ b/test/unit/math/prim/prob/bernoulli_logit_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -35,7 +35,7 @@ TEST(ProbDistributionsBernoulliLogit, distributionCheck) {
 }
 
 TEST(ProbDistributionsBernoulliLogit, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   EXPECT_NO_THROW(stan::math::bernoulli_logit_rng(-3.5, rng));
   EXPECT_THROW(
@@ -44,7 +44,7 @@ TEST(ProbDistributionsBernoulliLogit, error_check) {
 }
 
 TEST(ProbDistributionsBernoulliLogit, logitChiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   // number of samples
   int N = 10000;
 
diff --git a/test/unit/math/prim/prob/bernoulli_test.cpp b/test/unit/math/prim/prob/bernoulli_test.cpp
index a1e3b722d7a..df78b190439 100644
--- a/test/unit/math/prim/prob/bernoulli_test.cpp
+++ b/test/unit/math/prim/prob/bernoulli_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -35,7 +35,7 @@ TEST(ProbDistributionsBernoulli, distributionCheck) {
 }
 
 TEST(ProbDistributionsBernoulli, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::bernoulli_rng(0.6, rng));
 
   EXPECT_THROW(stan::math::bernoulli_rng(1.6, rng), std::domain_error);
@@ -45,7 +45,7 @@ TEST(ProbDistributionsBernoulli, error_check) {
 }
 
 TEST(ProbDistributionsBernoulli, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
 
   std::vector<double> expected;
diff --git a/test/unit/math/prim/prob/beta_binomial_test.cpp b/test/unit/math/prim/prob/beta_binomial_test.cpp
index 9adc12fe794..e0c8968e9c3 100644
--- a/test/unit/math/prim/prob/beta_binomial_test.cpp
+++ b/test/unit/math/prim/prob/beta_binomial_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -40,7 +40,7 @@ TEST(ProbDistributionsBetaBinomial, distributionCheck) {
 }
 
 TEST(ProbDistributionBetaBinomial, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::beta_binomial_rng(4, 0.6, 2.0, rng));
 
   EXPECT_THROW(stan::math::beta_binomial_rng(-4, 0.6, 2, rng),
diff --git a/test/unit/math/prim/prob/beta_neg_binomial_test.cpp b/test/unit/math/prim/prob/beta_neg_binomial_test.cpp
index 1a4a848daf3..e19b8b87322 100644
--- a/test/unit/math/prim/prob/beta_neg_binomial_test.cpp
+++ b/test/unit/math/prim/prob/beta_neg_binomial_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -37,7 +37,7 @@ TEST(ProbDistributionsBetaNegBinomial, distributionCheck) {
 }
 
 TEST(ProbDistributionBetaNegBinomial, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::beta_neg_binomial_rng(4, 0.6, 2.0, rng));
 
   EXPECT_THROW(stan::math::beta_neg_binomial_rng(-4, 0.6, 2, rng),
diff --git a/test/unit/math/prim/prob/beta_proportion_test.cpp b/test/unit/math/prim/prob/beta_proportion_test.cpp
index f38f7e71e16..81fe6e865d2 100644
--- a/test/unit/math/prim/prob/beta_proportion_test.cpp
+++ b/test/unit/math/prim/prob/beta_proportion_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -50,7 +50,7 @@ TEST(ProbDistributionsBetaProportion, distributionTest) {
 }
 
 TEST(ProbDistributionsBetaProportion, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::beta_proportion_rng(0.5, 3.0, rng));
   EXPECT_NO_THROW(stan::math::beta_proportion_rng(1e-10, 1e-10, rng));
   EXPECT_THROW(stan::math::beta_proportion_rng(2.0, -1.0, rng),
@@ -68,7 +68,7 @@ TEST(ProbDistributionsBetaProportion, error_check) {
 }
 
 TEST(ProbDistributionsBetaProportion, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
@@ -99,7 +99,7 @@ TEST(ProbDistributionsBetaProportion, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsBetaProportion, chiSquareGoodnessFitTest2) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/beta_test.cpp b/test/unit/math/prim/prob/beta_test.cpp
index 2fa451a1288..9aa7efe7efd 100644
--- a/test/unit/math/prim/prob/beta_test.cpp
+++ b/test/unit/math/prim/prob/beta_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsBeta, distributionTest) {
 }
 
 TEST(ProbDistributionsBeta, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::beta_rng(2.0, 1.0, rng));
   EXPECT_NO_THROW(stan::math::beta_rng(1e-10, 1e-10, rng));
 
@@ -60,7 +60,7 @@ TEST(ProbDistributionsBeta, error_check) {
 }
 
 TEST(ProbDistributionsBeta, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
@@ -83,7 +83,7 @@ TEST(ProbDistributionsBeta, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsBeta, chiSquareGoodnessFitTestSmallParameters) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/binomial_test.cpp b/test/unit/math/prim/prob/binomial_test.cpp
index c12507db534..0d4d2ae24bb 100644
--- a/test/unit/math/prim/prob/binomial_test.cpp
+++ b/test/unit/math/prim/prob/binomial_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -39,7 +39,7 @@ TEST(ProbDistributionsBinomial, distributionCheck) {
 }
 
 TEST(ProbDistributionBinomiali, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::binomial_rng(4, 0.6, rng));
   EXPECT_THROW(stan::math::binomial_rng(-4, 0.6, rng), std::domain_error);
   EXPECT_THROW(stan::math::binomial_rng(4, -0.6, rng), std::domain_error);
@@ -50,7 +50,7 @@ TEST(ProbDistributionBinomiali, error_check) {
 }
 
 TEST(ProbDistributionsBinomial, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::binomial_distribution<> dist(100, 0.6);
diff --git a/test/unit/math/prim/prob/categorical_logit_rng_test.cpp b/test/unit/math/prim/prob/categorical_logit_rng_test.cpp
index b60603f28e7..50933863981 100644
--- a/test/unit/math/prim/prob/categorical_logit_rng_test.cpp
+++ b/test/unit/math/prim/prob/categorical_logit_rng_test.cpp
@@ -1,13 +1,13 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 
 TEST(ProbDistributionsCategoricalLogit, error_check) {
   using Eigen::VectorXd;
   using stan::math::categorical_logit_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   VectorXd beta(3);
 
@@ -27,7 +27,7 @@ TEST(ProbDistributionsCategoricalLogit, error_check) {
 TEST(ProbDistributionsCategoricalLogit, chiSquareGoodnessFitTest) {
   using Eigen::VectorXd;
   using stan::math::softmax;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = 3;
   VectorXd beta(K);
diff --git a/test/unit/math/prim/prob/categorical_test.cpp b/test/unit/math/prim/prob/categorical_test.cpp
index 6923d506870..0755b8a6902 100644
--- a/test/unit/math/prim/prob/categorical_test.cpp
+++ b/test/unit/math/prim/prob/categorical_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <vector>
 #include <limits>
@@ -88,7 +88,7 @@ TEST(ProbDistributionsCategorical, error_check) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using stan::math::categorical_lpmf;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   Matrix<double, Dynamic, Dynamic> theta(3, 1);
   theta << 0.15, 0.45, 0.50;
@@ -100,7 +100,7 @@ TEST(ProbDistributionsCategorical, chiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using stan::math::categorical_lpmf;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   int N = 10000;
   Matrix<double, Dynamic, Dynamic> theta(3, 1);
diff --git a/test/unit/math/prim/prob/cauchy_test.cpp b/test/unit/math/prim/prob/cauchy_test.cpp
index 7409f7bcc73..a99a5dc4131 100644
--- a/test/unit/math/prim/prob/cauchy_test.cpp
+++ b/test/unit/math/prim/prob/cauchy_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsCauchy, distributionTest) {
 }
 
 TEST(ProbDistributionsCauchy, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::cauchy_rng(2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::cauchy_rng(2.0, -1.0, rng), std::domain_error);
@@ -59,7 +59,7 @@ TEST(ProbDistributionsCauchy, error_check) {
 }
 
 TEST(ProbDistributionsCauchy, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/chi_square_test.cpp b/test/unit/math/prim/prob/chi_square_test.cpp
index 1e19844a18f..fbeae073aa1 100644
--- a/test/unit/math/prim/prob/chi_square_test.cpp
+++ b/test/unit/math/prim/prob/chi_square_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -42,7 +42,7 @@ TEST(ProbDistributionsChiSquare, distributionTest) {
 }
 
 TEST(ProbDistributionsChiSquare, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::chi_square_rng(2.0, rng));
 
   EXPECT_THROW(stan::math::chi_square_rng(-2.0, rng), std::domain_error);
@@ -51,7 +51,7 @@ TEST(ProbDistributionsChiSquare, error_check) {
 }
 
 TEST(ProbDistributionsChiSquare, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/dirichlet_multinomial_test.cpp b/test/unit/math/prim/prob/dirichlet_multinomial_test.cpp
index 74de88fe44b..61e1679ecf5 100644
--- a/test/unit/math/prim/prob/dirichlet_multinomial_test.cpp
+++ b/test/unit/math/prim/prob/dirichlet_multinomial_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -90,7 +90,7 @@ TEST(ProbDistributionsDirichletMultinomial, zeros) {
 TEST(ProbDistributionsDirichletMultinomial, rng) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   Matrix<double, Dynamic, 1> theta(3);
   theta << 0.5, 5.5, 12.0;
@@ -113,7 +113,7 @@ TEST(ProbDistributionsDirichletMultinomial, rng) {
 TEST(ProbDistributionsDirichletMultinomial, rngError) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   Matrix<double, Dynamic, 1> theta(3);
   theta << 0.5, 1.5, 4.0;
@@ -135,7 +135,7 @@ TEST(ProbDistributionsDirichletMultinomial, rngError) {
 TEST(ProbDistributionsDirichletMultinomial, chiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int M = 100;
   int trials = 100;
   int N = M * trials;
@@ -171,7 +171,7 @@ TEST(ProbDistributionsDirichletMultinomial, equivBetaBinomial) {
   using Eigen::Matrix;
   using stan::math::beta_binomial_lpmf;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   std::vector<int> ns;
   Matrix<double, Dynamic, 1> theta(2, 1);
   double alpha = 4.0;
diff --git a/test/unit/math/prim/prob/dirichlet_test.cpp b/test/unit/math/prim/prob/dirichlet_test.cpp
index e2386ed08a7..ee710a131db 100644
--- a/test/unit/math/prim/prob/dirichlet_test.cpp
+++ b/test/unit/math/prim/prob/dirichlet_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
@@ -139,7 +139,7 @@ inline double chi_square(std::vector<int> bin, std::vector<double> expect) {
 }
 
 inline void test_dirichlet3_1(Eigen::VectorXd alpha) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
@@ -168,7 +168,7 @@ inline void test_dirichlet3_2(Eigen::VectorXd alpha) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using Eigen::VectorXd;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::beta_distribution<> dist(alpha(1), alpha(0) + alpha(2));
@@ -213,7 +213,7 @@ TEST(ProbDistributionsDirichlet, random) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using Eigen::VectorXd;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   VectorXd alpha(3);
   alpha << 2.0, 3.0, 11.0;
   EXPECT_NO_THROW(stan::math::dirichlet_rng(alpha, rng));
diff --git a/test/unit/math/prim/prob/discrete_range_test.cpp b/test/unit/math/prim/prob/discrete_range_test.cpp
index d6535d384f3..2e4558d3ebf 100644
--- a/test/unit/math/prim/prob/discrete_range_test.cpp
+++ b/test/unit/math/prim/prob/discrete_range_test.cpp
@@ -1,12 +1,12 @@
 #include <stan/math/prim.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <vector>
 
 TEST(ProbDistributionsDiscreteRange, error_check) {
   using stan::math::discrete_range_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   std::vector<int> lower{5, 11, -15};
   std::vector<int> upper{7, 15, -10};
@@ -22,7 +22,7 @@ TEST(ProbDistributionsDiscreteRange, error_check) {
 
 TEST(ProbDistributionsDiscreteRange, boundary_values) {
   using stan::math::discrete_range_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   std::vector<int> lower{-5, 11, 17};
   EXPECT_EQ(lower, discrete_range_rng(lower, lower, rng));
@@ -37,7 +37,7 @@ TEST(ProbDistributionsDiscreteRange, boundary_values) {
 }
 
 TEST(ProbDistributionsDiscreteRange, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   int N = 10000;
   int lower = -3;
diff --git a/test/unit/math/prim/prob/double_exponential_test.cpp b/test/unit/math/prim/prob/double_exponential_test.cpp
index 8e8d677e226..ae0ab8675a6 100644
--- a/test/unit/math/prim/prob/double_exponential_test.cpp
+++ b/test/unit/math/prim/prob/double_exponential_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsDoubleExponential, distributionTest) {
 }
 
 TEST(ProbDistributionsDoubleExponential, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::double_exponential_rng(2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::double_exponential_rng(2.0, -1.0, rng),
@@ -60,7 +60,7 @@ TEST(ProbDistributionsDoubleExponential, error_check) {
 }
 
 TEST(ProbDistributionsDoubleExponential, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/exp_mod_normal_test.cpp b/test/unit/math/prim/prob/exp_mod_normal_test.cpp
index d817dd10342..33179176436 100644
--- a/test/unit/math/prim/prob/exp_mod_normal_test.cpp
+++ b/test/unit/math/prim/prob/exp_mod_normal_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <vector>
 
@@ -30,7 +30,7 @@ TEST(ProbDistributionsExpModNormal, errorCheck) {
  * the distributions manually
  */
 TEST(ProbDistributionsExpModNormal, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
 
   double mu = 2.0;
@@ -93,7 +93,7 @@ TEST(ProbDistributionsExpModNormal, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsExpModNormal, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::exp_mod_normal_rng(10.0, 2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::exp_mod_normal_rng(10.0, 2.0, -1.0, rng),
diff --git a/test/unit/math/prim/prob/exponential_test.cpp b/test/unit/math/prim/prob/exponential_test.cpp
index c581dd64508..615d9432fdd 100644
--- a/test/unit/math/prim/prob/exponential_test.cpp
+++ b/test/unit/math/prim/prob/exponential_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -43,7 +43,7 @@ TEST(ProbDistributionsExponential, distributionTest) {
 }
 
 TEST(ProbDistributionsExponential, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::exponential_rng(2.0, rng));
 
   EXPECT_THROW(stan::math::exponential_rng(-2.0, rng), std::domain_error);
@@ -53,7 +53,7 @@ TEST(ProbDistributionsExponential, error_check) {
 }
 
 TEST(ProbDistributionsExponential, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/frechet_test.cpp b/test/unit/math/prim/prob/frechet_test.cpp
index 8dd29a48bfb..796cc8f67df 100644
--- a/test/unit/math/prim/prob/frechet_test.cpp
+++ b/test/unit/math/prim/prob/frechet_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -51,7 +51,7 @@ TEST(ProbDistributionsFrechet, distributionTest) {
 }
 
 TEST(ProbDistributionsFrechet, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::frechet_rng(2.0, 3.0, rng));
 
   EXPECT_THROW(stan::math::frechet_rng(-2.0, 3.0, rng), std::domain_error);
@@ -62,7 +62,7 @@ TEST(ProbDistributionsFrechet, error_check) {
 }
 
 TEST(ProbDistributionsFrechet, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/gamma_test.cpp b/test/unit/math/prim/prob/gamma_test.cpp
index 63f7d396ee5..408dc28ae70 100644
--- a/test/unit/math/prim/prob/gamma_test.cpp
+++ b/test/unit/math/prim/prob/gamma_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsGamma, distributionTest) {
 }
 
 TEST(ProbDistributionGamma, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::gamma_rng(2.0, 3.0, rng));
 
   EXPECT_THROW(stan::math::gamma_rng(-2.0, 3.0, rng), std::domain_error);
@@ -58,7 +58,7 @@ TEST(ProbDistributionGamma, error_check) {
 }
 
 TEST(ProbDistributionGamma, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/gaussian_dlm_obs_rng_test.cpp b/test/unit/math/prim/prob/gaussian_dlm_obs_rng_test.cpp
index d49ae9a4d47..f0fc005ee64 100644
--- a/test/unit/math/prim/prob/gaussian_dlm_obs_rng_test.cpp
+++ b/test/unit/math/prim/prob/gaussian_dlm_obs_rng_test.cpp
@@ -1,7 +1,7 @@
 #include <gtest/gtest.h>
 #include <stan/math/prim.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <exception>
 #include <limits>
@@ -36,7 +36,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesF) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   MatrixXd FF_sz1 = MatrixXd::Random(4, 3);
   EXPECT_THROW(gaussian_dlm_obs_rng(FF_sz1, GG, V, W, m0, C0, T, rng),
                std::invalid_argument);
@@ -66,7 +66,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesG) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   // size
   MatrixXd GG_sz1 = MatrixXd::Random(3, 3);
   EXPECT_THROW(gaussian_dlm_obs_rng(FF, GG_sz1, V, W, m0, C0, T, rng),
@@ -97,7 +97,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesW) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   // Not symmetric
   MatrixXd W_asym = W;
   W_asym(0, 1) = 1;
@@ -155,7 +155,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesVMatrix) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   // Not symmetric
   MatrixXd V_asym = V;
   V_asym(0, 2) = 1;
@@ -198,7 +198,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesVVector) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   // negative
   MatrixXd V_neg = V_vec;
   V_neg(0) = -1;
@@ -230,7 +230,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, Policiesm0) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   // size
   VectorXd m0_sz = VectorXd::Zero(4, 1);
   EXPECT_THROW(gaussian_dlm_obs_rng(FF, GG, V, W, m0_sz, C0, T, rng),
@@ -256,7 +256,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesC0) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   // size
   MatrixXd C0_sz = MatrixXd::Identity(3, 3);
   EXPECT_THROW(gaussian_dlm_obs_rng(FF, GG, V, W, m0, C0_sz, T, rng),
@@ -293,7 +293,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, PoliciesT) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   // Must be positive.
   EXPECT_THROW(gaussian_dlm_obs_rng(FF, GG, V, W, m0, C0, 0, rng),
                std::domain_error);
@@ -305,7 +305,7 @@ TEST_F(ProbDistributionsGaussianDLMInputsRng, chiSquaredGoodnessOfFit) {
   using Eigen::MatrixXd;
   using Eigen::VectorXd;
   using stan::math::gaussian_dlm_obs_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   // With identity state transition and initial mean 0 and identity
   // covariance matrices:
diff --git a/test/unit/math/prim/prob/gumbel_test.cpp b/test/unit/math/prim/prob/gumbel_test.cpp
index 9de4ec2c765..aa81038a743 100644
--- a/test/unit/math/prim/prob/gumbel_test.cpp
+++ b/test/unit/math/prim/prob/gumbel_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsGumbel, distributionTest) {
 }
 
 TEST(ProbDistributionsGumbel, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::gumbel_rng(10.0, 2.0, rng));
 
   EXPECT_THROW(
@@ -56,7 +56,7 @@ TEST(ProbDistributionsGumbel, error_check) {
 }
 
 TEST(ProbDistributionsGumbel, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp b/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp
index 747784ab5b9..c15db51d8d3 100644
--- a/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp
+++ b/test/unit/math/prim/prob/hmm_hidden_state_prob_test.cpp
@@ -1,7 +1,6 @@
 #include <test/unit/math/prim/prob/hmm_util.hpp>
 #include <stan/math/prim/prob/hmm_hidden_state_prob.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random.hpp>
 #include <test/unit/math/test_ad.hpp>
 #include <test/unit/util.hpp>
 #include <gtest/gtest.h>
diff --git a/test/unit/math/prim/prob/hmm_latent_rng_test.cpp b/test/unit/math/prim/prob/hmm_latent_rng_test.cpp
index 0a5606e9abd..6d690f8a152 100644
--- a/test/unit/math/prim/prob/hmm_latent_rng_test.cpp
+++ b/test/unit/math/prim/prob/hmm_latent_rng_test.cpp
@@ -3,7 +3,7 @@
 #include <stan/math/prim/prob/hmm_latent_rng.hpp>
 #include <stan/math/prim/prob/chi_square_lcdf.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random.hpp>
+#include <random>
 #include <test/unit/math/test_ad.hpp>
 #include <test/unit/util.hpp>
 #include <gtest/gtest.h>
@@ -26,7 +26,7 @@ TEST(hmm_rng_test, chiSquareGoodnessFitTest) {
   Eigen::MatrixXd log_omegas
       = Eigen::MatrixXd::Ones(n_states, n_transitions + 1);
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
 
   std::vector<double> expected;
@@ -68,7 +68,7 @@ TEST(hmm_rng_test, chiSquareGoodnessFitTest_symmetric) {
   Eigen::MatrixXd log_omegas
       = Eigen::MatrixXd::Ones(n_states, n_transitions + 1);
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
 
   std::vector<double> expected_0;
diff --git a/test/unit/math/prim/prob/hmm_marginal_test.cpp b/test/unit/math/prim/prob/hmm_marginal_test.cpp
index 3f9b653eb31..80e18a7d3bd 100644
--- a/test/unit/math/prim/prob/hmm_marginal_test.cpp
+++ b/test/unit/math/prim/prob/hmm_marginal_test.cpp
@@ -1,7 +1,6 @@
 #include <test/unit/math/prim/prob/hmm_util.hpp>
 #include <stan/math/prim/prob/hmm_marginal.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random.hpp>
 #include <test/unit/math/test_ad.hpp>
 #include <test/unit/util.hpp>
 #include <gtest/gtest.h>
diff --git a/test/unit/math/prim/prob/hmm_util.hpp b/test/unit/math/prim/prob/hmm_util.hpp
index fe114281107..8d347887489 100644
--- a/test/unit/math/prim/prob/hmm_util.hpp
+++ b/test/unit/math/prim/prob/hmm_util.hpp
@@ -2,7 +2,6 @@
 #define TEST_UNIT_MATH_PRIM_PROB_HMM_UTIL
 #include <stan/math/prim/prob/hmm_marginal.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random.hpp>
 #include <test/unit/math/test_ad.hpp>
 #include <test/unit/util.hpp>
 #include <gtest/gtest.h>
diff --git a/test/unit/math/prim/prob/hypergeometric_test.cpp b/test/unit/math/prim/prob/hypergeometric_test.cpp
index d84fc39c781..748dfff8e9e 100644
--- a/test/unit/math/prim/prob/hypergeometric_test.cpp
+++ b/test/unit/math/prim/prob/hypergeometric_test.cpp
@@ -1,11 +1,11 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
 TEST(ProbDistributionsHypergeometric, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::hypergeometric_rng(10, 10, 15, rng));
 
   EXPECT_THROW(stan::math::hypergeometric_rng(30, 10, 15, rng),
@@ -19,7 +19,7 @@ TEST(ProbDistributionsHypergeometric, error_check) {
 }
 
 TEST(ProbDistributionsHypergeometric, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int num_draws = 10;
   int K = num_draws;
diff --git a/test/unit/math/prim/prob/inv_chi_square_test.cpp b/test/unit/math/prim/prob/inv_chi_square_test.cpp
index f9fd797e1c6..41f8b5b9e55 100644
--- a/test/unit/math/prim/prob/inv_chi_square_test.cpp
+++ b/test/unit/math/prim/prob/inv_chi_square_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -42,7 +42,7 @@ TEST(ProbDistributionsInvChiSquare, distributionTest) {
 }
 
 TEST(ProbDistributionsInvChiSquare, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::inv_chi_square_rng(4.0, rng));
 
   EXPECT_THROW(stan::math::inv_chi_square_rng(-4.0, rng), std::domain_error);
@@ -52,7 +52,7 @@ TEST(ProbDistributionsInvChiSquare, error_check) {
 }
 
 TEST(ProbDistributionsInvChiSquare, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/inv_gamma_test.cpp b/test/unit/math/prim/prob/inv_gamma_test.cpp
index 4687a03df55..becc8dab16d 100644
--- a/test/unit/math/prim/prob/inv_gamma_test.cpp
+++ b/test/unit/math/prim/prob/inv_gamma_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsInvGamma, distributionTest) {
 }
 
 TEST(ProbDistributionsInvGamma, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::inv_gamma_rng(4.0, 3.0, rng));
 
   EXPECT_THROW(stan::math::inv_gamma_rng(-4.0, 3.0, rng), std::domain_error);
@@ -60,7 +60,7 @@ TEST(ProbDistributionsInvGamma, error_check) {
 }
 
 TEST(ProbDistributionsInvGamma, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp
index 95369fffa58..3db158f8768 100644
--- a/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp
+++ b/test/unit/math/prim/prob/inv_wishart_cholesky_rng_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions/chi_squared.hpp>
 #include <gtest/gtest.h>
 #include <stdexcept>
@@ -13,7 +13,7 @@ TEST(ProbDistributionsInvWishartCholesky, rng) {
 
   using stan::math::inv_wishart_cholesky_rng;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   MatrixXd omega(3, 4);
   EXPECT_THROW(inv_wishart_cholesky_rng(3.0, omega, rng),
@@ -35,7 +35,7 @@ TEST(ProbDistributionsInvWishartCholesky, rng_pos_def) {
   using Eigen::MatrixXd;
   using stan::math::inv_wishart_cholesky_rng;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   MatrixXd Sigma(2, 2);
   MatrixXd Sigma_non_pos_def(2, 2);
@@ -58,7 +58,7 @@ TEST(ProbDistributionsInvWishartCholesky, marginalTwoChiSquareGoodnessFitTest) {
   using stan::math::inv_wishart_cholesky_rng;
   using std::log;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 2.0, -3.0, 4.0, 0.0, 2.0, 0.0, 3.0;
 
@@ -91,7 +91,7 @@ TEST(ProbDistributionsInvWishartCholesky, SpecialRNGTest) {
   using stan::math::inv_wishart_cholesky_rng;
   using stan::math::multiply_lower_tri_self_transpose;
 
-  boost::random::mt19937 rng(92343U);
+  std::mt19937 rng(92343U);
   int N = 1e5;
   double tol = 0.1;
   for (int k = 1; k < 5; k++) {
@@ -123,7 +123,7 @@ TEST(ProbDistributionsInvWishartCholesky, compareToInvWishart) {
   using stan::math::multiply_lower_tri_self_transpose;
   using stan::math::qr_thin_Q;
 
-  boost::random::mt19937 rng(92343U);
+  std::mt19937 rng(92343U);
   int N = 1e4;
   double tol = 0.05;
   for (int k = 1; k < 4; k++) {
diff --git a/test/unit/math/prim/prob/inv_wishart_cholesky_test.cpp b/test/unit/math/prim/prob/inv_wishart_cholesky_test.cpp
index 062de6d554c..6a3203316ef 100644
--- a/test/unit/math/prim/prob/inv_wishart_cholesky_test.cpp
+++ b/test/unit/math/prim/prob/inv_wishart_cholesky_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 
diff --git a/test/unit/math/prim/prob/inv_wishart_rng_test.cpp b/test/unit/math/prim/prob/inv_wishart_rng_test.cpp
index 77c953e771a..77d3a3acca8 100644
--- a/test/unit/math/prim/prob/inv_wishart_rng_test.cpp
+++ b/test/unit/math/prim/prob/inv_wishart_rng_test.cpp
@@ -1,14 +1,14 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/util.hpp>
 #include <boost/math/distributions/chi_squared.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <stdexcept>
 
 TEST(ProbDistributionsInvWishart, rng) {
   using Eigen::MatrixXd;
   using stan::math::inv_wishart_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   MatrixXd omega(3, 4);
   EXPECT_THROW(inv_wishart_rng(3.0, omega, rng), std::invalid_argument);
@@ -26,7 +26,7 @@ TEST(ProbDistributionsInvWishart, rng_pos_def) {
   using Eigen::MatrixXd;
   using stan::math::inv_wishart_rng;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   MatrixXd Sigma(2, 2);
   MatrixXd Sigma_non_pos_def(2, 2);
@@ -46,7 +46,7 @@ TEST(ProbDistributionsInvWishart, rng_symmetry) {
   using stan::test::unit::expect_symmetric;
   using stan::test::unit::spd_rng;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   for (int k = 1; k < 20; ++k)
     for (double nu = k - 0.5; nu < k + 20; ++nu)
       for (int n = 0; n < 10; ++n)
@@ -61,7 +61,7 @@ TEST(ProbDistributionsInvWishart, chiSquareGoodnessFitTest) {
   using stan::math::inv_wishart_rng;
   using std::log;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
   int N = 10000;
@@ -92,7 +92,7 @@ TEST(ProbDistributionsInvWishart, SpecialRNGTest) {
   using Eigen::MatrixXd;
   using stan::math::inv_wishart_rng;
 
-  boost::random::mt19937 rng(1234U);
+  std::mt19937 rng(1234U);
   int N = 1e5;
   double tol = 0.1;
   for (int k = 1; k < 5; k++) {
diff --git a/test/unit/math/prim/prob/inv_wishart_test.cpp b/test/unit/math/prim/prob/inv_wishart_test.cpp
index 58f00da356d..60e15d2f05d 100644
--- a/test/unit/math/prim/prob/inv_wishart_test.cpp
+++ b/test/unit/math/prim/prob/inv_wishart_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 
diff --git a/test/unit/math/prim/prob/lkj_corr_test.cpp b/test/unit/math/prim/prob/lkj_corr_test.cpp
index 220f19eb828..3c655c888ce 100644
--- a/test/unit/math/prim/prob/lkj_corr_test.cpp
+++ b/test/unit/math/prim/prob/lkj_corr_test.cpp
@@ -1,10 +1,10 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsLkjCorr, testIdentity) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setZero();
@@ -18,7 +18,7 @@ TEST(ProbDistributionsLkjCorr, testIdentity) {
 }
 
 TEST(ProbDistributionsLkjCorr, testHalf) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setConstant(0.5);
@@ -33,7 +33,7 @@ TEST(ProbDistributionsLkjCorr, testHalf) {
 }
 
 TEST(ProbDistributionsLkjCorr, Sigma) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setZero();
@@ -55,7 +55,7 @@ TEST(ProbDistributionsLkjCorr, Sigma) {
 }
 
 TEST(ProbDistributionsLKJCorr, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::lkj_corr_cholesky_rng(5, 1.0, rng));
   EXPECT_NO_THROW(stan::math::lkj_corr_rng(5, 1.0, rng));
 
@@ -65,7 +65,7 @@ TEST(ProbDistributionsLKJCorr, error_check) {
 }
 
 TEST(ProbDistributionsLKJCorr, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::beta_distribution<> dist(2.5, 2.5);
@@ -101,7 +101,7 @@ TEST(ProbDistributionsLKJCorr, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsLkjCorrCholesky, testIdentity) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setZero();
@@ -115,7 +115,7 @@ TEST(ProbDistributionsLkjCorrCholesky, testIdentity) {
 }
 
 TEST(ProbDistributionsLkjCorrCholesky, testHalf) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setConstant(0.5);
diff --git a/test/unit/math/prim/prob/lkj_cov_test.cpp b/test/unit/math/prim/prob/lkj_cov_test.cpp
index 5e10b5d9ce5..9c8c05910f4 100644
--- a/test/unit/math/prim/prob/lkj_cov_test.cpp
+++ b/test/unit/math/prim/prob/lkj_cov_test.cpp
@@ -1,10 +1,10 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbDistributionsLkjCov, testIdentity) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setZero();
@@ -26,7 +26,7 @@ TEST(ProbDistributionsLkjCov, testIdentity) {
 }
 
 TEST(ProbDistributionsLkjCov, testHalf) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setConstant(0.5);
@@ -49,7 +49,7 @@ TEST(ProbDistributionsLkjCov, testHalf) {
 }
 
 TEST(ProbDistributionsLkjCov, ErrorChecks) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   unsigned int K = 4;
   Eigen::MatrixXd Sigma(K, K);
   Sigma.setZero();
diff --git a/test/unit/math/prim/prob/logistic_test.cpp b/test/unit/math/prim/prob/logistic_test.cpp
index 4af72d4c321..0ddbd9acd86 100644
--- a/test/unit/math/prim/prob/logistic_test.cpp
+++ b/test/unit/math/prim/prob/logistic_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <limits>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsLogistic, distributionTest) {
 }
 
 TEST(ProbDistributionsLogistic, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::logistic_rng(4.0, 3.0, rng));
 
   EXPECT_THROW(stan::math::logistic_rng(4.0, -3.0, rng), std::domain_error);
@@ -59,7 +59,7 @@ TEST(ProbDistributionsLogistic, error_check) {
 }
 
 TEST(ProbDistributionsLogistic, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/loglogistic_test.cpp b/test/unit/math/prim/prob/loglogistic_test.cpp
index 38422afd7ee..95a7123e719 100644
--- a/test/unit/math/prim/prob/loglogistic_test.cpp
+++ b/test/unit/math/prim/prob/loglogistic_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <limits>
 #include <vector>
 
@@ -32,7 +32,7 @@ TEST(ProbDistributionsLoglogistic, errorCheck) {
 }
 
 TEST(ProbDistributionsLoglogistic, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::loglogistic_rng(10.0, 2.0, rng));
 
   EXPECT_THROW(stan::math::loglogistic_rng(2.0, -1.0, rng), std::domain_error);
@@ -57,7 +57,7 @@ TEST(ProbDistributionsLoglogistic, test_sampling_icdf) {
 }
 
 TEST(ProbDistributionsLoglogistic, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/lognormal_test.cpp b/test/unit/math/prim/prob/lognormal_test.cpp
index bfc16aaa881..7373fd38887 100644
--- a/test/unit/math/prim/prob/lognormal_test.cpp
+++ b/test/unit/math/prim/prob/lognormal_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsLogNormalMat, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsLogNormalPrim, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::lognormal_rng(2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::lognormal_rng(2.0, -1.0, rng), std::domain_error);
@@ -59,7 +59,7 @@ TEST(ProbDistributionsLogNormalPrim, error_check) {
 }
 
 TEST(ProbDistributionsLogNormalPrim, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/matrix_normal_prec_rng_test.cpp b/test/unit/math/prim/prob/matrix_normal_prec_rng_test.cpp
index 18fad92196e..25f5c8b79f2 100644
--- a/test/unit/math/prim/prob/matrix_normal_prec_rng_test.cpp
+++ b/test/unit/math/prim/prob/matrix_normal_prec_rng_test.cpp
@@ -1,6 +1,6 @@
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <exception>
@@ -10,7 +10,7 @@
 TEST(ProbDistributionsMatrixNormalPrecRng, ErrorSigma) {
   using Eigen::MatrixXd;
   using stan::math::matrix_normal_prec_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   double nan = std::numeric_limits<double>::quiet_NaN();
   double inf = std::numeric_limits<double>::infinity();
   double ninf = -std::numeric_limits<double>::infinity();
@@ -58,7 +58,7 @@ TEST(ProbDistributionsMatrixNormalPrecRng, ErrorSigma) {
 TEST(ProbDistributionsMatrixNormalPrecRng, ErrorD) {
   using Eigen::MatrixXd;
   using stan::math::matrix_normal_prec_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   double nan = std::numeric_limits<double>::quiet_NaN();
   double inf = std::numeric_limits<double>::infinity();
   double ninf = -std::numeric_limits<double>::infinity();
@@ -106,7 +106,7 @@ TEST(ProbDistributionsMatrixNormalPrecRng, ErrorD) {
 TEST(ProbDistributionsMatrixNormalPrecRng, ErrorSize) {
   using Eigen::MatrixXd;
   using stan::math::matrix_normal_prec_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   MatrixXd Mu = MatrixXd::Zero(3, 5);
 
@@ -171,7 +171,7 @@ inline std::vector<double> extract_sum_of_entries(
 TEST(ProbDistributionsMatrixNormalPrecRng, marginalChiSquareGoodnessFitTest) {
   using Eigen::MatrixXd;
   using stan::math::matrix_normal_prec_rng;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   MatrixXd Mu(2, 3);
   Mu << 1, 2, 3, 4, 5, 6;
diff --git a/test/unit/math/prim/prob/multi_normal_cholesky_test.cpp b/test/unit/math/prim/prob/multi_normal_cholesky_test.cpp
index 03e5c9e089a..01e6e6cb5ee 100644
--- a/test/unit/math/prim/prob/multi_normal_cholesky_test.cpp
+++ b/test/unit/math/prim/prob/multi_normal_cholesky_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <vector>
 
@@ -8,7 +8,7 @@ TEST(ProbDistributionsMultiNormalCholesky, NotVectorized) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -23,7 +23,7 @@ TEST(ProbDistributionsMultiNormalCholesky, Vectorized) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   vector<Matrix<double, Dynamic, 1> > vec_y(2);
   vector<Matrix<double, 1, Dynamic> > vec_y_t(2);
   Matrix<double, Dynamic, 1> y(3);
@@ -89,7 +89,7 @@ TEST(ProbDistributionsMultiNormalCholesky, MultiNormalOneRow) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, 1, Dynamic> y(3);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3);
@@ -105,7 +105,7 @@ TEST(ProbDistributionsMultiNormalCholesky, error_check) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> mu(3, 1);
   mu << 2.0, -2.0, 11.0;
 
@@ -124,7 +124,7 @@ TEST(ProbDistributionsMultiNormalCholesky,
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
   Matrix<double, Dynamic, Dynamic> L = sigma.llt().matrixL();
@@ -173,7 +173,7 @@ TEST(ProbDistributionsMultiNormalCholesky,
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
   Matrix<double, Dynamic, Dynamic> L = sigma.llt().matrixL();
@@ -222,7 +222,7 @@ TEST(ProbDistributionsMultiNormalCholesky,
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 16.0;
   Matrix<double, Dynamic, Dynamic> L = sigma.llt().matrixL();
diff --git a/test/unit/math/prim/prob/multi_normal_prec_rng_test.cpp b/test/unit/math/prim/prob/multi_normal_prec_rng_test.cpp
index 33830607980..475d345df5d 100644
--- a/test/unit/math/prim/prob/multi_normal_prec_rng_test.cpp
+++ b/test/unit/math/prim/prob/multi_normal_prec_rng_test.cpp
@@ -1,6 +1,6 @@
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <limits>
@@ -8,7 +8,7 @@
 
 TEST(ProbDistributionsMultiNormalPrec, vectorized) {
   // Test scalar/vector combinations.
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   Eigen::VectorXd mu(3);
   Eigen::RowVectorXd mu_t(3);
@@ -33,7 +33,7 @@ TEST(ProbDistributionsMultiNormalPrec, vectorized) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, policiesSigma) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   double nan = std::numeric_limits<double>::quiet_NaN();
   double inf = std::numeric_limits<double>::infinity();
@@ -80,7 +80,7 @@ TEST(ProbDistributionsMultiNormalPrec, policiesSigma) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, policiesMu) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   double nan = std::numeric_limits<double>::quiet_NaN();
   double inf = std::numeric_limits<double>::infinity();
@@ -109,7 +109,7 @@ TEST(ProbDistributionsMultiNormalPrec, policiesMu) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, SizeMismatch) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Eigen::VectorXd mu(2);
   mu << 1.0, -1.0;
   Eigen::MatrixXd Sigma(3, 3);
@@ -119,7 +119,7 @@ TEST(ProbDistributionsMultiNormalPrec, SizeMismatch) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, marginalOneChiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   Eigen::MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
@@ -153,7 +153,7 @@ TEST(ProbDistributionsMultiNormalPrec, marginalOneChiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, marginalTwoChiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   Eigen::MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
@@ -187,7 +187,7 @@ TEST(ProbDistributionsMultiNormalPrec, marginalTwoChiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsMultiNormalPrec, marginalThreeChiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   Eigen::MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 16.0;
diff --git a/test/unit/math/prim/prob/multi_normal_test.cpp b/test/unit/math/prim/prob/multi_normal_test.cpp
index c8a56e62aa0..e0d5007d481 100644
--- a/test/unit/math/prim/prob/multi_normal_test.cpp
+++ b/test/unit/math/prim/prob/multi_normal_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -10,7 +10,7 @@ TEST(ProbDistributionsMultiNormal, NotVectorized) {
   using Eigen::Matrix;
   using std::vector;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -25,7 +25,7 @@ TEST(ProbDistributionsMultiNormal, Vectorized) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   vector<Matrix<double, Dynamic, 1> > vec_y(2);
   vector<Matrix<double, 1, Dynamic> > vec_y_t(2);
   Matrix<double, Dynamic, 1> y(3);
@@ -89,7 +89,7 @@ TEST(ProbDistributionsMultiNormal, Sigma) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(2, 1);
   y << 2.0, -2.0;
   Matrix<double, Dynamic, 1> mu(2, 1);
@@ -108,7 +108,7 @@ TEST(ProbDistributionsMultiNormal, Mu) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -132,7 +132,7 @@ TEST(ProbDistributionsMultiNormal, MultiNormalOneRow) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, 1, Dynamic> y(1, 3);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -147,7 +147,7 @@ TEST(ProbDistributionsMultiNormal, SigmaMultiRow) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, 1, Dynamic> y(1, 2);
   y << 2.0, -2.0;
   Matrix<double, Dynamic, 1> mu(2, 1);
@@ -172,7 +172,7 @@ TEST(ProbDistributionsMultiNormal, MuMultiRow) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, 1, Dynamic> y(1, 3);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -196,7 +196,7 @@ TEST(ProbDistributionsMultiNormal, SizeMismatch) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, 1, Dynamic> y(1, 3);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(2, 1);
@@ -213,7 +213,7 @@ TEST(ProbDistributionsMultiNormal, marginalOneChiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
   std::vector<Matrix<double, Dynamic, 1> > mu(3);
@@ -260,7 +260,7 @@ TEST(ProbDistributionsMultiNormal, marginalTwoChiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 3.0;
   std::vector<Matrix<double, 1, Dynamic> > mu(3);
@@ -307,7 +307,7 @@ TEST(ProbDistributionsMultiNormal, marginalThreeChiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, Dynamic> sigma(3, 3);
   sigma << 9.0, -3.0, 0.0, -3.0, 4.0, 1.0, 0.0, 1.0, 16.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -354,7 +354,7 @@ TEST(multiNormalRng, nonPosDefErrorTest) {
   S << 0, 1, 1, 0;  // not pos definite
   Eigen::VectorXd mu(2);
   mu << 1, 2;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_THROW(multi_normal_rng(mu, S, rng), std::domain_error);
 }
 
diff --git a/test/unit/math/prim/prob/multi_student_t_cholesky_test.cpp b/test/unit/math/prim/prob/multi_student_t_cholesky_test.cpp
index 5111af3312b..8febe06003a 100644
--- a/test/unit/math/prim/prob/multi_student_t_cholesky_test.cpp
+++ b/test/unit/math/prim/prob/multi_student_t_cholesky_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <vector>
 #include <limits>
@@ -11,7 +11,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, NotVectorized) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -34,7 +34,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, Vectorized) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   vector<Matrix<double, Dynamic, 1> > vec_y(2);
   vector<Matrix<double, 1, Dynamic> > vec_y_t(2);
   Matrix<double, Dynamic, 1> y(3);
@@ -115,7 +115,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, Sigma) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -139,7 +139,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, Mu) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -195,7 +195,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, Nu) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -253,7 +253,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, ErrorSize0) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y_empty(0, 1);
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
@@ -293,7 +293,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, ErrorSize2) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -315,7 +315,7 @@ TEST(ProbDistributionsMultiStudentTCholesky, ErrorSize3) {
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -356,7 +356,7 @@ TEST(ProbDistributionsMultiStudentTCholesky,
   using Eigen::Matrix;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> mu(3, 1);
   mu << 2.0, 3.0, 11.0;
 
@@ -401,7 +401,7 @@ TEST(ProbDistributionsMultiStudentTCholesky,
   using stan::math::multi_student_t_cholesky_lpdf;
   using stan::math::multi_student_t_cholesky_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> mu(3, 1);
   mu << 2.0, 3.0, 11.0;
 
diff --git a/test/unit/math/prim/prob/multi_student_t_test.cpp b/test/unit/math/prim/prob/multi_student_t_test.cpp
index 5f656450225..72273942bb3 100644
--- a/test/unit/math/prim/prob/multi_student_t_test.cpp
+++ b/test/unit/math/prim/prob/multi_student_t_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <vector>
 #include <limits>
@@ -11,7 +11,7 @@ TEST(ProbDistributionsMultiStudentT, NotVectorized) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -30,7 +30,7 @@ TEST(ProbDistributionsMultiStudentT, Vectorized) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   vector<Matrix<double, Dynamic, 1> > vec_y(2);
   vector<Matrix<double, 1, Dynamic> > vec_y_t(2);
   Matrix<double, Dynamic, 1> y(3);
@@ -100,7 +100,7 @@ TEST(ProbDistributionsMultiStudentT, Sigma) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -123,7 +123,7 @@ TEST(ProbDistributionsMultiStudentT, Mu) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -178,7 +178,7 @@ TEST(ProbDistributionsMultiStudentT, Nu) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -234,7 +234,7 @@ TEST(ProbDistributionsMultiStudentT, ErrorSize2) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -253,7 +253,7 @@ TEST(ProbDistributionsMultiStudentT, ErrorSize3) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -272,7 +272,7 @@ TEST(ProbDistributionsMultiStudentT, ErrorSizeSigma) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> y(3, 1);
   y << 2.0, -2.0, 11.0;
   Matrix<double, Dynamic, 1> mu(3, 1);
@@ -308,7 +308,7 @@ TEST(ProbDistributionsMultiStudentT, marginalOneChiSquareGoodnessFitTest) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> mu(3, 1);
   mu << 2.0, 3.0, 11.0;
 
@@ -351,7 +351,7 @@ TEST(ProbDistributionsMultiStudentT, marginalTwoChiSquareGoodnessFitTest) {
   using stan::math::multi_student_t_lpdf;
   using stan::math::multi_student_t_rng;
   using std::vector;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> mu(3, 1);
   mu << 2.0, 3.0, 11.0;
 
diff --git a/test/unit/math/prim/prob/multinomial_logit_test.cpp b/test/unit/math/prim/prob/multinomial_logit_test.cpp
index dd4aa7cff2b..9f0e49db98e 100644
--- a/test/unit/math/prim/prob/multinomial_logit_test.cpp
+++ b/test/unit/math/prim/prob/multinomial_logit_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -9,7 +9,7 @@ using Eigen::Dynamic;
 using Eigen::Matrix;
 
 TEST(ProbDistributionsMultinomialLogit, RNGZero) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> beta(3);
   beta << 1.3, 0.1, -2.6;
   // bug in 4.8.1: RNG does not allow a zero total count
@@ -22,7 +22,7 @@ TEST(ProbDistributionsMultinomialLogit, RNGZero) {
 }
 
 TEST(ProbDistributionsMultinomialLogit, RNGSize) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> beta(5);
   beta << log(0.3), log(0.1), log(0.2), log(0.2), log(0.2);
   std::vector<int> sample = stan::math::multinomial_logit_rng(beta, 10, rng);
@@ -106,7 +106,7 @@ TEST(ProbDistributionsMultinomialLogit, zeros) {
 }
 
 TEST(ProbDistributionsMultinomialLogit, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int M = 10;
   int trials = 1000;
   int N = M * trials;
diff --git a/test/unit/math/prim/prob/multinomial_test.cpp b/test/unit/math/prim/prob/multinomial_test.cpp
index 1ce01a4b87b..cd0c1f92e17 100644
--- a/test/unit/math/prim/prob/multinomial_test.cpp
+++ b/test/unit/math/prim/prob/multinomial_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -8,7 +8,7 @@
 TEST(ProbDistributionsMultinomial, RNGZero) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> theta(3);
   theta << 0.3, 0.1, 0.6;
   // bug in 4.8.1: RNG does not allow a zero total count
@@ -23,7 +23,7 @@ TEST(ProbDistributionsMultinomial, RNGZero) {
 TEST(ProbDistributionsMultinomial, RNGSize) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   Matrix<double, Dynamic, 1> theta(5);
   // error in 2.1.0 due to overflow in binomial call due to division
   theta << 0.3, 0.1, 0.2, 0.2, 0.2;
@@ -124,7 +124,7 @@ TEST(ProbDistributionsMultinomial, zeros) {
 TEST(ProbDistributionsMultinomial, error_check) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   Matrix<double, Dynamic, 1> theta(3);
   theta << 0.15, 0.45, 0.40;
@@ -138,7 +138,7 @@ TEST(ProbDistributionsMultinomial, error_check) {
 TEST(ProbDistributionsMultinomial, chiSquareGoodnessFitTest) {
   using Eigen::Dynamic;
   using Eigen::Matrix;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int M = 10;
   int trials = 1000;
   int N = M * trials;
diff --git a/test/unit/math/prim/prob/neg_binomial_2_ccdf_log_test.cpp b/test/unit/math/prim/prob/neg_binomial_2_ccdf_log_test.cpp
index 309a32c1bdc..6fa49af0058 100644
--- a/test/unit/math/prim/prob/neg_binomial_2_ccdf_log_test.cpp
+++ b/test/unit/math/prim/prob/neg_binomial_2_ccdf_log_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 
 TEST(ProbNegBinomial2, ccdf_log_matches_lccdf) {
diff --git a/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp b/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp
index b1ae1b2d5d1..f7e737aeb48 100644
--- a/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp
+++ b/test/unit/math/prim/prob/neg_binomial_2_log_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/NegativeBinomial2LogTestRig.hpp>
 #include <test/unit/math/expect_near_rel.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <algorithm>
 #include <limits>
@@ -17,7 +17,7 @@ TEST(ProbDistributionsNegativeBinomial2Log, errorCheck) {
 TEST(ProbDistributionsNegBinomial2Log, error_check) {
   using std::log;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::neg_binomial_2_log_rng(6, 2, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_2_log_rng(-0.5, 1, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_2_log_rng(log(1e8), 1, rng));
@@ -71,7 +71,7 @@ TEST(ProbDistributionsNegBinomial2Log, error_check) {
 }
 
 TEST(ProbDistributionsNegBinomial2Log, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(1.1,
@@ -109,7 +109,7 @@ TEST(ProbDistributionsNegBinomial2Log, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsNegBinomial2Log, chiSquareGoodnessFitTest2) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(0.6,
@@ -147,7 +147,7 @@ TEST(ProbDistributionsNegBinomial2Log, chiSquareGoodnessFitTest2) {
 }
 
 TEST(ProbDistributionsNegBinomial2Log, chiSquareGoodnessFitTest3) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(121,
diff --git a/test/unit/math/prim/prob/neg_binomial_2_test.cpp b/test/unit/math/prim/prob/neg_binomial_2_test.cpp
index 232e8b45289..33b96ed7d32 100644
--- a/test/unit/math/prim/prob/neg_binomial_2_test.cpp
+++ b/test/unit/math/prim/prob/neg_binomial_2_test.cpp
@@ -4,7 +4,7 @@
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <test/unit/math/expect_near_rel.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -39,7 +39,7 @@ TEST(ProbDistributionsNegBinomial2, distributionCheck) {
 }
 
 TEST(ProbDistributionsNegBinomial2, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::neg_binomial_2_rng(6, 2, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_2_rng(0.5, 1, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_2_rng(1e8, 1, rng));
@@ -91,7 +91,7 @@ TEST(ProbDistributionsNegBinomial2, error_check) {
 }
 
 TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(1.1, 1.1 / (1.1 + 2.4));
@@ -128,7 +128,7 @@ TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest2) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(0.6, 0.6 / (0.6 + 2.4));
@@ -165,7 +165,7 @@ TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest2) {
 }
 
 TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest3) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(30, 30 / (30 + 60.4));
@@ -202,7 +202,7 @@ TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest3) {
 }
 
 TEST(ProbDistributionsNegBinomial2, chiSquareGoodnessFitTest4) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::negative_binomial_distribution<> dist(80, 80 / (80 + 30.4));
diff --git a/test/unit/math/prim/prob/neg_binomial_test.cpp b/test/unit/math/prim/prob/neg_binomial_test.cpp
index 7d1240ccd85..ed5e316a537 100644
--- a/test/unit/math/prim/prob/neg_binomial_test.cpp
+++ b/test/unit/math/prim/prob/neg_binomial_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -37,7 +37,7 @@ TEST(ProbDistributionsNegativeBinomial, distributionCheck) {
 }
 
 TEST(ProbDistributionsNegBinomial, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::neg_binomial_rng(6, 2, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_rng(0.5, 1, rng));
   EXPECT_NO_THROW(stan::math::neg_binomial_rng(1e9, 1, rng));
@@ -82,7 +82,7 @@ inline void expected_bin_sizes(double* expect, const int K, const int N,
 }
 
 TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   double p = 0.6;
   double alpha = 5;
   double beta = p / (1 - p);
@@ -120,7 +120,7 @@ TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest2) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   double p = 0.8;
   double alpha = 2.4;
   double beta = p / (1 - p);
@@ -158,13 +158,12 @@ TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest2) {
 }
 
 TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest3) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   double p = 0.2;
   double alpha = 0.4;
   double beta = p / (1 - p);
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, (1 - p)));
-  boost::math::chi_squared mydist(K - 1);
 
   int loc[K - 1];
   for (int i = 1; i < K; i++)
@@ -178,6 +177,30 @@ TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest3) {
     bin[i] = 0;
   expected_bin_sizes(expect, K, N, alpha, beta);
 
+  // alpha < 1 makes this distribution's tail very sparse, so most of the K
+  // bins have expected counts far below N; a single sample landing in one
+  // of them can dominate the statistic. Pool bins so each group has an
+  // expected count of at least 5, per standard chi-square goodness-of-fit
+  // practice. Grouping is based only on the theoretical expected counts
+  // (not the realized samples), so it doesn't bias the test.
+  const double min_expected_count = 5.0;
+  std::vector<int> group_end;  // exclusive upper bound of each group
+  double acc_exp = 0;
+  for (int j = 0; j < K; j++) {
+    acc_exp += expect[j];
+    if (acc_exp >= min_expected_count) {
+      group_end.push_back(j + 1);
+      acc_exp = 0;
+    }
+  }
+  if (acc_exp > 0) {
+    if (!group_end.empty()) {
+      group_end.back() = K;
+    } else {
+      group_end.push_back(K);
+    }
+  }
+
   while (count < N) {
     int a = stan::math::neg_binomial_rng(alpha, beta, rng);
     int i = 0;
@@ -187,12 +210,20 @@ TEST(ProbDistributionsNegBinomial, chiSquareGoodnessFitTest3) {
     count++;
   }
 
-  double chi = 0;
+  int num_groups = group_end.size();
+  boost::math::chi_squared mydist(num_groups - 1);
 
-  for (int j = 0; j < K; j++) {
-    if (expect[j] != 0) {
-      chi += ((bin[j] - expect[j]) * (bin[j] - expect[j]) / expect[j]);
+  double chi = 0;
+  int start = 0;
+  for (int g = 0; g < num_groups; g++) {
+    double obs = 0;
+    double exp = 0;
+    for (int j = start; j < group_end[g]; j++) {
+      obs += bin[j];
+      exp += expect[j];
     }
+    chi += (obs - exp) * (obs - exp) / exp;
+    start = group_end[g];
   }
 
   EXPECT_LT(chi, boost::math::quantile(boost::math::complement(mydist, 1e-6)));
diff --git a/test/unit/math/prim/prob/normal_test.cpp b/test/unit/math/prim/prob/normal_test.cpp
index 1e4a2860ede..a2462211330 100644
--- a/test/unit/math/prim/prob/normal_test.cpp
+++ b/test/unit/math/prim/prob/normal_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -77,7 +77,7 @@ TEST(ProbDistributionsNormal, distributionTest) {
 }
 
 TEST(ProbDistributionsNormal, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::normal_rng(10.0, 2.0, rng));
 
   EXPECT_THROW(stan::math::normal_rng(10.0, -2.0, rng), std::domain_error);
@@ -90,7 +90,7 @@ TEST(ProbDistributionsNormal, error_check) {
 }
 
 TEST(ProbDistributionsNormal, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/ordered_logistic_test.cpp b/test/unit/math/prim/prob/ordered_logistic_test.cpp
index 423d211ea6a..9caafa1ed46 100644
--- a/test/unit/math/prim/prob/ordered_logistic_test.cpp
+++ b/test/unit/math/prim/prob/ordered_logistic_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -169,7 +169,7 @@ TEST(ProbDistributions, ordered_logistic) {
 }
 
 TEST(ProbDistributionOrderedLogistic, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   double inf = std::numeric_limits<double>::infinity();
   Eigen::VectorXd c(4);
   c << -2, 2.0, 5, 10;
@@ -189,7 +189,7 @@ TEST(ProbDistributionOrderedLogistic, error_check) {
 
 TEST(ProbDistributionOrderedLogistic, chiSquareGoodnessFitTest) {
   using stan::math::inv_logit;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   double eta = 1.0;
   Eigen::VectorXd theta(3);
diff --git a/test/unit/math/prim/prob/ordered_probit_test.cpp b/test/unit/math/prim/prob/ordered_probit_test.cpp
index c9c5094afcc..55018cfe46c 100644
--- a/test/unit/math/prim/prob/ordered_probit_test.cpp
+++ b/test/unit/math/prim/prob/ordered_probit_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 
@@ -140,7 +140,7 @@ TEST(ProbDistributions, ordered_probit) {
 }
 
 TEST(ProbDistributionOrderedProbit, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   double inf = std::numeric_limits<double>::infinity();
   Eigen::VectorXd c(4);
   c << -2, 2.0, 5, 10;
@@ -155,7 +155,7 @@ TEST(ProbDistributionOrderedProbit, error_check) {
 
 TEST(ProbDistributionOrderedProbit, chiSquareGoodnessFitTest) {
   using stan::math::Phi;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   double eta = 1.0;
   Eigen::VectorXd theta(3);
diff --git a/test/unit/math/prim/prob/pareto_test.cpp b/test/unit/math/prim/prob/pareto_test.cpp
index df1c57c64ca..1f726381fb9 100644
--- a/test/unit/math/prim/prob/pareto_test.cpp
+++ b/test/unit/math/prim/prob/pareto_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsPareto, distributionTest) {
 }
 
 TEST(ProbDistributionsPareto, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::pareto_rng(2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::pareto_rng(2.0, -1.0, rng), std::domain_error);
@@ -59,7 +59,7 @@ TEST(ProbDistributionsPareto, error_check) {
 }
 
 TEST(ProbDistributionsPareto, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/pareto_type_2_test.cpp b/test/unit/math/prim/prob/pareto_type_2_test.cpp
index 5af75607cf6..5e4100551fa 100644
--- a/test/unit/math/prim/prob/pareto_type_2_test.cpp
+++ b/test/unit/math/prim/prob/pareto_type_2_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 #include <vector>
@@ -31,7 +31,7 @@ TEST(ProbDistributionsParetoType2, errorCheck) {
  * the distributions manually
  */
 TEST(ProbDistributionsParetoType2Mat, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
 
   double mu = 3.0;
@@ -108,7 +108,7 @@ TEST(ProbDistributionsParetoType2Mat, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsParetoType2Prim, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::chi_squared mydist(K - 1);
diff --git a/test/unit/math/prim/prob/poisson_binomial_test.cpp b/test/unit/math/prim/prob/poisson_binomial_test.cpp
index 72726bd9f3c..2dcb6269427 100644
--- a/test/unit/math/prim/prob/poisson_binomial_test.cpp
+++ b/test/unit/math/prim/prob/poisson_binomial_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -215,7 +215,7 @@ TEST(ProbDistributionsPoissonBinomial,
 TEST(ProbDistributionsPoissonBinomial, chiSquareGoodnessFitTest) {
   using vec = Eigen::Matrix<double, Eigen::Dynamic, 1>;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::binomial_distribution<> dist(100, 0.6);
diff --git a/test/unit/math/prim/prob/poisson_log_test.cpp b/test/unit/math/prim/prob/poisson_log_test.cpp
index 0399d2bb712..dbf9f566448 100644
--- a/test/unit/math/prim/prob/poisson_log_test.cpp
+++ b/test/unit/math/prim/prob/poisson_log_test.cpp
@@ -3,7 +3,6 @@
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
 #include <limits>
 #include <vector>
 
diff --git a/test/unit/math/prim/prob/poisson_test.cpp b/test/unit/math/prim/prob/poisson_test.cpp
index 709395c481a..9ee2dcfe719 100644
--- a/test/unit/math/prim/prob/poisson_test.cpp
+++ b/test/unit/math/prim/prob/poisson_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -36,7 +36,7 @@ TEST(ProbDistributionsPoisson, distributionCheck) {
 TEST(ProbDistributionsPoisson, error_check) {
   using std::log;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::poisson_rng(6, rng));
 
   EXPECT_THROW(stan::math::poisson_rng(-6, rng), std::domain_error);
@@ -56,7 +56,7 @@ TEST(ProbDistributionsPoisson, error_check) {
 }
 
 TEST(ProbDistributionsPoisson, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::poisson_distribution<> dist(5);
@@ -95,7 +95,7 @@ TEST(ProbDistributionsPoisson, chiSquareGoodnessFitTest) {
 TEST(ProbDistributionsPoisson, chiSquareGoodnessFitTest2) {
   using std::log;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 1000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
   boost::math::poisson_distribution<> dist(5);
diff --git a/test/unit/math/prim/prob/rayleigh_test.cpp b/test/unit/math/prim/prob/rayleigh_test.cpp
index 4ca66c9825f..43c076d2dfa 100644
--- a/test/unit/math/prim/prob/rayleigh_test.cpp
+++ b/test/unit/math/prim/prob/rayleigh_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <limits>
 #include <vector>
 
@@ -43,14 +43,14 @@ TEST(ProbDistributionsRayleigh, distributionTest) {
 }
 
 TEST(ProbDistributionsRayleigh, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::rayleigh_rng(2.0, rng));
 
   EXPECT_THROW(stan::math::rayleigh_rng(-2.0, rng), std::domain_error);
 }
 
 TEST(ProbDistributionsRayleigh, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/scaled_inv_chi_square_test.cpp b/test/unit/math/prim/prob/scaled_inv_chi_square_test.cpp
index 0999d35b88c..c2509067f72 100644
--- a/test/unit/math/prim/prob/scaled_inv_chi_square_test.cpp
+++ b/test/unit/math/prim/prob/scaled_inv_chi_square_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <limits>
 #include <vector>
 
@@ -46,7 +46,7 @@ TEST(ProbDistributionsScaledInvChiSquare, distributionTest) {
 }
 
 TEST(ProbDistributionsScaledInvChiSquare, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::scaled_inv_chi_square_rng(2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::scaled_inv_chi_square_rng(-2.0, 1.0, rng),
@@ -62,7 +62,7 @@ TEST(ProbDistributionsScaledInvChiSquare, error_check) {
 }
 
 TEST(ProbDistributionsScaledInvChiSquare, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/skew_double_exponential_ccdf_log_test.cpp b/test/unit/math/prim/prob/skew_double_exponential_ccdf_log_test.cpp
index f1769ee046c..9c1974dad84 100644
--- a/test/unit/math/prim/prob/skew_double_exponential_ccdf_log_test.cpp
+++ b/test/unit/math/prim/prob/skew_double_exponential_ccdf_log_test.cpp
@@ -3,7 +3,6 @@
 #include <stan/math/prim/fun/fabs.hpp>
 #include <stan/math/prim/fun/log1m.hpp>
 
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/tools/promotion.hpp>
 #include <gtest/gtest.h>
diff --git a/test/unit/math/prim/prob/skew_double_exponential_cdf_log_test.cpp b/test/unit/math/prim/prob/skew_double_exponential_cdf_log_test.cpp
index ab096906b49..1027ecacaf7 100644
--- a/test/unit/math/prim/prob/skew_double_exponential_cdf_log_test.cpp
+++ b/test/unit/math/prim/prob/skew_double_exponential_cdf_log_test.cpp
@@ -2,7 +2,6 @@
 #include <vector>
 #include <gtest/gtest.h>
 
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/tools/promotion.hpp>
 
diff --git a/test/unit/math/prim/prob/skew_double_exponential_test.cpp b/test/unit/math/prim/prob/skew_double_exponential_test.cpp
index 9132a2f188c..056d6496984 100644
--- a/test/unit/math/prim/prob/skew_double_exponential_test.cpp
+++ b/test/unit/math/prim/prob/skew_double_exponential_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <limits>
 #include <vector>
 
@@ -35,7 +35,7 @@ TEST(ProbDistributionsSkewedDoubleExponential, errorCheck) {
 }
 
 TEST(ProbDistributionsSkewedDoubleExponential, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::skew_double_exponential_rng(10.0, 2.0, .1, rng));
 
   EXPECT_THROW(stan::math::skew_double_exponential_rng(10.0, 2.0, -1.0, rng),
@@ -60,7 +60,7 @@ TEST(ProbDistributionsSkewedDoubleExponential, test_sampling_icdf) {
 }
 
 TEST(ProbDistributionsSkewedDoubleExponential, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/skew_normal_test.cpp b/test/unit/math/prim/prob/skew_normal_test.cpp
index 1fe061f3392..45616dddddd 100644
--- a/test/unit/math/prim/prob/skew_normal_test.cpp
+++ b/test/unit/math/prim/prob/skew_normal_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -47,7 +47,7 @@ TEST(ProbDistributionsSkewNormal, distributionTest) {
 }
 
 TEST(ProbDistributionsSkewNormal, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::skew_normal_rng(10.0, 2.0, 1.0, rng));
 
   EXPECT_THROW(stan::math::skew_normal_rng(10.0, -2.0, 1.0, rng),
@@ -61,7 +61,7 @@ TEST(ProbDistributionsSkewNormal, error_check) {
 }
 
 TEST(ProbDistributionsSkewNormal, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/std_normal_test.cpp b/test/unit/math/prim/prob/std_normal_test.cpp
index 3ed4aa6500f..f7c8a31c1d7 100644
--- a/test/unit/math/prim/prob/std_normal_test.cpp
+++ b/test/unit/math/prim/prob/std_normal_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -68,12 +68,12 @@ TEST(ProbDistributionsStdNormal, distributionTest) {
 }
 
 TEST(ProbDistributionsStdNormal, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::std_normal_rng(rng));
 }
 
 TEST(ProbDistributionsStdNormal, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/student_t_test.cpp b/test/unit/math/prim/prob/student_t_test.cpp
index 6d95ebbf347..4ffae5fbf3b 100644
--- a/test/unit/math/prim/prob/student_t_test.cpp
+++ b/test/unit/math/prim/prob/student_t_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -48,7 +48,7 @@ TEST(ProbDistributionsStudentT, distributionTest) {
 }
 
 TEST(ProbDistributionsStudentT, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::student_t_rng(3.0, 2.0, 2.0, rng));
 
   EXPECT_THROW(stan::math::student_t_rng(3.0, 2.0, -2.0, rng),
@@ -67,7 +67,7 @@ TEST(ProbDistributionsStudentT, error_check) {
 }
 
 TEST(ProbDistributionsStudentT, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/uniform_test.cpp b/test/unit/math/prim/prob/uniform_test.cpp
index e0ff8592551..9904f0986e4 100644
--- a/test/unit/math/prim/prob/uniform_test.cpp
+++ b/test/unit/math/prim/prob/uniform_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <limits>
 #include <vector>
 
@@ -63,7 +63,7 @@ TEST(ProbDistributionsUniform, distributionTest) {
 }
 
 TEST(ProbDistributionsUniform, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::uniform_rng(1.0, 2.0, rng));
 
   EXPECT_THROW(
@@ -75,7 +75,7 @@ TEST(ProbDistributionsUniform, error_check) {
 }
 
 TEST(ProbDistributionsUniform, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/util.hpp b/test/unit/math/prim/prob/util.hpp
index 553d3749fe6..62fd47dbfb2 100644
--- a/test/unit/math/prim/prob/util.hpp
+++ b/test/unit/math/prim/prob/util.hpp
@@ -16,12 +16,45 @@ inline void assert_chi_squared(const std::vector<int>& counts,
   int bins = counts.size();
   EXPECT_EQ(bins, expected.size());
 
-  double chi = 0;
+  // Merge adjacent bins until each merged bin clears the threshold.
+  constexpr double kMinExpectedCount = 5.0;
+  std::vector<double> merged_expected;
+  std::vector<int> merged_counts;
+  double acc_expected = 0;
+  int acc_counts = 0;
   for (int i = 0; i < bins; ++i) {
-    double discrepancy = expected[i] - counts[i];
-    chi += discrepancy * discrepancy / expected[i];
+    acc_expected += expected[i];
+    acc_counts += counts[i];
+    if (acc_expected >= kMinExpectedCount) {
+      merged_expected.push_back(acc_expected);
+      merged_counts.push_back(acc_counts);
+      acc_expected = 0;
+      acc_counts = 0;
+    }
+  }
+  if (acc_expected > 0) {
+    if (merged_expected.empty()) {
+      merged_expected.push_back(acc_expected);
+      merged_counts.push_back(acc_counts);
+    } else {
+      merged_expected.back() += acc_expected;
+      merged_counts.back() += acc_counts;
+    }
+  }
+
+  int merged_bins = merged_expected.size();
+  if (merged_bins < 2) {
+    // Not enough expected mass spread across bins to run a meaningful
+    // goodness-of-fit test.
+    return;
+  }
+
+  double chi = 0;
+  for (int i = 0; i < merged_bins; ++i) {
+    double discrepancy = merged_expected[i] - merged_counts[i];
+    chi += discrepancy * discrepancy / merged_expected[i];
   }
-  boost::math::chi_squared dist(bins - 1);
+  boost::math::chi_squared dist(merged_bins - 1);
   double chi_threshold = quantile(complement(dist, tolerance));
 
   EXPECT_TRUE(chi < chi_threshold);
diff --git a/test/unit/math/prim/prob/vector_rng_test_helper.hpp b/test/unit/math/prim/prob/vector_rng_test_helper.hpp
index 717ae8a7ec0..777ac42fb67 100644
--- a/test/unit/math/prim/prob/vector_rng_test_helper.hpp
+++ b/test/unit/math/prim/prob/vector_rng_test_helper.hpp
@@ -3,7 +3,6 @@
 
 #include <gtest/gtest.h>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/meta/apply_template_permutations.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
@@ -167,7 +166,7 @@ struct check_dist_throws {
   template <typename T_param1, typename T_param2, typename T_param3,
             typename T_rig>
   void operator()(const T_rig& rig) const {
-    boost::random::mt19937 rng;
+    std::mt19937 rng;
 
     T_param1 p1;
     T_param2 p2;
@@ -373,7 +372,7 @@ struct check_quantiles {
   template <typename T_param1, typename T_param2, typename T_param3,
             typename T_rig>
   void operator()(const T_rig& rig) const {
-    boost::random::mt19937 rng;
+    std::mt19937 rng;
     T_param1 p1;
     T_param2 p2;
     T_param3 p3;
@@ -533,7 +532,7 @@ struct check_counts {
   template <typename T_param1, typename T_param2, typename T_param3,
             typename T_rig>
   void operator()(const T_rig& rig) const {
-    boost::random::mt19937 rng;
+    std::mt19937 rng;
     T_param1 p1;
     T_param2 p2;
     T_param3 p3;
diff --git a/test/unit/math/prim/prob/von_mises_test.cpp b/test/unit/math/prim/prob/von_mises_test.cpp
index d5c98f2d176..ef28059144a 100644
--- a/test/unit/math/prim/prob/von_mises_test.cpp
+++ b/test/unit/math/prim/prob/von_mises_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/util.hpp>
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <vector>
 
@@ -30,7 +30,7 @@ TEST(ProbDistributionsVonMises, errorCheck) {
  * the distributions manually
  */
 TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
 
   std::vector<double> loc
@@ -75,7 +75,7 @@ TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest) {
 }
 
 TEST(ProbDistributionsVonMises, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::von_mises_rng(1.0, 2.0, rng));
   EXPECT_NO_THROW(stan::math::von_mises_rng(1.0, 0.0, rng));
 
@@ -90,7 +90,7 @@ TEST(ProbDistributionsVonMises, error_check) {
 }
 
 TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest1) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = 80;
   boost::math::chi_squared mydist(K - 1);
@@ -138,7 +138,7 @@ TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest1) {
 }
 
 TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest2) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = 80;
   boost::math::chi_squared mydist(K - 1);
@@ -181,7 +181,7 @@ TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest2) {
 }
 
 TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest3) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = 80;
   boost::math::chi_squared mydist(K - 1);
@@ -226,7 +226,7 @@ TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest3) {
 }
 
 TEST(ProbDistributionsVonMises, chiSquareGoodnessFitTest4) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/weibull_test.cpp b/test/unit/math/prim/prob/weibull_test.cpp
index 19308aea778..9b8d64ff6be 100644
--- a/test/unit/math/prim/prob/weibull_test.cpp
+++ b/test/unit/math/prim/prob/weibull_test.cpp
@@ -2,7 +2,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/util.hpp>
 #include <boost/math/distributions.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <gtest/gtest.h>
 #include <limits>
 #include <vector>
@@ -46,7 +46,7 @@ TEST(ProbDistributionsWeibull, distributionTest) {
 }
 
 TEST(ProbDistributionsWeibull, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   EXPECT_NO_THROW(stan::math::weibull_rng(2.0, 3.0, rng));
 
   EXPECT_THROW(stan::math::weibull_rng(-2.0, 3.0, rng), std::domain_error);
@@ -57,7 +57,7 @@ TEST(ProbDistributionsWeibull, error_check) {
 }
 
 TEST(ProbDistributionsWeibull, chiSquareGoodnessFitTest) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int N = 10000;
   int K = stan::math::round(2 * std::pow(N, 0.4));
 
diff --git a/test/unit/math/prim/prob/wishart_cholesky_rng_test.cpp b/test/unit/math/prim/prob/wishart_cholesky_rng_test.cpp
index dab129f4ab9..f9bbf58fe5f 100644
--- a/test/unit/math/prim/prob/wishart_cholesky_rng_test.cpp
+++ b/test/unit/math/prim/prob/wishart_cholesky_rng_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions/chi_squared.hpp>
 #include <gtest/gtest.h>
 #include <stdexcept>
@@ -13,7 +13,7 @@ TEST(ProbDistributionsWishartCholesky, rng) {
 
   using stan::math::wishart_cholesky_rng;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   MatrixXd omega(3, 4);
   EXPECT_THROW(wishart_cholesky_rng(3.0, omega, rng), std::invalid_argument);
@@ -34,7 +34,7 @@ TEST(ProbDistributionsWishartCholesky, rng_pos_def) {
   using Eigen::MatrixXd;
   using stan::math::wishart_cholesky_rng;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   MatrixXd Sigma(2, 2);
   MatrixXd Sigma_non_pos_def(2, 2);
@@ -57,7 +57,7 @@ TEST(ProbDistributionsWishartCholesky, marginalTwoChiSquareGoodnessFitTest) {
   using stan::math::wishart_cholesky_rng;
   using std::log;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 2.0, -3.0, 4.0, 0.0, 2.0, 0.0, 3.0;
   int N = 10000;
@@ -91,7 +91,7 @@ TEST(ProbDistributionsWishartCholesky, SpecialRNGTest) {
   using stan::math::multiply_lower_tri_self_transpose;
   using stan::math::wishart_cholesky_rng;
 
-  boost::random::mt19937 rng(1234);
+  std::mt19937 rng(1234);
 
   MatrixXd sigma(3, 3);
 
diff --git a/test/unit/math/prim/prob/wishart_cholesky_test.cpp b/test/unit/math/prim/prob/wishart_cholesky_test.cpp
index 3a9d6d4f769..563054651a5 100644
--- a/test/unit/math/prim/prob/wishart_cholesky_test.cpp
+++ b/test/unit/math/prim/prob/wishart_cholesky_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 #include <boost/math/distributions.hpp>
 
diff --git a/test/unit/math/prim/prob/wishart_rng_test.cpp b/test/unit/math/prim/prob/wishart_rng_test.cpp
index 28da632500a..d7c0593a7db 100644
--- a/test/unit/math/prim/prob/wishart_rng_test.cpp
+++ b/test/unit/math/prim/prob/wishart_rng_test.cpp
@@ -1,6 +1,6 @@
 #include <stan/math/prim.hpp>
 #include <test/unit/math/prim/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions/chi_squared.hpp>
 #include <gtest/gtest.h>
 #include <stdexcept>
@@ -10,7 +10,7 @@ TEST(ProbDistributionsWishart, rng) {
   using Eigen::MatrixXd;
   using stan::math::wishart_rng;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   MatrixXd omega(3, 4);
   EXPECT_THROW(wishart_rng(3.0, omega, rng), std::invalid_argument);
@@ -28,7 +28,7 @@ TEST(ProbDistributionsWishart, rng_pos_def) {
   using Eigen::MatrixXd;
   using stan::math::wishart_rng;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   MatrixXd Sigma(2, 2);
   MatrixXd Sigma_non_pos_def(2, 2);
@@ -48,7 +48,7 @@ TEST(ProbDistributionsWishart, rng_symmetry) {
   using stan::test::unit::expect_symmetric;
   using stan::test::unit::spd_rng;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   for (int k = 1; k < 20; ++k)
     for (double nu = k - 0.9; nu < k + 10; ++nu)
       for (int n = 0; n < 10; ++n)
@@ -63,7 +63,7 @@ TEST(ProbDistributionsWishart, marginalTwoChiSquareGoodnessFitTest) {
   using stan::math::wishart_rng;
   using std::log;
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   MatrixXd sigma(3, 3);
   sigma << 9.0, -3.0, 2.0, -3.0, 4.0, 0.0, 2.0, 0.0, 3.0;
   int N = 10000;
@@ -94,7 +94,7 @@ TEST(ProbDistributionsWishart, SpecialRNGTest) {
   using Eigen::VectorXd;
   using stan::math::wishart_rng;
 
-  boost::random::mt19937 rng(1234);
+  std::mt19937 rng(1234);
 
   MatrixXd sigma(3, 3);
 
diff --git a/test/unit/math/prim/prob/wishart_test.cpp b/test/unit/math/prim/prob/wishart_test.cpp
index 85df75934a4..d8f7c6a7fd6 100644
--- a/test/unit/math/prim/prob/wishart_test.cpp
+++ b/test/unit/math/prim/prob/wishart_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/special_functions/digamma.hpp>
 #include <boost/math/distributions.hpp>
 
diff --git a/test/unit/math/prim/prob/yule_simon_test.cpp b/test/unit/math/prim/prob/yule_simon_test.cpp
index 8935d088494..ffb3825d1d2 100644
--- a/test/unit/math/prim/prob/yule_simon_test.cpp
+++ b/test/unit/math/prim/prob/yule_simon_test.cpp
@@ -3,7 +3,7 @@
 #include <test/unit/math/prim/prob/vector_rng_test_helper.hpp>
 #include <test/unit/math/prim/prob/VectorIntRNGTestRig.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <limits>
 #include <vector>
@@ -36,7 +36,7 @@ TEST(ProbDistributionsYuleSimon, distributionCheck) {
 }
 
 TEST(ProbDistributionsYuleSimon, error_check) {
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
 
   EXPECT_NO_THROW(stan::math::yule_simon_rng(1.0, rng));
   EXPECT_NO_THROW(stan::math::yule_simon_rng(2.0, rng));
diff --git a/test/unit/math/rev/fun/cholesky_decompose_test.cpp b/test/unit/math/rev/fun/cholesky_decompose_test.cpp
index 1cec6737840..7313aeabb23 100644
--- a/test/unit/math/rev/fun/cholesky_decompose_test.cpp
+++ b/test/unit/math/rev/fun/cholesky_decompose_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/rev.hpp>
 #include <test/unit/math/rev/util.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <vector>
 
 template <typename T_x>
@@ -184,7 +184,7 @@ inline void test_gradients_simple(int size, double prec) {
 
   stan::math::welford_covar_estimator estimator(size);
 
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   for (int i = 0; i < 1000; ++i) {
     Eigen::VectorXd q(size);
     for (int j = 0; j < size; ++j)
@@ -228,7 +228,7 @@ inline void test_gp_grad(int mat_size, double prec) {
   test_vals[1] = 1.5;
   test_vals[2] = 10;
 
-  boost::random::mt19937 rng(2);
+  std::mt19937 rng(2);
 
   for (int i = 0; i < mat_size; ++i) {
     test_vec(i) = stan::math::normal_rng(0.0, 0.1, rng);
@@ -266,7 +266,7 @@ inline void test_chol_mult(int mat_size, double prec) {
   Eigen::VectorXd test_vec(mat_size);
   Eigen::VectorXd test_vals(vec_size);
 
-  boost::random::mt19937 rng(2);
+  std::mt19937 rng(2);
 
   for (int i = 0; i < test_vals.size(); ++i) {
     if (i < test_vec.size()) {
@@ -292,7 +292,7 @@ inline void test_chol_mult(int mat_size, double prec) {
 
 inline void test_simple_vec_mult(int size, double prec) {
   Eigen::VectorXd test_vec(size);
-  boost::random::mt19937 rng(2);
+  std::mt19937 rng(2);
 
   for (int i = 0; i < test_vec.size(); ++i)
     test_vec(i) = stan::math::normal_rng(0.0, 0.1, rng);
diff --git a/test/unit/math/rev/fun/promote_elements_test.cpp b/test/unit/math/rev/fun/promote_elements_test.cpp
index c60201fca8e..d18b3506f35 100644
--- a/test/unit/math/rev/fun/promote_elements_test.cpp
+++ b/test/unit/math/rev/fun/promote_elements_test.cpp
@@ -1,7 +1,6 @@
 #include <stan/math/rev.hpp>
 #include <test/unit/math/rev/util.hpp>
 #include <stan/math/prim.hpp>
-#include <boost/typeof/typeof.hpp>
 #include <gtest/gtest.h>
 #include <type_traits>
 #include <vector>
@@ -14,7 +13,7 @@ using std::vector;
 TEST_F(AgradRev, MathFunctionsScalPromote_Elements_double2var) {
   double from;
   promote_elements<var, double> p;
-  typedef BOOST_TYPEOF(p.promote(from)) result_t;
+  typedef decltype(p.promote(from)) result_t;
   bool same = std::is_same<var, result_t>::value;
   EXPECT_TRUE(same);
 }
@@ -25,7 +24,7 @@ TEST_F(AgradRev, MathFunctionsArrPromote_Elements_doubleVec2varVec) {
   from.push_back(2);
   from.push_back(3);
   promote_elements<vector<var>, vector<double> > p;
-  typedef BOOST_TYPEOF(p.promote(from)) result_t;
+  typedef decltype(p.promote(from)) result_t;
   bool same = std::is_same<vector<var>, result_t>::value;
   EXPECT_TRUE(same);
 }
@@ -34,7 +33,7 @@ TEST_F(AgradRev, MathFunctionsMatPromote_Elements_doubleMat2varMat) {
   stan::math::matrix_d m1(2, 3);
   m1 << 1, 2, 3, 4, 5, 6;
   promote_elements<Matrix<var, 2, 3>, Matrix<double, 2, 3> > p;
-  typedef BOOST_TYPEOF(p.promote(m1)) result_t;
+  typedef decltype(p.promote(m1)) result_t;
   bool same = std::is_same<Matrix<var, 2, 3>, result_t>::value;
   EXPECT_TRUE(same);
 }
diff --git a/test/unit/math/rev/prob/lkj_corr_test.cpp b/test/unit/math/rev/prob/lkj_corr_test.cpp
index 6c912bd949a..229709e4c30 100644
--- a/test/unit/math/rev/prob/lkj_corr_test.cpp
+++ b/test/unit/math/rev/prob/lkj_corr_test.cpp
@@ -2,13 +2,13 @@
 #include <test/unit/math/rev/prob/lkj_corr_cholesky_test_functors.hpp>
 #include <test/unit/math/rev/prob/test_gradients.hpp>
 #include <test/unit/math/rev/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
+#include <random>
 #include <boost/math/distributions.hpp>
 #include <gtest/gtest.h>
 
 TEST_F(AgradRev, ProbDistributionsLkjCorr_var) {
   using stan::math::var;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int K = 4;
   Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Sigma_d(K, K);
@@ -29,7 +29,7 @@ TEST_F(AgradRev, ProbDistributionsLkjCorr_var) {
 
 TEST_F(AgradRev, ProbDistributionsLkjCorrCholesky_var) {
   using stan::math::var;
-  boost::random::mt19937 rng;
+  std::mt19937 rng;
   int K = 4;
   Eigen::Matrix<var, Eigen::Dynamic, Eigen::Dynamic> Sigma(K, K);
   Sigma.setZero();
diff --git a/test/unit/math/rev/prob/skew_double_exponential_ccdf_log_test.cpp b/test/unit/math/rev/prob/skew_double_exponential_ccdf_log_test.cpp
index 64641e76981..18cc10416c3 100644
--- a/test/unit/math/rev/prob/skew_double_exponential_ccdf_log_test.cpp
+++ b/test/unit/math/rev/prob/skew_double_exponential_ccdf_log_test.cpp
@@ -1,7 +1,6 @@
 
 #include <stan/math/rev.hpp>
 #include <test/unit/math/rev/util.hpp>
-#include <boost/random/mersenne_twister.hpp>
 #include <boost/math/distributions.hpp>
 #include <boost/math/tools/promotion.hpp>
 #include <limits>
diff --git a/test/unit/util.hpp b/test/unit/util.hpp
index 6d13989d579..a01fb696103 100644
--- a/test/unit/util.hpp
+++ b/test/unit/util.hpp
@@ -2,9 +2,7 @@
 #define TEST_UNIT_UTIL_HPP
 
 #include <stan/math/prim/fun/Eigen.hpp>
-#include <boost/random/mersenne_twister.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
-#include <boost/typeof/typeof.hpp>
+#include <random>
 #include <test/unit/pretty_print_types.hpp>
 #include <gtest/gtest.h>
 #include <type_traits>
@@ -214,8 +212,8 @@ namespace test {
 
 auto make_sparse_matrix_random(int rows, int cols) {
   using eigen_triplet = Eigen::Triplet<double>;
-  boost::mt19937 gen;
-  boost::random::uniform_real_distribution<double> dist(0.0, 1.0);
+  std::mt19937 gen;
+  std::uniform_real_distribution<double> dist(0.0, 1.0);
   std::vector<eigen_triplet> tripletList;
   for (int i = 0; i < rows; ++i) {
     for (int j = 0; j < cols; ++j) {