diff --git a/src/stan/analyze/mcmc/autocovariance.hpp b/src/stan/analyze/mcmc/autocovariance.hpp
index ae360141d42..ab427662873 100644
--- a/src/stan/analyze/mcmc/autocovariance.hpp
+++ b/src/stan/analyze/mcmc/autocovariance.hpp
@@ -2,9 +2,6 @@
 #define STAN_ANALYZE_MCMC_AUTOCOVARIANCE_HPP
 
 #include <stan/math/prim.hpp>
-#include <boost/accumulators/accumulators.hpp>
-#include <boost/accumulators/statistics/stats.hpp>
-#include <boost/accumulators/statistics/variance.hpp>
 #include <unsupported/Eigen/FFT>
 #include <complex>
 #include <vector>
@@ -87,16 +84,9 @@ void autocovariance(const Eigen::MatrixBase<DerivedA>& y,
   Eigen::FFT<T> fft;
   autocorrelation(y, acov, fft);
 
-  using boost::accumulators::accumulator_set;
-  using boost::accumulators::stats;
-  using boost::accumulators::tag::variance;
+  double variance = (y.array() - y.mean()).matrix().squaredNorm() / y.size();
 
-  accumulator_set<double, stats<variance>> acc;
-  for (int n = 0; n < y.size(); ++n) {
-    acc(y(n));
-  }
-
-  acov = acov.array() * boost::accumulators::variance(acc);
+  acov = acov.array() * variance;
 }
 
 /**
diff --git a/src/stan/analyze/mcmc/compute_potential_scale_reduction.hpp b/src/stan/analyze/mcmc/compute_potential_scale_reduction.hpp
index baa686e63d8..901206c394e 100644
--- a/src/stan/analyze/mcmc/compute_potential_scale_reduction.hpp
+++ b/src/stan/analyze/mcmc/compute_potential_scale_reduction.hpp
@@ -4,11 +4,6 @@
 #include <stan/math/prim.hpp>
 #include <stan/analyze/mcmc/autocovariance.hpp>
 #include <stan/analyze/mcmc/split_chains.hpp>
-#include <boost/accumulators/accumulators.hpp>
-#include <boost/accumulators/statistics/stats.hpp>
-#include <boost/accumulators/statistics/mean.hpp>
-#include <boost/accumulators/statistics/variance.hpp>
-#include <boost/math/distributions/normal.hpp>
 #include <algorithm>
 #include <cmath>
 #include <vector>
@@ -74,30 +69,22 @@ inline double compute_potential_scale_reduction(
     }
   }
 
-  using boost::accumulators::accumulator_set;
-  using boost::accumulators::stats;
-  using boost::accumulators::tag::mean;
-  using boost::accumulators::tag::variance;
-
   Eigen::VectorXd chain_mean(num_chains);
-  accumulator_set<double, stats<variance>> acc_chain_mean;
   Eigen::VectorXd chain_var(num_chains);
-  double unbiased_var_scale = num_draws / (num_draws - 1.0);
 
   for (int chain = 0; chain < num_chains; ++chain) {
-    accumulator_set<double, stats<mean, variance>> acc_draw;
-    for (int n = 0; n < num_draws; ++n) {
-      acc_draw(draws[chain][n]);
-    }
+    Eigen::Map<const Eigen::Matrix<double, Eigen::Dynamic, 1>> draw(
+        draws[chain], num_draws);
 
-    chain_mean(chain) = boost::accumulators::mean(acc_draw);
-    acc_chain_mean(chain_mean(chain));
-    chain_var(chain)
-        = boost::accumulators::variance(acc_draw) * unbiased_var_scale;
+    chain_mean(chain) = draw.mean();
+    chain_var(chain) = (draw.array() - chain_mean(chain)).matrix().squaredNorm()
+                       / (num_draws - 1.0);
   }
 
-  double var_between = num_draws * boost::accumulators::variance(acc_chain_mean)
-                       * num_chains / (num_chains - 1);
+  double var_between
+      = num_draws
+        * (chain_mean.array() - chain_mean.mean()).matrix().squaredNorm()
+        / (num_chains - 1.0);
   double var_within = chain_var.mean();
 
   return sqrt((var_between / var_within + num_draws - 1) / num_draws);
diff --git a/src/stan/io/json/json_data_handler.hpp b/src/stan/io/json/json_data_handler.hpp
index 3e0adeb279d..fb45d7d4b0d 100644
--- a/src/stan/io/json/json_data_handler.hpp
+++ b/src/stan/io/json/json_data_handler.hpp
@@ -4,7 +4,9 @@
 #include <stan/io/json/json_error.hpp>
 #include <stan/io/json/json_handler.hpp>
 #include <stan/io/json/rapidjson_parser.hpp>
+#include <stan/io/string_utils.hpp>
 #include <stan/io/var_context.hpp>
+#include <algorithm>
 #include <cctype>
 #include <iostream>
 #include <ostream>
@@ -15,8 +17,6 @@
 #include <string>
 #include <utility>
 #include <vector>
-#include <boost/algorithm/string.hpp>
-#include <boost/regex.hpp>
 
 namespace stan {
 
@@ -146,9 +146,7 @@ class json_data_handler : public stan::json::json_handler {
     array_start_r = 0;
   }
 
-  inline std::string key_str() {
-    return boost::algorithm::join(key_stack, ".");
-  }
+  inline std::string key_str() { return stan::io::join(key_stack, "."); }
 
   std::string outer_key_str() {
     std::string result;
@@ -172,8 +170,11 @@ class json_data_handler : public stan::json::json_handler {
    *  and contain only letters, numbers, or an underscore.
    */
   bool valid_varname(const std::string& name) {
-    static const boost::regex re("[a-zA-Z][a-zA-Z0-9_]*");
-    return boost::regex_match(name, re);
+    if (name.empty() || !std::isalpha(static_cast<unsigned char>(name[0])))
+      return false;
+    return std::all_of(name.begin() + 1, name.end(), [](unsigned char c) {
+      return std::isalnum(c) || c == '_';
+    });
   }
 
   bool is_array_tuples(const std::vector<std::string>& keys) {
@@ -181,7 +182,7 @@ class json_data_handler : public stan::json::json_handler {
     std::string key;
     stack.pop_back();
     while (!stack.empty()) {
-      key = boost::algorithm::join(stack, ".");
+      key = stan::io::join(stack, ".");
       if (slot_types_map[key] == meta_type::ARRAY_OF_TUPLES)
         return true;
       stack.pop_back();
@@ -194,12 +195,12 @@ class json_data_handler : public stan::json::json_handler {
     std::string key;
     stack.pop_back();
     while (!stack.empty()) {
-      key = boost::algorithm::join(stack, ".");
+      key = stan::io::join(stack, ".");
       if (slot_dims_map.count(key) == 1)
         return slot_dims_map[key];
       stack.pop_back();
     }
-    key = boost::algorithm::join(keys, ".");
+    key = stan::io::join(keys, ".");
     if (slot_dims_map.count(key) != 1)
       unexpected_error(key, "not an array");
     return slot_dims_map[key];
@@ -210,13 +211,13 @@ class json_data_handler : public stan::json::json_handler {
     std::string key;
     stack.pop_back();
     while (!stack.empty()) {
-      key = boost::algorithm::join(stack, ".");
+      key = stan::io::join(stack, ".");
       if (slot_dims_map.count(key) == 1)
         break;
       stack.pop_back();
     }
     if (stack.empty()) {
-      key = boost::algorithm::join(key_stack, ".");
+      key = stan::io::join(key_stack, ".");
       unexpected_error(key, "ill-formed array");
     }
     slot_dims_map[key] = update;
@@ -354,8 +355,7 @@ class json_data_handler : public stan::json::json_handler {
         continue;
       }
       std::vector<size_t> all_dims;
-      std::vector<std::string> slots;
-      split(slots, var.first, boost::is_any_of("."), boost::token_compress_on);
+      std::vector<std::string> slots = stan::io::split(var.first, ".", true);
       std::string slot;
       for (size_t i = 0; i < slots.size(); ++i) {
         slot.append(slots[i]);
diff --git a/src/stan/io/random_var_context.hpp b/src/stan/io/random_var_context.hpp
index 6f1a40cd122..9875dc0d303 100644
--- a/src/stan/io/random_var_context.hpp
+++ b/src/stan/io/random_var_context.hpp
@@ -3,9 +3,9 @@
 
 #include <stan/io/var_context.hpp>
 #include <stan/io/validate_dims.hpp>
-#include <boost/random/uniform_real_distribution.hpp>
 #include <algorithm>
 #include <limits>
+#include <random>
 #include <string>
 #include <vector>
 
@@ -50,8 +50,7 @@ class random_var_context : public var_context {
       for (size_t n = 0; n < num_unconstrained_; ++n)
         unconstrained_params_[n] = 0.0;
     } else {
-      boost::random::uniform_real_distribution<double> unif(-init_radius,
-                                                            init_radius);
+      std::uniform_real_distribution<double> unif(-init_radius, init_radius);
       for (size_t n = 0; n < num_unconstrained_; ++n)
         unconstrained_params_[n] = unif(rng);
     }
diff --git a/src/stan/io/stan_csv_reader.hpp b/src/stan/io/stan_csv_reader.hpp
index b4e14bcc702..ecf40e250e4 100644
--- a/src/stan/io/stan_csv_reader.hpp
+++ b/src/stan/io/stan_csv_reader.hpp
@@ -1,8 +1,9 @@
 #ifndef STAN_IO_STAN_CSV_READER_HPP
 #define STAN_IO_STAN_CSV_READER_HPP
 
-#include <boost/algorithm/string.hpp>
+#include <stan/io/string_utils.hpp>
 #include <stan/math/prim.hpp>
+#include <algorithm>
 #include <cctype>
 #include <istream>
 #include <iostream>
@@ -15,8 +16,7 @@ namespace io {
 
 inline void prettify_stan_csv_name(std::string& variable) {
   if (variable.find_first_of(":.") != std::string::npos) {
-    std::vector<std::string> parts;
-    boost::split(parts, variable, boost::is_any_of(":"));
+    std::vector<std::string> parts = split(variable, ":");
     for (auto& part : parts) {
       int pos = part.find('.');
       if (pos > 0) {
@@ -25,7 +25,7 @@ inline void prettify_stan_csv_name(std::string& variable) {
         part += "]";
       }
     }
-    variable = boost::algorithm::join(parts, ".");
+    variable = join(parts, ".");
   }
 }
 
@@ -126,10 +126,10 @@ class stan_csv_reader {
       size_t equal = lhs.find("=");
       if (equal != std::string::npos) {
         name = lhs.substr(0, equal);
-        boost::trim(name);
+        trim(name);
         value = lhs.substr(equal + 1, lhs.size());
-        boost::trim(value);
-        boost::replace_first(value, " (Default)", "");
+        trim(value);
+        replace_first(value, " (Default)", "");
       } else {
         if (lhs.compare(" data") == 0) {
           ss >> comment;
@@ -138,9 +138,9 @@ class stan_csv_reader {
           size_t equal = lhs.find("=");
           if (equal != std::string::npos) {
             name = lhs.substr(0, equal);
-            boost::trim(name);
+            trim(name);
             value = lhs.substr(equal + 2, lhs.size());
-            boost::replace_first(value, " (Default)", "");
+            replace_first(value, " (Default)", "");
           }
 
           if (name.compare("file") == 0)
@@ -176,7 +176,7 @@ class stan_csv_reader {
         std::stringstream(value) >> metadata.chain_id;
       } else if (name.compare("init") == 0) {
         metadata.init = value;
-        boost::trim(metadata.init);
+        trim(metadata.init);
       } else if (name.compare("seed") == 0) {
         std::stringstream(value) >> metadata.seed;
         metadata.random_seed = false;
@@ -209,7 +209,7 @@ class stan_csv_reader {
     while (ss.good()) {
       std::string token;
       std::getline(ss, token, ',');
-      boost::trim(token);
+      trim(token);
 
       if (prettify_name) {
         prettify_stan_csv_name(token);
@@ -239,7 +239,7 @@ class stan_csv_reader {
 
     // parse stepsize
     std::getline(ss, line, '=');  // stepsize
-    boost::trim(line);
+    trim(line);
     ss >> adaptation.step_size;
     if (lines == 2)  // ADVI reports stepsize, no metric
       return;
@@ -265,7 +265,7 @@ class stan_csv_reader {
       for (int col = 0; col < cols; col++) {
         std::string token;
         std::getline(line_ss, token, ',');
-        boost::trim(token);
+        trim(token);
         std::stringstream(token) >> adaptation.metric(row, col);
       }
       std::getline(ss, line);
@@ -335,7 +335,7 @@ class stan_csv_reader {
         std::stringstream ls(line);
         for (int col = 0; col < cols; col++) {
           std::getline(ls, line, ',');
-          boost::trim(line);
+          trim(line);
           try {
             samples(row, col) = static_cast<double>(std::stold(line));
             // If the value read is out of the range of representable values by
diff --git a/src/stan/io/string_utils.hpp b/src/stan/io/string_utils.hpp
new file mode 100644
index 00000000000..523fda460fb
--- /dev/null
+++ b/src/stan/io/string_utils.hpp
@@ -0,0 +1,103 @@
+#ifndef STAN_IO_STRING_UTILS_HPP
+#define STAN_IO_STRING_UTILS_HPP
+
+#include <algorithm>
+#include <cctype>
+#include <string>
+#include <vector>
+
+namespace stan {
+namespace io {
+
+/**
+ * Splits a string on any character found in `delims`.
+ *
+ * If `compress_delims` is false, every delimiter produces a split, so
+ * adjacent delimiters yield an empty token between them (e.g. splitting
+ * "a,,b" on "," gives {"a", "", "b"}).
+ *
+ * If `compress_delims` is true, runs of adjacent delimiters are treated
+ * as a single delimiter, except that a leading or trailing run still
+ * yields a single empty token at that end (matching the behavior of
+ * boost::split with token_compress_on).
+ *
+ * @param s string to split
+ * @param delims set of delimiter characters
+ * @param compress_delims whether to collapse adjacent delimiters
+ * @return the tokens found in `s`
+ */
+inline std::vector<std::string> split(const std::string& s,
+                                      const std::string& delims,
+                                      bool compress_delims = false) {
+  std::vector<std::string> tokens;
+  size_t start = 0;
+  while (true) {
+    size_t pos = s.find_first_of(delims, start);
+    if (pos == std::string::npos) {
+      tokens.push_back(s.substr(start));
+      break;
+    }
+    tokens.push_back(s.substr(start, pos - start));
+    start = pos + 1;
+  }
+  if (compress_delims) {
+    std::vector<std::string> compressed;
+    for (size_t i = 0; i < tokens.size(); ++i) {
+      // drop empty tokens produced by interior runs of delimiters,
+      // keeping a single empty token at either end
+      if (tokens[i].empty() && i > 0 && i + 1 < tokens.size())
+        continue;
+      compressed.push_back(std::move(tokens[i]));
+    }
+    tokens.swap(compressed);
+  }
+  return tokens;
+}
+
+/**
+ * Joins a vector of strings into a single string, separated by `sep`.
+ *
+ * @param parts strings to join
+ * @param sep separator inserted between parts
+ * @return the joined string
+ */
+inline std::string join(const std::vector<std::string>& parts,
+                        const std::string& sep) {
+  std::string result;
+  for (size_t i = 0; i < parts.size(); ++i) {
+    if (i > 0)
+      result += sep;
+    result += parts[i];
+  }
+  return result;
+}
+
+/**
+ * Trims leading and trailing whitespace from a string, in place.
+ *
+ * @param s string to trim
+ */
+inline void trim(std::string& s) {
+  auto not_space = [](unsigned char c) { return std::isspace(c) == 0; };
+  s.erase(s.begin(), std::find_if(s.begin(), s.end(), not_space));
+  s.erase(std::find_if(s.rbegin(), s.rend(), not_space).base(), s.end());
+}
+
+/**
+ * Replaces the first occurrence of `target` in `s` with `replacement`,
+ * in place. Does nothing if `target` is not found.
+ *
+ * @param s string to modify
+ * @param target substring to replace
+ * @param replacement replacement text
+ */
+inline void replace_first(std::string& s, const std::string& target,
+                          const std::string& replacement) {
+  size_t pos = s.find(target);
+  if (pos != std::string::npos)
+    s.replace(pos, target.size(), replacement);
+}
+
+}  // namespace io
+}  // namespace stan
+#endif
diff --git a/src/stan/mcmc/chains.hpp b/src/stan/mcmc/chains.hpp
index fdccc53cf5d..3232bcd883f 100644
--- a/src/stan/mcmc/chains.hpp
+++ b/src/stan/mcmc/chains.hpp
@@ -6,14 +6,6 @@
 #include <stan/analyze/mcmc/compute_effective_sample_size.hpp>
 #include <stan/analyze/mcmc/compute_potential_scale_reduction.hpp>
 #include <stan/analyze/mcmc/split_rank_normalized_rhat.hpp>
-#include <boost/accumulators/accumulators.hpp>
-#include <boost/accumulators/statistics/stats.hpp>
-#include <boost/accumulators/statistics/mean.hpp>
-#include <boost/accumulators/statistics/tail_quantile.hpp>
-#include <boost/accumulators/statistics/p_square_quantile.hpp>
-#include <boost/accumulators/statistics/variance.hpp>
-#include <boost/accumulators/statistics/covariance.hpp>
-#include <boost/accumulators/statistics/variates/covariate.hpp>
 #include <algorithm>
 #include <cmath>
 #include <iostream>
@@ -65,109 +57,35 @@ class chains {
                            std::ostream* err = 0) {
     if (x.rows() != y.rows() && err)
       *err << "warning: covariance of different length chains";
-    using boost::accumulators::accumulator_set;
-    using boost::accumulators::stats;
-    using boost::accumulators::tag::covariance;
-    using boost::accumulators::tag::covariate1;
-    using boost::accumulators::tag::variance;
-
-    accumulator_set<double, stats<covariance<double, covariate1> > > acc;
 
     int M = std::min(x.size(), y.size());
-    for (int i = 0; i < M; i++)
-      acc(x(i), boost::accumulators::covariate1 = y(i));
-
-    return boost::accumulators::covariance(acc) * M / (M - 1);
+    double mx = x.head(M).mean();
+    double my = y.head(M).mean();
+    return ((x.head(M).array() - mx) * (y.head(M).array() - my)).sum()
+           / (M - 1.0);
   }
 
   static double correlation(const Eigen::VectorXd& x, const Eigen::VectorXd& y,
                             std::ostream* err = 0) {
     if (x.rows() != y.rows() && err)
       *err << "warning: covariance of different length chains";
-    using boost::accumulators::accumulator_set;
-    using boost::accumulators::stats;
-    using boost::accumulators::tag::covariance;
-    using boost::accumulators::tag::covariate1;
-    using boost::accumulators::tag::variance;
-
-    accumulator_set<double, stats<variance, covariance<double, covariate1> > >
-        acc_xy;
-    accumulator_set<double, stats<variance> > acc_y;
 
     int M = std::min(x.size(), y.size());
-    for (int i = 0; i < M; i++) {
-      acc_xy(x(i), boost::accumulators::covariate1 = y(i));
-      acc_y(y(i));
-    }
-
-    double cov = boost::accumulators::covariance(acc_xy);
+    double cov = covariance(x, y);
     if (cov > -1e-8 && cov < 1e-8)
       return cov;
-    return cov
-           / std::sqrt(boost::accumulators::variance(acc_xy)
-                       * boost::accumulators::variance(acc_y));
+    return cov / std::sqrt(variance(x.head(M)) * variance(y.head(M)));
   }
 
   static double quantile(const Eigen::VectorXd& x, const double prob) {
-    using boost::accumulators::accumulator_set;
-    using boost::accumulators::left;
-    using boost::accumulators::quantile;
-    using boost::accumulators::quantile_probability;
-    using boost::accumulators::right;
-    using boost::accumulators::stats;
-    using boost::accumulators::tag::tail;
-    using boost::accumulators::tag::tail_quantile;
-    double M = x.rows();
-    // size_t cache_size = std::min(prob, 1-prob)*M + 2;
-    size_t cache_size = M;
-
-    if (prob < 0.5) {
-      accumulator_set<double, stats<tail_quantile<left> > > acc(
-          tail<left>::cache_size = cache_size);
-      for (int i = 0; i < M; i++)
-        acc(x(i));
-      return quantile(acc, quantile_probability = prob);
-    }
-    accumulator_set<double, stats<tail_quantile<right> > > acc(
-        tail<right>::cache_size = cache_size);
-    for (int i = 0; i < M; i++)
-      acc(x(i));
-    return quantile(acc, quantile_probability = prob);
+    return stan::math::quantile(x, prob);
   }
 
   static Eigen::VectorXd quantiles(const Eigen::VectorXd& x,
                                    const Eigen::VectorXd& probs) {
-    using boost::accumulators::accumulator_set;
-    using boost::accumulators::left;
-    using boost::accumulators::quantile;
-    using boost::accumulators::quantile_probability;
-    using boost::accumulators::right;
-    using boost::accumulators::stats;
-    using boost::accumulators::tag::tail;
-    using boost::accumulators::tag::tail_quantile;
-    double M = x.rows();
-
-    // size_t cache_size = M/2 + 2;
-    size_t cache_size = M;  // 2 + 2;
-
-    accumulator_set<double, stats<tail_quantile<left> > > acc_left(
-        tail<left>::cache_size = cache_size);
-    accumulator_set<double, stats<tail_quantile<right> > > acc_right(
-        tail<right>::cache_size = cache_size);
-
-    for (int i = 0; i < M; i++) {
-      acc_left(x(i));
-      acc_right(x(i));
-    }
-
-    Eigen::VectorXd q(probs.size());
-    for (int i = 0; i < probs.size(); i++) {
-      if (probs(i) < 0.5)
-        q(i) = quantile(acc_left, quantile_probability = probs(i));
-      else
-        q(i) = quantile(acc_right, quantile_probability = probs(i));
-    }
-    return q;
+    std::vector<double> probs_vec(probs.data(), probs.data() + probs.size());
+    std::vector<double> q = stan::math::quantile(x, probs_vec);
+    return Eigen::Map<const Eigen::VectorXd>(q.data(), q.size());
   }
 
   static Eigen::VectorXd autocorrelation(const Eigen::VectorXd& x) {
diff --git a/src/stan/mcmc/hmc/base_hmc.hpp b/src/stan/mcmc/hmc/base_hmc.hpp
index 0d93b81a301..a2d700cfbf2 100644
--- a/src/stan/mcmc/hmc/base_hmc.hpp
+++ b/src/stan/mcmc/hmc/base_hmc.hpp
@@ -6,9 +6,9 @@
 #include <stan/callbacks/structured_writer.hpp>
 #include <stan/mcmc/base_mcmc.hpp>
 #include <stan/mcmc/hmc/hamiltonians/ps_point.hpp>
-#include <boost/random/uniform_01.hpp>
 #include <cmath>
 #include <limits>
+#include <random>
 #include <stdexcept>
 #include <string>
 #include <vector>
@@ -35,7 +35,7 @@ class base_hmc : public base_mcmc {
         integrator_(),
         hamiltonian_(model),
         rand_int_(rng),
-        rand_uniform_(rand_int_),
+        rand_uniform_(0.0, 1.0),
         nom_epsilon_(0.1),
         epsilon_(nom_epsilon_),
         epsilon_jitter_(0.0) {}
@@ -196,7 +196,9 @@ class base_hmc : public base_mcmc {
     this->epsilon_ = this->nom_epsilon_;
     if (this->epsilon_jitter_)
       this->epsilon_
-          *= 1.0 + this->epsilon_jitter_ * (2.0 * this->rand_uniform_() - 1.0);
+          *= 1.0
+             + this->epsilon_jitter_
+                   * (2.0 * this->rand_uniform_(this->rand_int_) - 1.0);
   }
 
  protected:
@@ -207,7 +209,7 @@ class base_hmc : public base_mcmc {
   BaseRNG& rand_int_;
 
   // Uniform(0, 1) RNG
-  boost::uniform_01<BaseRNG&> rand_uniform_;
+  std::uniform_real_distribution<double> rand_uniform_;
 
   double nom_epsilon_;
   double epsilon_;
diff --git a/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp
index 308d5ab26a4..a33401236f7 100644
--- a/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp
+++ b/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp
@@ -5,8 +5,7 @@
 #include <stan/math/prim.hpp>
 #include <stan/mcmc/hmc/hamiltonians/base_hamiltonian.hpp>
 #include <stan/mcmc/hmc/hamiltonians/dense_e_point.hpp>
-#include <boost/random/variate_generator.hpp>
-#include <boost/random/normal_distribution.hpp>
+#include <random>
 
 namespace stan {
 namespace mcmc {
@@ -42,13 +41,12 @@ class dense_e_metric : public base_hamiltonian<Model, dense_e_point, BaseRNG> {
 
   void sample_p(dense_e_point& z, BaseRNG& rng) {
     typedef typename stan::math::index_type<Eigen::VectorXd>::type idx_t;
-    boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
-        rand_dense_gaus(rng, boost::normal_distribution<>());
+    std::normal_distribution<> rand_dense_gaus;
 
     Eigen::VectorXd u(z.p.size());
 
     for (idx_t i = 0; i < u.size(); ++i)
-      u(i) = rand_dense_gaus();
+      u(i) = rand_dense_gaus(rng);
 
     z.p = z.inv_e_metric_.llt().matrixU().solve(u);
   }
diff --git a/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp
index 98bfee84294..988625d9a9c 100644
--- a/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp
+++ b/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp
@@ -4,8 +4,7 @@
 #include <stan/callbacks/logger.hpp>
 #include <stan/mcmc/hmc/hamiltonians/base_hamiltonian.hpp>
 #include <stan/mcmc/hmc/hamiltonians/diag_e_point.hpp>
-#include <boost/random/variate_generator.hpp>
-#include <boost/random/normal_distribution.hpp>
+#include <random>
 
 namespace stan {
 namespace mcmc {
@@ -42,11 +41,10 @@ class diag_e_metric : public base_hamiltonian<Model, diag_e_point, BaseRNG> {
   }
 
   void sample_p(diag_e_point& z, BaseRNG& rng) {
-    boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
-        rand_diag_gaus(rng, boost::normal_distribution<>());
+    std::normal_distribution<> rand_diag_gaus;
 
     for (int i = 0; i < z.p.size(); ++i)
-      z.p(i) = rand_diag_gaus() / sqrt(z.inv_e_metric_(i));
+      z.p(i) = rand_diag_gaus(rng) / sqrt(z.inv_e_metric_(i));
   }
 };
 
diff --git a/src/stan/mcmc/hmc/hamiltonians/softabs_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/softabs_metric.hpp
index 39cd987d955..2ec6be2cdd9 100644
--- a/src/stan/mcmc/hmc/hamiltonians/softabs_metric.hpp
+++ b/src/stan/mcmc/hmc/hamiltonians/softabs_metric.hpp
@@ -4,8 +4,7 @@
 #include <stan/math/mix.hpp>
 #include <stan/mcmc/hmc/hamiltonians/base_hamiltonian.hpp>
 #include <stan/mcmc/hmc/hamiltonians/softabs_point.hpp>
-#include <boost/random/variate_generator.hpp>
-#include <boost/random/normal_distribution.hpp>
+#include <random>
 
 namespace stan {
 namespace mcmc {
@@ -83,13 +82,12 @@ class softabs_metric : public base_hamiltonian<Model, softabs_point, BaseRNG> {
   }
 
   void sample_p(softabs_point& z, BaseRNG& rng) {
-    boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
-        rand_unit_gaus(rng, boost::normal_distribution<>());
+    std::normal_distribution<> rand_unit_gaus;
 
     Eigen::VectorXd a(z.p.size());
 
     for (idx_t n = 0; n < z.p.size(); ++n)
-      a(n) = sqrt(z.softabs_lambda(n)) * rand_unit_gaus();
+      a(n) = sqrt(z.softabs_lambda(n)) * rand_unit_gaus(rng);
 
     z.p = z.eigen_deco.eigenvectors() * a;
   }
diff --git a/src/stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp
index 9d9ca5467cc..65365157a27 100644
--- a/src/stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp
+++ b/src/stan/mcmc/hmc/hamiltonians/unit_e_metric.hpp
@@ -3,8 +3,7 @@
 
 #include <stan/mcmc/hmc/hamiltonians/base_hamiltonian.hpp>
 #include <stan/mcmc/hmc/hamiltonians/unit_e_point.hpp>
-#include <boost/random/variate_generator.hpp>
-#include <boost/random/normal_distribution.hpp>
+#include <random>
 
 namespace stan {
 namespace mcmc {
@@ -37,11 +36,10 @@ class unit_e_metric : public base_hamiltonian<Model, unit_e_point, BaseRNG> {
   }
 
   void sample_p(unit_e_point& z, BaseRNG& rng) {
-    boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
-        rand_unit_gaus(rng, boost::normal_distribution<>());
+    std::normal_distribution<> rand_unit_gaus;
 
     for (int i = 0; i < z.p.size(); ++i)
-      z.p(i) = rand_unit_gaus();
+      z.p(i) = rand_unit_gaus(rng);
   }
 };
 
diff --git a/src/stan/mcmc/hmc/nuts/base_nuts.hpp b/src/stan/mcmc/hmc/nuts/base_nuts.hpp
index b8005fa9c13..3b446007e2c 100644
--- a/src/stan/mcmc/hmc/nuts/base_nuts.hpp
+++ b/src/stan/mcmc/hmc/nuts/base_nuts.hpp
@@ -128,7 +128,7 @@ class base_nuts : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
       bool valid_subtree = false;
       double log_sum_weight_subtree = -std::numeric_limits<double>::infinity();
 
-      if (this->rand_uniform_() > 0.5) {
+      if (this->rand_uniform_(this->rand_int_) > 0.5) {
         // Extend the current trajectory forward
         this->z_.ps_point::operator=(z_fwd);
         rho_bck = rho;
@@ -164,7 +164,7 @@ class base_nuts : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
         z_sample = z_propose;
       } else {
         double accept_prob = std::exp(log_sum_weight_subtree - log_sum_weight);
-        if (this->rand_uniform_() < accept_prob)
+        if (this->rand_uniform_(this->rand_int_) < accept_prob)
           z_sample = z_propose;
       }
 
@@ -328,7 +328,7 @@ class base_nuts : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
     } else {
       double accept_prob
           = std::exp(log_sum_weight_final - log_sum_weight_subtree);
-      if (this->rand_uniform_() < accept_prob)
+      if (this->rand_uniform_(this->rand_int_) < accept_prob)
         z_propose = z_propose_final;
     }
 
diff --git a/src/stan/mcmc/hmc/nuts_classic/base_nuts_classic.hpp b/src/stan/mcmc/hmc/nuts_classic/base_nuts_classic.hpp
index 1eee6d75d58..03f579f4aca 100644
--- a/src/stan/mcmc/hmc/nuts_classic/base_nuts_classic.hpp
+++ b/src/stan/mcmc/hmc/nuts_classic/base_nuts_classic.hpp
@@ -84,7 +84,7 @@ class base_nuts_classic
     util.H0 = this->hamiltonian_.H(this->z_);
 
     // Sample the slice variable
-    util.log_u = std::log(this->rand_uniform_());
+    util.log_u = std::log(this->rand_uniform_(this->rand_int_));
 
     // Build a balanced binary tree until the NUTS criterion fails
     util.criterion = true;
@@ -101,7 +101,7 @@ class base_nuts_classic
       ps_point* z = 0;
       Eigen::VectorXd* rho = 0;
 
-      if (this->rand_uniform_() > 0.5) {
+      if (this->rand_uniform_(this->rand_int_) > 0.5) {
         z = &z_plus;
         rho = &rho_plus;
         util.sign = 1;
@@ -133,7 +133,7 @@ class base_nuts_classic
         subtree_prob = n_valid_subtree ? 1 : 0;
       }
 
-      if (this->rand_uniform_() < subtree_prob)
+      if (this->rand_uniform_(this->rand_int_) < subtree_prob)
         z_sample = z_propose;
 
       n_valid += n_valid_subtree;
@@ -227,7 +227,8 @@ class base_nuts_classic
       double accept_prob
           = static_cast<double>(n2) / static_cast<double>(n1 + n2);
 
-      if (util.criterion && (this->rand_uniform_() < accept_prob))
+      if (util.criterion
+          && (this->rand_uniform_(this->rand_int_) < accept_prob))
         z_propose = z_propose_right;
 
       Eigen::VectorXd& subtree_rho = left_subtree_rho;
diff --git a/src/stan/mcmc/hmc/static/base_static_hmc.hpp b/src/stan/mcmc/hmc/static/base_static_hmc.hpp
index 1d027cb08b5..0265b270267 100644
--- a/src/stan/mcmc/hmc/static/base_static_hmc.hpp
+++ b/src/stan/mcmc/hmc/static/base_static_hmc.hpp
@@ -59,7 +59,7 @@ class base_static_hmc
 
     double acceptProb = std::exp(H0 - h);
 
-    if (acceptProb < 1 && this->rand_uniform_() > acceptProb)
+    if (acceptProb < 1 && this->rand_uniform_(this->rand_int_) > acceptProb)
       this->z_.ps_point::operator=(z_init);
 
     acceptProb = acceptProb > 1 ? 1 : acceptProb;
diff --git a/src/stan/mcmc/hmc/static_uniform/base_static_uniform.hpp b/src/stan/mcmc/hmc/static_uniform/base_static_uniform.hpp
index 1672aade433..f1e0d35e9c5 100644
--- a/src/stan/mcmc/hmc/static_uniform/base_static_uniform.hpp
+++ b/src/stan/mcmc/hmc/static_uniform/base_static_uniform.hpp
@@ -4,9 +4,9 @@
 #include <stan/callbacks/logger.hpp>
 #include <stan/mcmc/hmc/base_hmc.hpp>
 #include <stan/mcmc/hmc/hamiltonians/ps_point.hpp>
-#include <boost/random/uniform_int_distribution.hpp>
 #include <cmath>
 #include <limits>
+#include <random>
 #include <string>
 #include <vector>
 
@@ -46,7 +46,7 @@ class base_static_uniform
     double sum_prob = 1;
     double sum_metro_prob = 1;
 
-    boost::random::uniform_int_distribution<> uniform(0, L_ - 1);
+    std::uniform_int_distribution<> uniform(0, L_ - 1);
     int Lp = uniform(this->rand_int_);
 
     for (int l = 0; l < Lp; ++l) {
@@ -61,7 +61,7 @@ class base_static_uniform
       sum_prob += prob;
       sum_metro_prob += prob > 1 ? 1 : prob;
 
-      if (this->rand_uniform_() < prob / sum_prob)
+      if (this->rand_uniform_(this->rand_int_) < prob / sum_prob)
         z_sample = this->z_;
     }
 
@@ -79,7 +79,7 @@ class base_static_uniform
       sum_prob += prob;
       sum_metro_prob += prob > 1 ? 1 : prob;
 
-      if (this->rand_uniform_() < prob / sum_prob)
+      if (this->rand_uniform_(this->rand_int_) < prob / sum_prob)
         z_sample = this->z_;
     }
 
diff --git a/src/stan/mcmc/hmc/xhmc/base_xhmc.hpp b/src/stan/mcmc/hmc/xhmc/base_xhmc.hpp
index 2ca2d49b1c9..078e368bb09 100644
--- a/src/stan/mcmc/hmc/xhmc/base_xhmc.hpp
+++ b/src/stan/mcmc/hmc/xhmc/base_xhmc.hpp
@@ -82,7 +82,7 @@ class base_xhmc : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
       double ave_subtree = 0;
       double log_sum_weight_subtree = -std::numeric_limits<double>::infinity();
 
-      if (this->rand_uniform_() > 0.5) {
+      if (this->rand_uniform_(this->rand_int_) > 0.5) {
         this->z_.ps_point::operator=(z_plus);
         valid_subtree = build_tree(this->depth_, z_propose, ave_subtree,
                                    log_sum_weight_subtree, H0, 1, n_leapfrog,
@@ -105,7 +105,7 @@ class base_xhmc : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
       ++(this->depth_);
 
       double accept_prob = std::exp(log_sum_weight_subtree - log_sum_weight);
-      if (this->rand_uniform_() < accept_prob)
+      if (this->rand_uniform_(this->rand_int_) < accept_prob)
         z_sample = z_propose;
 
       // Break if exhaustion criterion is satisfied
@@ -224,7 +224,7 @@ class base_xhmc : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
 
     double accept_prob
         = std::exp(log_sum_weight_right - log_sum_weight_subtree);
-    if (this->rand_uniform_() < accept_prob)
+    if (this->rand_uniform_(this->rand_int_) < accept_prob)
       z_propose = z_propose_right;
 
     return std::abs(ave_subtree) >= x_delta_;
diff --git a/src/stan/optimization/lbfgs_update.hpp b/src/stan/optimization/lbfgs_update.hpp
index ba9d7afcb62..449578a1681 100644
--- a/src/stan/optimization/lbfgs_update.hpp
+++ b/src/stan/optimization/lbfgs_update.hpp
@@ -2,7 +2,7 @@
 #define STAN_OPTIMIZATION_LBFGS_UPDATE_HPP
 
 #include <stan/math/prim/fun/Eigen.hpp>
-#include <boost/circular_buffer.hpp>
+#include <stan/util/ring_buffer.hpp>
 #include <tuple>
 #include <vector>
 
@@ -53,11 +53,10 @@ class LBFGSUpdate {
       B0fact = 1.0;
     }
 
-    // New updates are pushed to the "back" of the circular buffer
+    // New updates are pushed to the "back" of the ring buffer
     Scalar invskyk = 1.0 / skyk;
     _gammak = skyk / yk.squaredNorm();
-    _buf.push_back();
-    _buf.back() = std::tie(invskyk, yk, sk);
+    _buf.push_back() = std::tie(invskyk, yk, sk);
 
     return B0fact;
   }
@@ -72,8 +71,8 @@ class LBFGSUpdate {
    **/
   inline void search_direction(VectorT &pk, const VectorT &gk) const {
     std::vector<Scalar> alphas(_buf.size());
-    typename boost::circular_buffer<UpdateT>::const_reverse_iterator buf_rit;
-    typename boost::circular_buffer<UpdateT>::const_iterator buf_it;
+    typename stan::util::ring_buffer<UpdateT>::const_reverse_iterator buf_rit;
+    typename stan::util::ring_buffer<UpdateT>::const_iterator buf_it;
     typename std::vector<Scalar>::const_iterator alpha_it;
     typename std::vector<Scalar>::reverse_iterator alpha_rit;
 
@@ -103,7 +102,7 @@ class LBFGSUpdate {
   }
 
  protected:
-  boost::circular_buffer<UpdateT> _buf;
+  stan::util::ring_buffer<UpdateT> _buf;
   Scalar _gammak;
 };
 }  // namespace optimization
diff --git a/src/stan/services/pathfinder/multi.hpp b/src/stan/services/pathfinder/multi.hpp
index 6eecba35bdc..7cdd4154b68 100644
--- a/src/stan/services/pathfinder/multi.hpp
+++ b/src/stan/services/pathfinder/multi.hpp
@@ -17,7 +17,7 @@
 #include <tbb/parallel_for.h>
 #include <tbb/concurrent_vector.h>
 #include <tbb/concurrent_queue.h>
-#include <boost/random/discrete_distribution.hpp>
+#include <random>
 #include <string>
 #include <vector>
 
@@ -230,14 +230,12 @@ inline int pathfinder_lbfgs_multi(
   Eigen::Array<double, Eigen::Dynamic, 1> weight_vals
       = stan::services::psis::psis_weights(lp_ratios, tail_len, logger);
   stan::rng_t rng = util::create_rng(random_seed, stride_id);
-  using discrete_dist_t
-      = boost::random::discrete_distribution<Eigen::Index, double>;
-  boost::variate_generator<stan::rng_t&, discrete_dist_t> rand_psis_idx(
-      rng, discrete_dist_t(boost::iterator_range<double*>(
-               weight_vals.data(), weight_vals.data() + weight_vals.size())));
+  using discrete_dist_t = std::discrete_distribution<Eigen::Index>;
+  discrete_dist_t rand_psis_idx(weight_vals.data(),
+                                weight_vals.data() + weight_vals.size());
   Eigen::Matrix<Eigen::Index, -1, 1> psis_draw_idxs(num_multi_draws);
   for (size_t i = 0; i <= num_multi_draws - 1; ++i) {
-    psis_draw_idxs.coeffRef(i) = rand_psis_idx();
+    psis_draw_idxs.coeffRef(i) = rand_psis_idx(rng);
   }
   /**
    * The sort helps two main things
diff --git a/src/stan/services/pathfinder/single.hpp b/src/stan/services/pathfinder/single.hpp
index 7202fa6dcf2..8b644321488 100644
--- a/src/stan/services/pathfinder/single.hpp
+++ b/src/stan/services/pathfinder/single.hpp
@@ -12,10 +12,11 @@
 #include <stan/services/util/initialize.hpp>
 #include <stan/services/util/create_rng.hpp>
 #include <stan/services/util/duration_diff.hpp>
-#include <boost/circular_buffer.hpp>
+#include <stan/util/ring_buffer.hpp>
 #include <tbb/parallel_for.h>
 #include <tbb/concurrent_queue.h>
 #include <tbb/task_group.h>
+#include <random>
 #include <string>
 #include <vector>
 #include <atomic>
@@ -213,8 +214,9 @@ inline elbo_est_t est_approx_draws(LPF&& lp_fun, RNG&& rng,
                                    size_t num_samples, const EigVec& alpha,
                                    const std::string& iter_msg, Logger&& logger,
                                    bool calculate_lp = true) {
-  boost::variate_generator<stan::rng_t&, boost::normal_distribution<>>
-      rand_unit_gaus(rng, boost::normal_distribution<>());
+  std::normal_distribution<> unit_gaus_dist;
+  auto rand_unit_gaus
+      = [&rng, &unit_gaus_dist]() { return unit_gaus_dist(rng); };
   const auto num_params = taylor_approx.x_center.size();
   size_t lp_fun_calls = 0;
   Eigen::MatrixXd unit_samps
@@ -659,8 +661,8 @@ inline auto pathfinder_lbfgs_single(
                 + std::to_string(lbfgs.logp()));
   }
   int ret = 0;
-  boost::circular_buffer<Eigen::VectorXd> param_buff(max_history_size);
-  boost::circular_buffer<Eigen::VectorXd> grad_buff(max_history_size);
+  stan::util::ring_buffer<Eigen::VectorXd> param_buff(max_history_size);
+  stan::util::ring_buffer<Eigen::VectorXd> grad_buff(max_history_size);
   Eigen::VectorXd prev_params
       = Eigen::Map<Eigen::VectorXd>(cont_vector.data(), cont_vector.size());
   std::size_t history_size = 0;
diff --git a/src/stan/util/ring_buffer.hpp b/src/stan/util/ring_buffer.hpp
new file mode 100644
index 00000000000..b26d2fffedf
--- /dev/null
+++ b/src/stan/util/ring_buffer.hpp
@@ -0,0 +1,174 @@
+#ifndef STAN_UTIL_RING_BUFFER_HPP
+#define STAN_UTIL_RING_BUFFER_HPP
+
+#include <algorithm>
+#include <cstddef>
+#include <iterator>
+#include <utility>
+#include <vector>
+
+namespace stan {
+namespace util {
+
+/**
+ * A fixed-capacity ring buffer. Once full, each push evicts the oldest
+ * element. The backing storage is allocated once, up front, and never
+ * grows or shrinks on push/evict: `push_back()` hands back a reference to
+ * an already-constructed (and, once the buffer has wrapped around at least
+ * once, already appropriately-sized) element for the caller to assign into,
+ * so pushing an `Eigen`-typed element reuses its existing heap buffer
+ * instead of allocating a new one.
+ *
+ * @tparam T element type
+ */
+template <typename T>
+class ring_buffer {
+ public:
+  /**
+   * Construct a ring buffer with a fixed capacity.
+   *
+   * @param capacity maximum number of elements retained at once
+   */
+  explicit ring_buffer(size_t capacity)
+      : buf_(capacity),
+        capacity_(capacity),
+        start_(0),
+        size_(0),
+        last_idx_(0) {}
+
+  size_t size() const { return size_; }
+  size_t capacity() const { return capacity_; }
+
+  void clear() {
+    start_ = 0;
+    size_ = 0;
+  }
+
+  /**
+   * Make the next slot available for writing, evicting the oldest element
+   * if the buffer is already full, and return a reference to it. The
+   * returned slot holds whatever was last stored there (or a
+   * default-constructed `T` if this capacity has never been filled), ready
+   * to be overwritten by assignment.
+   */
+  T& push_back() {
+    size_t idx;
+    if (size_ < capacity_) {
+      idx = (start_ + size_) % capacity_;
+      ++size_;
+    } else {
+      idx = start_;
+      start_ = (start_ + 1 == capacity_) ? 0 : start_ + 1;
+    }
+    last_idx_ = idx;
+    return buf_[idx];
+  }
+
+  /**
+   * Push a value, evicting the oldest element if the buffer is full.
+   * Assigns directly into the reused slot rather than constructing a
+   * temporary, so an rvalue expression (e.g. an Eigen expression template)
+   * is evaluated straight into the slot's existing storage.
+   */
+  template <typename U>
+  void push_back(U&& value) {
+    push_back() = std::forward<U>(value);
+  }
+
+  T& back() { return buf_[last_idx_]; }
+  const T& back() const { return buf_[last_idx_]; }
+
+  T& operator[](size_t i) { return buf_[(start_ + i) % capacity_]; }
+  const T& operator[](size_t i) const { return buf_[(start_ + i) % capacity_]; }
+
+  /**
+   * Change the capacity, keeping the most-recently-pushed
+   * `min(size(), new_capacity)` elements.
+   */
+  void rset_capacity(size_t new_capacity) {
+    std::vector<T> new_buf(new_capacity);
+    size_t keep = std::min(size_, new_capacity);
+    size_t old_first = (start_ + (size_ - keep)) % capacity_;
+    for (size_t i = 0; i < keep; ++i)
+      new_buf[i] = buf_[(old_first + i) % capacity_];
+    buf_.swap(new_buf);
+    capacity_ = new_capacity;
+    start_ = 0;
+    size_ = keep;
+    last_idx_ = keep == 0 ? 0 : keep - 1;
+  }
+
+  class const_iterator {
+   public:
+    using iterator_category = std::bidirectional_iterator_tag;
+    using value_type = T;
+    using difference_type = std::ptrdiff_t;
+    using pointer = const T*;
+    using reference = const T&;
+
+    const_iterator() : buf_(nullptr), capacity_(0), idx_(0), pos_(0) {}
+    const_iterator(const std::vector<T>* buf, size_t capacity, size_t idx,
+                   size_t pos)
+        : buf_(buf), capacity_(capacity), idx_(idx), pos_(pos) {}
+
+    reference operator*() const { return (*buf_)[idx_]; }
+    pointer operator->() const { return &(*buf_)[idx_]; }
+
+    const_iterator& operator++() {
+      ++pos_;
+      idx_ = (idx_ + 1 == capacity_) ? 0 : idx_ + 1;
+      return *this;
+    }
+    const_iterator operator++(int) {
+      const_iterator tmp = *this;
+      ++*this;
+      return tmp;
+    }
+    const_iterator& operator--() {
+      --pos_;
+      idx_ = (idx_ == 0) ? capacity_ - 1 : idx_ - 1;
+      return *this;
+    }
+    const_iterator operator--(int) {
+      const_iterator tmp = *this;
+      --*this;
+      return tmp;
+    }
+
+    bool operator==(const const_iterator& o) const { return pos_ == o.pos_; }
+    bool operator!=(const const_iterator& o) const { return pos_ != o.pos_; }
+
+   private:
+    const std::vector<T>* buf_;
+    size_t capacity_;
+    size_t idx_;
+    size_t pos_;
+  };
+
+  const_iterator begin() const {
+    return const_iterator(&buf_, capacity_, start_, 0);
+  }
+  const_iterator end() const {
+    return const_iterator(&buf_, capacity_, (start_ + size_) % capacity_,
+                          size_);
+  }
+
+  using const_reverse_iterator = std::reverse_iterator<const_iterator>;
+  const_reverse_iterator rbegin() const {
+    return const_reverse_iterator(end());
+  }
+  const_reverse_iterator rend() const {
+    return const_reverse_iterator(begin());
+  }
+
+ private:
+  std::vector<T> buf_;
+  size_t capacity_;
+  size_t start_;
+  size_t size_;
+  size_t last_idx_;
+};
+
+}  // namespace util
+}  // namespace stan
+#endif
diff --git a/src/stan/variational/advi.hpp b/src/stan/variational/advi.hpp
index 681ce65e82c..579e222b41d 100644
--- a/src/stan/variational/advi.hpp
+++ b/src/stan/variational/advi.hpp
@@ -9,7 +9,7 @@
 #include <stan/variational/print_progress.hpp>
 #include <stan/variational/families/normal_fullrank.hpp>
 #include <stan/variational/families/normal_meanfield.hpp>
-#include <boost/circular_buffer.hpp>
+#include <stan/util/ring_buffer.hpp>
 #include <algorithm>
 #include <chrono>
 #include <limits>
@@ -335,7 +335,7 @@ class advi {
     // Heuristic to estimate how far to look back in rolling window
     int cb_size
         = static_cast<int>(std::max(0.1 * max_iterations / eval_elbo_, 2.0));
-    boost::circular_buffer<double> elbo_diff(cb_size);
+    stan::util::ring_buffer<double> elbo_diff(cb_size);
 
     logger.info("Begin stochastic gradient ascent.");
     logger.info(
@@ -528,10 +528,10 @@ class advi {
    * @param[in] cb circular buffer with some number of values in it.
    * @return median of values in circular buffer.
    */
-  double circ_buff_median(const boost::circular_buffer<double>& cb) const {
+  double circ_buff_median(const stan::util::ring_buffer<double>& cb) const {
     // FIXME: naive implementation; creates a copy as a vector
     std::vector<double> v;
-    for (boost::circular_buffer<double>::const_iterator i = cb.begin();
+    for (stan::util::ring_buffer<double>::const_iterator i = cb.begin();
          i != cb.end(); ++i) {
       v.push_back(*i);
     }
diff --git a/src/test/unit/callbacks/stream_logger_test.cpp b/src/test/unit/callbacks/stream_logger_test.cpp
index 706b5d65469..015d326aa95 100644
--- a/src/test/unit/callbacks/stream_logger_test.cpp
+++ b/src/test/unit/callbacks/stream_logger_test.cpp
@@ -1,5 +1,4 @@
 #include <gtest/gtest.h>
-#include <boost/lexical_cast.hpp>
 #include <sstream>
 #include <stan/callbacks/stream_logger.hpp>
 
diff --git a/src/test/unit/callbacks/stream_writer_test.cpp b/src/test/unit/callbacks/stream_writer_test.cpp
index e0bc01bb9be..cfcbf765949 100644
--- a/src/test/unit/callbacks/stream_writer_test.cpp
+++ b/src/test/unit/callbacks/stream_writer_test.cpp
@@ -1,5 +1,4 @@
 #include <gtest/gtest.h>
-#include <boost/lexical_cast.hpp>
 #include <stan/callbacks/stream_writer.hpp>
 
 class StanInterfaceCallbacksStreamWriter : public ::testing::Test {
@@ -32,7 +31,7 @@ TEST_F(StanInterfaceCallbacksStreamWriter, string_vector) {
   const int N = 5;
   std::vector<std::string> x;
   for (int n = 0; n < N; ++n)
-    x.push_back(boost::lexical_cast<std::string>(n));
+    x.push_back(std::to_string(n));
 
   EXPECT_NO_THROW(writer(x));
   EXPECT_EQ("0,1,2,3,4\n", ss.str());
diff --git a/src/test/unit/callbacks/unique_stream_writer_test.cpp b/src/test/unit/callbacks/unique_stream_writer_test.cpp
index 0727313b4c1..f80b68964e9 100644
--- a/src/test/unit/callbacks/unique_stream_writer_test.cpp
+++ b/src/test/unit/callbacks/unique_stream_writer_test.cpp
@@ -1,5 +1,4 @@
 #include <gtest/gtest.h>
-#include <boost/lexical_cast.hpp>
 #include <stan/callbacks/unique_stream_writer.hpp>
 
 struct deleter_noop {
@@ -51,7 +50,7 @@ TEST_F(StanInterfaceCallbacksStreamWriter, string_vector) {
   const int N = 5;
   std::vector<std::string> x;
   for (int n = 0; n < N; ++n)
-    x.push_back(boost::lexical_cast<std::string>(n));
+    x.push_back(std::to_string(n));
 
   EXPECT_NO_THROW(writer(x));
   EXPECT_EQ("0,1,2,3,4\n", ss.str());
diff --git a/src/test/unit/callbacks/writer_test.cpp b/src/test/unit/callbacks/writer_test.cpp
index 730f95f1e19..9934d9f6bae 100644
--- a/src/test/unit/callbacks/writer_test.cpp
+++ b/src/test/unit/callbacks/writer_test.cpp
@@ -1,5 +1,4 @@
 #include <gtest/gtest.h>
-#include <boost/lexical_cast.hpp>
 #include <stan/callbacks/writer.hpp>
 
 class StanInterfaceCallbacksWriter : public ::testing::Test {
@@ -22,7 +21,7 @@ TEST_F(StanInterfaceCallbacksWriter, string_vector) {
   const int N = 5;
   std::vector<std::string> x;
   for (int n = 0; n < N; ++n)
-    x.push_back(boost::lexical_cast<std::string>(n));
+    x.push_back(std::to_string(n));
 
   EXPECT_NO_THROW(writer(x));
 }
diff --git a/src/test/unit/io/json/json_data_test.cpp b/src/test/unit/io/json/json_data_test.cpp
index b8ff6a9188a..b6cae2dfe18 100644
--- a/src/test/unit/io/json/json_data_test.cpp
+++ b/src/test/unit/io/json/json_data_test.cpp
@@ -7,12 +7,12 @@
 #include <test/unit/util.hpp>
 #include <test/unit/io/json/util.hpp>
 
-#include <boost/limits.hpp>
 #include <boost/math/concepts/real_concept.hpp>
 #include <boost/math/special_functions/fpclassify.hpp>
 #include <gtest/gtest.h>
 
 #include <complex>
+#include <limits>
 
 TEST(ioJson, jsonData_scalar_int) {
   std::string txt = "{ \"foo\" : 1 }";
diff --git a/src/test/unit/io/random_var_context_test.cpp b/src/test/unit/io/random_var_context_test.cpp
index 64ea370f17a..f39ca681059 100644
--- a/src/test/unit/io/random_var_context_test.cpp
+++ b/src/test/unit/io/random_var_context_test.cpp
@@ -2,7 +2,6 @@
 #include <stan/io/empty_var_context.hpp>
 #include <stan/services/util/create_rng.hpp>
 #include <gtest/gtest.h>
-#include <boost/random/uniform_real_distribution.hpp>
 #include <test/test-models/good/services/test_lp.hpp>
 #include <test/unit/util.hpp>
 
diff --git a/src/test/unit/io/stan_csv_reader_test.cpp b/src/test/unit/io/stan_csv_reader_test.cpp
index f913615732d..4f77e3a03e1 100644
--- a/src/test/unit/io/stan_csv_reader_test.cpp
+++ b/src/test/unit/io/stan_csv_reader_test.cpp
@@ -1,6 +1,5 @@
 #include <stan/io/stan_csv_reader.hpp>
 #include <test/unit/util.hpp>
-#include <boost/algorithm/string/predicate.hpp>
 #include <gtest/gtest.h>
 #include <fstream>
 #include <sstream>
diff --git a/src/test/unit/mcmc/hmc/base_hmc_test.cpp b/src/test/unit/mcmc/hmc/base_hmc_test.cpp
index ce3d8f51f45..e176fadeac7 100644
--- a/src/test/unit/mcmc/hmc/base_hmc_test.cpp
+++ b/src/test/unit/mcmc/hmc/base_hmc_test.cpp
@@ -2,7 +2,6 @@
 #include <stan/callbacks/stream_writer.hpp>
 #include <stan/mcmc/hmc/base_hmc.hpp>
 #include <stan/services/util/create_rng.hpp>
-#include <boost/algorithm/string/split.hpp>
 #include <test/unit/util.hpp>
 #include <gtest/gtest.h>
 
diff --git a/src/test/unit/mcmc/hmc/nuts/softabs_nuts_test.cpp b/src/test/unit/mcmc/hmc/nuts/softabs_nuts_test.cpp
index 8eaa84df6b3..22fdae99477 100644
--- a/src/test/unit/mcmc/hmc/nuts/softabs_nuts_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts/softabs_nuts_test.cpp
@@ -338,15 +338,15 @@ TEST(McmcSoftAbsNuts, transition_test) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_EQ(3, sampler.depth_);
-  EXPECT_EQ((2 << 3) - 1, sampler.n_leapfrog_);
+  EXPECT_EQ(5, sampler.depth_);
+  EXPECT_EQ(31, sampler.n_leapfrog_);
   EXPECT_FALSE(sampler.divergent_);
 
-  EXPECT_FLOAT_EQ(0.74693149, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.74414188, s.cont_params()(1));
-  EXPECT_FLOAT_EQ(0.60859376, s.cont_params()(2));
-  EXPECT_FLOAT_EQ(-0.74102008, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99934167, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.2057635, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(0.87303215, s.cont_params()(1));
+  EXPECT_FLOAT_EQ(0.21624902, s.cont_params()(2));
+  EXPECT_FLOAT_EQ(-0.42564371, s.log_prob());
+  EXPECT_FLOAT_EQ(0.99925238, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
diff --git a/src/test/unit/mcmc/hmc/nuts/unit_e_nuts_test.cpp b/src/test/unit/mcmc/hmc/nuts/unit_e_nuts_test.cpp
index 6e5bade65af..a2281c64497 100644
--- a/src/test/unit/mcmc/hmc/nuts/unit_e_nuts_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts/unit_e_nuts_test.cpp
@@ -338,15 +338,15 @@ TEST(McmcUnitENuts, transition_test) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_EQ(3, sampler.depth_);
-  EXPECT_EQ((2 << 3) - 1, sampler.n_leapfrog_);
+  EXPECT_EQ(4, sampler.depth_);
+  EXPECT_EQ(31, sampler.n_leapfrog_);
   EXPECT_FALSE(sampler.divergent_);
 
-  EXPECT_FLOAT_EQ(0.70149082, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.69831347, s.cont_params()(1));
-  EXPECT_FLOAT_EQ(0.54392564, s.cont_params()(2));
-  EXPECT_FLOAT_EQ(-0.63779306, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99912512, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.1622753, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(0.94628012, s.cont_params()(1));
+  EXPECT_FLOAT_EQ(0.17305008, s.cont_params()(2));
+  EXPECT_FLOAT_EQ(-0.47586286, s.log_prob());
+  EXPECT_FLOAT_EQ(0.99870592, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
diff --git a/src/test/unit/mcmc/hmc/static_uniform/derived_static_uniform_test.cpp b/src/test/unit/mcmc/hmc/static_uniform/derived_static_uniform_test.cpp
index 15527eaa8bd..f014b5831e8 100644
--- a/src/test/unit/mcmc/hmc/static_uniform/derived_static_uniform_test.cpp
+++ b/src/test/unit/mcmc/hmc/static_uniform/derived_static_uniform_test.cpp
@@ -41,9 +41,9 @@ TEST(McmcStaticUniform, unit_e_transition) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_FLOAT_EQ(1.0920367, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.59627211, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99985325, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.9501853, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(-0.45142606, s.log_prob());
+  EXPECT_FLOAT_EQ(1, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
@@ -78,9 +78,9 @@ TEST(McmcStaticUniform, diag_e_transition) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_FLOAT_EQ(1.0920367, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.59627211, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99985325, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.9501853, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(-0.45142606, s.log_prob());
+  EXPECT_FLOAT_EQ(1, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
@@ -115,9 +115,9 @@ TEST(McmcStaticUniform, dense_e_transition) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_FLOAT_EQ(1.0920367, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.59627211, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99985325, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.9501853, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(-0.45142606, s.log_prob());
+  EXPECT_FLOAT_EQ(1, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
@@ -152,9 +152,9 @@ TEST(McmcStaticUniform, softabs_transition) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_FLOAT_EQ(1.0826443, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.58605933, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99989599, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.95708251, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(-0.45800349, s.log_prob());
+  EXPECT_FLOAT_EQ(1, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
@@ -189,9 +189,9 @@ TEST(McmcStaticUniform, adapt_unit_e_transition) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_FLOAT_EQ(1.0920367, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.59627211, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99985325, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.9501853, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(-0.45142606, s.log_prob());
+  EXPECT_FLOAT_EQ(1, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
@@ -226,9 +226,9 @@ TEST(McmcStaticUniform, adapt_diag_e_transition) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_FLOAT_EQ(1.0920367, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.59627211, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99985325, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.9501853, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(-0.45142606, s.log_prob());
+  EXPECT_FLOAT_EQ(1, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
@@ -263,9 +263,9 @@ TEST(McmcStaticUniform, adapt_dense_e_transition) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_FLOAT_EQ(1.0920367, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.59627211, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99985325, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.9501853, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(-0.45142606, s.log_prob());
+  EXPECT_FLOAT_EQ(1, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
@@ -300,9 +300,9 @@ TEST(McmcStaticUniform, adapt_softabs_e_transition) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_FLOAT_EQ(1.0826443, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-0.58605933, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99989599, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.95708251, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(-0.45800349, s.log_prob());
+  EXPECT_FLOAT_EQ(1, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
diff --git a/src/test/unit/mcmc/hmc/xhmc/softabs_xhmc_test.cpp b/src/test/unit/mcmc/hmc/xhmc/softabs_xhmc_test.cpp
index faa96200538..459ff598a9f 100644
--- a/src/test/unit/mcmc/hmc/xhmc/softabs_xhmc_test.cpp
+++ b/src/test/unit/mcmc/hmc/xhmc/softabs_xhmc_test.cpp
@@ -112,7 +112,7 @@ TEST(McmcUnitEXHMC, transition) {
   EXPECT_FLOAT_EQ(-1, s.cont_params()(1));
   EXPECT_FLOAT_EQ(1, s.cont_params()(2));
   EXPECT_FLOAT_EQ(-1.5, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99870497, s.accept_stat());
+  EXPECT_FLOAT_EQ(0.99980384, s.accept_stat());
 
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
diff --git a/src/test/unit/mcmc/hmc/xhmc/unit_e_xhmc_test.cpp b/src/test/unit/mcmc/hmc/xhmc/unit_e_xhmc_test.cpp
index 61639c32475..533e84f0998 100644
--- a/src/test/unit/mcmc/hmc/xhmc/unit_e_xhmc_test.cpp
+++ b/src/test/unit/mcmc/hmc/xhmc/unit_e_xhmc_test.cpp
@@ -108,11 +108,11 @@ TEST(McmcUnitEXHMC, transition) {
 
   stan::mcmc::sample s = sampler.transition(init_sample, logger);
 
-  EXPECT_FLOAT_EQ(1, s.cont_params()(0));
-  EXPECT_FLOAT_EQ(-1, s.cont_params()(1));
-  EXPECT_FLOAT_EQ(1, s.cont_params()(2));
-  EXPECT_FLOAT_EQ(-1.5, s.log_prob());
-  EXPECT_FLOAT_EQ(0.99870926, s.accept_stat());
+  EXPECT_FLOAT_EQ(1.7558961, s.cont_params()(0));
+  EXPECT_FLOAT_EQ(-0.99812794, s.cont_params()(1));
+  EXPECT_FLOAT_EQ(0.62731504, s.cont_params()(2));
+  EXPECT_FLOAT_EQ(-2.2364774, s.log_prob());
+  EXPECT_FLOAT_EQ(0.99912089, s.accept_stat());
   EXPECT_EQ("", debug.str());
   EXPECT_EQ("", info.str());
   EXPECT_EQ("", warn.str());
diff --git a/src/test/unit/services/check_adaptation.hpp b/src/test/unit/services/check_adaptation.hpp
index af863ea8139..a242be61003 100644
--- a/src/test/unit/services/check_adaptation.hpp
+++ b/src/test/unit/services/check_adaptation.hpp
@@ -3,11 +3,11 @@
 
 #include <test/unit/services/instrumented_callbacks.hpp>
 #include <gtest/gtest.h>
+#include <stan/io/string_utils.hpp>
 #include <stdlib.h>
 #include <stdio.h>
 #include <string>
 #include <vector>
-#include <boost/algorithm/string.hpp>
 
 namespace stan {
 namespace test {
@@ -28,9 +28,8 @@ void check_adaptation(const size_t& num_params,
       break;
     }
   }
-  std::vector<std::string> strs;
-  boost::split(strs, param_strings[offset], boost::is_any_of(", "),
-               boost::token_compress_on);
+  std::vector<std::string> strs
+      = stan::io::split(param_strings[offset], ", ", true);
   EXPECT_EQ(num_params, strs.size());
   for (size_t i = 0; i < num_params; i++) {
     ASSERT_NEAR(param_vals[i], test::unit::stod(strs[i]), err_margin);
@@ -51,9 +50,8 @@ void check_adaptation(const size_t& num_rows, const size_t& num_cols,
     }
   }
   for (size_t i = 0, ij = 0; i < num_rows; i++) {
-    std::vector<std::string> strs;
-    boost::split(strs, param_strings[offset + i], boost::is_any_of(", "),
-                 boost::token_compress_on);
+    std::vector<std::string> strs
+        = stan::io::split(param_strings[offset + i], ", ", true);
     EXPECT_EQ(num_cols, strs.size());
     for (size_t j = 0; j < num_cols; j++, ij++) {
       ASSERT_NEAR(param_vals[ij], test::unit::stod(strs[j]), err_margin);
@@ -74,9 +72,8 @@ void check_different(const size_t& num_params,
       break;
     }
   }
-  std::vector<std::string> strs;
-  boost::split(strs, param_strings[offset], boost::is_any_of(", "),
-               boost::token_compress_on);
+  std::vector<std::string> strs
+      = stan::io::split(param_strings[offset], ", ", true);
   EXPECT_EQ(num_params, strs.size());
   for (size_t i = 0; i < num_params; i++) {
     ASSERT_GT(fabs(param_vals[i] - test::unit::stod(strs[i])), margin);
@@ -97,9 +94,8 @@ void check_different(const size_t& num_rows, const size_t& num_cols,
     }
   }
   for (size_t i = 0, ij = 0; i < num_rows; i++) {
-    std::vector<std::string> strs;
-    boost::split(strs, param_strings[offset + i], boost::is_any_of(", "),
-                 boost::token_compress_on);
+    std::vector<std::string> strs
+        = stan::io::split(param_strings[offset + i], ", ", true);
     EXPECT_EQ(num_cols, strs.size());
     for (size_t j = 0; j < num_cols; j++, ij++) {
       ASSERT_GT(fabs(param_vals[ij] - test::unit::stod(strs[j])), margin);
diff --git a/src/test/unit/services/optimize/laplace_jacobian_test.cpp b/src/test/unit/services/optimize/laplace_jacobian_test.cpp
index 7abeadc9455..874bef71b21 100644
--- a/src/test/unit/services/optimize/laplace_jacobian_test.cpp
+++ b/src/test/unit/services/optimize/laplace_jacobian_test.cpp
@@ -1,4 +1,3 @@
-#include <boost/algorithm/string.hpp>
 #include <gtest/gtest.h>
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/callbacks/stream_writer.hpp>
diff --git a/src/test/unit/services/optimize/laplace_sample_test.cpp b/src/test/unit/services/optimize/laplace_sample_test.cpp
index e292305cb70..a90c133f99f 100644
--- a/src/test/unit/services/optimize/laplace_sample_test.cpp
+++ b/src/test/unit/services/optimize/laplace_sample_test.cpp
@@ -1,4 +1,3 @@
-#include <boost/algorithm/string.hpp>
 #include <gtest/gtest.h>
 #include <stan/callbacks/stream_logger.hpp>
 #include <stan/callbacks/stream_writer.hpp>
diff --git a/src/test/unit/services/pathfinder/eight_schools_test.cpp b/src/test/unit/services/pathfinder/eight_schools_test.cpp
index cd901623040..3740b9a8bd4 100644
--- a/src/test/unit/services/pathfinder/eight_schools_test.cpp
+++ b/src/test/unit/services/pathfinder/eight_schools_test.cpp
@@ -142,7 +142,7 @@ TEST_F(ServicesPathfinderEightSchools, multi) {
   all_mean_vals.row(2) = mean_vals - r_mean_vals;
   // This samples badly, but is a known issue with initialization.
   for (Eigen::Index i = 0; i < all_mean_vals.cols(); i++) {
-    EXPECT_NEAR(0, all_mean_vals(2, i), 1);
+    EXPECT_NEAR(0, all_mean_vals(2, i), 1.5);
   }
 
   Eigen::MatrixXd all_sd_vals(3, 20);
diff --git a/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_test.cpp b/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_test.cpp
index bfcdc14d6a7..740496d85eb 100644
--- a/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_test.cpp
+++ b/src/test/unit/services/sample/hmc_nuts_diag_e_adapt_test.cpp
@@ -3,7 +3,7 @@
 #include <stan/io/empty_var_context.hpp>
 #include <test/test-models/good/optimization/rosenbrock.hpp>
 #include <test/unit/services/instrumented_callbacks.hpp>
-#include <boost/algorithm/string.hpp>
+#include <stan/io/string_utils.hpp>
 #include <iostream>
 
 class ServicesSampleHmcNutsDiagEAdapt : public testing::Test {
@@ -274,8 +274,7 @@ TEST_F(ServicesSampleHmcNutsDiagEAdapt, term_buffer_1) {
   std::vector<std::string> messages = parameter.string_values();
   for (auto msg : messages) {
     if (msg.find("Step size") != std::string::npos) {
-      std::vector<std::string> toks;
-      boost::split(toks, msg, boost::is_any_of(" "));
+      std::vector<std::string> toks = stan::io::split(msg, " ");
       auto adapted = std::stod(toks[toks.size() - 1]);
       EXPECT_NEAR(draw[2], adapted, 1e-5);
     }
@@ -360,8 +359,7 @@ TEST_F(ServicesSampleHmcNutsDiagEAdapt, schedule_a) {
   std::vector<std::string> messages = parameter.string_values();
   for (auto msg : messages) {
     if (msg.find("Step size") != std::string::npos) {
-      std::vector<std::string> toks;
-      boost::split(toks, msg, boost::is_any_of(" "));
+      std::vector<std::string> toks = stan::io::split(msg, " ");
       auto adapted = std::stod(toks[toks.size() - 1]);
       EXPECT_NEAR(draw[2], adapted, 1e-5);
     }
@@ -406,8 +404,7 @@ TEST_F(ServicesSampleHmcNutsDiagEAdapt, schedule_b) {
   std::vector<std::string> messages = parameter.string_values();
   for (auto msg : messages) {
     if (msg.find("Step size") != std::string::npos) {
-      std::vector<std::string> toks;
-      boost::split(toks, msg, boost::is_any_of(" "));
+      std::vector<std::string> toks = stan::io::split(msg, " ");
       auto adapted = std::stod(toks[toks.size() - 1]);
       EXPECT_NEAR(draw[2], adapted, 1e-5);
     }
@@ -452,8 +449,7 @@ TEST_F(ServicesSampleHmcNutsDiagEAdapt, schedule_c) {
   std::vector<std::string> messages = parameter.string_values();
   for (auto msg : messages) {
     if (msg.find("Step size") != std::string::npos) {
-      std::vector<std::string> toks;
-      boost::split(toks, msg, boost::is_any_of(" "));
+      std::vector<std::string> toks = stan::io::split(msg, " ");
       auto adapted = std::stod(toks[toks.size() - 1]);
       EXPECT_NEAR(draw[2], adapted, 1e-5);
     }
diff --git a/src/test/unit/services/sample/standalone_gqs_parallel_test.cpp b/src/test/unit/services/sample/standalone_gqs_parallel_test.cpp
index 5bdc8a3f21e..517d1e41236 100644
--- a/src/test/unit/services/sample/standalone_gqs_parallel_test.cpp
+++ b/src/test/unit/services/sample/standalone_gqs_parallel_test.cpp
@@ -1,4 +1,3 @@
-#include <boost/algorithm/string.hpp>
 #include <gtest/gtest.h>
 #include <iostream>
 #include <stan/callbacks/stream_logger.hpp>
diff --git a/src/test/unit/services/sample/standalone_gqs_test.cpp b/src/test/unit/services/sample/standalone_gqs_test.cpp
index 653b063e840..d8354d55b3b 100644
--- a/src/test/unit/services/sample/standalone_gqs_test.cpp
+++ b/src/test/unit/services/sample/standalone_gqs_test.cpp
@@ -1,4 +1,3 @@
-#include <boost/algorithm/string.hpp>
 #include <gtest/gtest.h>
 #include <iostream>
 #include <stan/callbacks/stream_logger.hpp>
diff --git a/src/test/unit/util.hpp b/src/test/unit/util.hpp
index f47ddcad39a..550a41fffbb 100644
--- a/src/test/unit/util.hpp
+++ b/src/test/unit/util.hpp
@@ -2,8 +2,8 @@
 #define TEST_UNIT_UTIL_HPP
 
 #include <stan/io/stan_csv_reader.hpp>
+#include <stan/io/string_utils.hpp>
 
-#include <boost/algorithm/string.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <gtest/gtest.h>
 #include <rapidjson/document.h>
@@ -48,8 +48,7 @@ void match_csv_columns(const Eigen::MatrixXd& samples,
     if (row == num_rows + 1) {
       break;
     }
-    cells.clear();
-    boost::algorithm::split(cells, line, boost::is_any_of(","));
+    cells = stan::io::split(line, ",");
     for (size_t i = 0; i < num_columns; ++i) {
       cell_ss.str(std::string());
       cell_ss.clear();