diff --git a/src/stan/analyze/mcmc/compute_effective_sample_size.hpp b/src/stan/analyze/mcmc/compute_effective_sample_size.hpp
index 39eb636c71..bdfbbfd142 100644
--- a/src/stan/analyze/mcmc/compute_effective_sample_size.hpp
+++ b/src/stan/analyze/mcmc/compute_effective_sample_size.hpp
@@ -12,6 +12,8 @@
 namespace stan {
 namespace analyze {
 /**
+ * \deprecated use split_rank_normalized_ess instead
+ *
  * Computes the effective sample size (ESS) for the specified
  * parameter across all kept samples.  The value returned is the
  * minimum of ESS and the number_total_draws *
@@ -29,6 +31,11 @@ namespace analyze {
  * @param sizes stores sizes of chains
  * @return effective sample size for the specified parameter
  */
+#if defined(__GNUC__) || defined(__clang__)
+__attribute__((deprecated))
+#elif defined(_MSC_VER)
+__declspec(deprecated)
+#endif
 inline double compute_effective_sample_size(std::vector<const double*> draws,
                                             std::vector<size_t> sizes) {
   int num_chains = sizes.size();
@@ -138,6 +145,8 @@ inline double compute_effective_sample_size(std::vector<const double*> draws,
 }
 
 /**
+ * \deprecated use split_rank_normalized_ess instead
+ *
  * Computes the effective sample size (ESS) for the specified
  * parameter across all kept samples.  The value returned is the
  * minimum of ESS and the number_total_draws *
@@ -156,6 +165,11 @@ inline double compute_effective_sample_size(std::vector<const double*> draws,
  * @param size size of chains
  * @return effective sample size for the specified parameter
  */
+#if defined(__GNUC__) || defined(__clang__)
+__attribute__((deprecated))
+#elif defined(_MSC_VER)
+__declspec(deprecated)
+#endif
 inline double compute_effective_sample_size(std::vector<const double*> draws,
                                             size_t size) {
   int num_chains = draws.size();
@@ -164,6 +178,8 @@ inline double compute_effective_sample_size(std::vector<const double*> draws,
 }
 
 /**
+ * \deprecated use split_rank_normalized_ess instead
+ *
  * Computes the split effective sample size (ESS) for the specified
  * parameter across all kept samples.  The value returned is the
  * minimum of ESS and the number_total_draws *
@@ -182,6 +198,11 @@ inline double compute_effective_sample_size(std::vector<const double*> draws,
  * @param sizes stores sizes of chains
  * @return effective sample size for the specified parameter
  */
+#if defined(__GNUC__) || defined(__clang__)
+__attribute__((deprecated))
+#elif defined(_MSC_VER)
+__declspec(deprecated)
+#endif
 inline double compute_split_effective_sample_size(
     std::vector<const double*> draws, std::vector<size_t> sizes) {
   int num_chains = sizes.size();
@@ -199,6 +220,8 @@ inline double compute_split_effective_sample_size(
 }
 
 /**
+ * \deprecated use split_rank_normalized_ess instead
+ *
  * Computes the split effective sample size (ESS) for the specified
  * parameter across all kept samples.  The value returned is the
  * minimum of ESS and the number_total_draws *
@@ -218,6 +241,11 @@ inline double compute_split_effective_sample_size(
  * @param size size of chains
  * @return effective sample size for the specified parameter
  */
+#if defined(__GNUC__) || defined(__clang__)
+__attribute__((deprecated))
+#elif defined(_MSC_VER)
+__declspec(deprecated)
+#endif
 inline double compute_split_effective_sample_size(
     std::vector<const double*> draws, size_t size) {
   int num_chains = draws.size();
diff --git a/src/stan/analyze/mcmc/compute_potential_scale_reduction.hpp b/src/stan/analyze/mcmc/compute_potential_scale_reduction.hpp
index e40253ebb4..ff2363e643 100644
--- a/src/stan/analyze/mcmc/compute_potential_scale_reduction.hpp
+++ b/src/stan/analyze/mcmc/compute_potential_scale_reduction.hpp
@@ -18,6 +18,8 @@ namespace stan {
 namespace analyze {
 
 /**
+ * DEPRECATED - re-implemented as function `rhat`.
+ *
  * Computes the potential scale reduction (Rhat) for the specified
  * parameter across all kept samples.
  *
@@ -32,6 +34,11 @@ namespace analyze {
  * @param sizes stores sizes of chains
  * @return potential scale reduction for the specified parameter
  */
+#if defined(__GNUC__) || defined(__clang__)
+__attribute__((deprecated))
+#elif defined(_MSC_VER)
+__declspec(deprecated)
+#endif
 inline double compute_potential_scale_reduction(
     std::vector<const double*> draws, std::vector<size_t> sizes) {
   int num_chains = sizes.size();
@@ -102,6 +109,8 @@ inline double compute_potential_scale_reduction(
 }
 
 /**
+ * \deprecated use split_rank_normalized_rhat instead
+ *
  * Computes the potential scale reduction (Rhat) for the specified
  * parameter across all kept samples.
  *
@@ -117,6 +126,11 @@ inline double compute_potential_scale_reduction(
  * @param size stores sizes of chains
  * @return potential scale reduction for the specified parameter
  */
+#if defined(__GNUC__) || defined(__clang__)
+__attribute__((deprecated))
+#elif defined(_MSC_VER)
+__declspec(deprecated)
+#endif
 inline double compute_potential_scale_reduction(
     std::vector<const double*> draws, size_t size) {
   int num_chains = draws.size();
@@ -125,6 +139,8 @@ inline double compute_potential_scale_reduction(
 }
 
 /**
+ * \deprecated use split_rank_normalized_rhat instead
+ *
  * Computes the split potential scale reduction (Rhat) for the
  * specified parameter across all kept samples.  When the number of
  * total draws N is odd, the (N+1)/2th draw is ignored.
@@ -140,6 +156,11 @@ inline double compute_potential_scale_reduction(
  * @param sizes stores sizes of chains
  * @return potential scale reduction for the specified parameter
  */
+#if defined(__GNUC__) || defined(__clang__)
+__attribute__((deprecated))
+#elif defined(_MSC_VER)
+__declspec(deprecated)
+#endif
 inline double compute_split_potential_scale_reduction(
     std::vector<const double*> draws, std::vector<size_t> sizes) {
   int num_chains = sizes.size();
@@ -157,6 +178,8 @@ inline double compute_split_potential_scale_reduction(
 }
 
 /**
+ * \deprecated use split_rank_normalized_rhat instead
+ *
  * Computes the split potential scale reduction (Rhat) for the
  * specified parameter across all kept samples.  When the number of
  * total draws N is odd, the (N+1)/2th draw is ignored.
@@ -173,6 +196,11 @@ inline double compute_split_potential_scale_reduction(
  * @param size stores sizes of chains
  * @return potential scale reduction for the specified parameter
  */
+#if defined(__GNUC__) || defined(__clang__)
+__attribute__((deprecated))
+#elif defined(_MSC_VER)
+__declspec(deprecated)
+#endif
 inline double compute_split_potential_scale_reduction(
     std::vector<const double*> draws, size_t size) {
   int num_chains = draws.size();
diff --git a/src/stan/analyze/mcmc/ess.hpp b/src/stan/analyze/mcmc/ess.hpp
index 49070894ba..dedc0f7b50 100644
--- a/src/stan/analyze/mcmc/ess.hpp
+++ b/src/stan/analyze/mcmc/ess.hpp
@@ -15,10 +15,13 @@ namespace analyze {
  * Computes the effective sample size (ESS) for the specified
  * parameter across all chains.  The number of draws per chain must be > 3,
  * and the values across all draws must be finite and not constant.
- * The value returned is the minimum of ESS and (sample_sz * log10(sample_sz).
- * Sample autocovariance is computed using Stan math library implmentation.
  * See https://arxiv.org/abs/1903.08008, section 3.2 for discussion.
  *
+ * Sample autocovariance is computed using the implementation in this namespace
+ * which normalizes lag-k autocorrelation estimators by N instead of (N - k),
+ * yielding biased but more stable estimators as discussed in Geyer (1992); see
+ * https://projecteuclid.org/euclid.ss/1177011137.
+ * 
  * @param chains matrix of draws across all chains
  * @return effective sample size for the specified parameter
  */
diff --git a/src/stan/mcmc/chainset.hpp b/src/stan/mcmc/chainset.hpp
index b1f735d344..9eb9369621 100644
--- a/src/stan/mcmc/chainset.hpp
+++ b/src/stan/mcmc/chainset.hpp
@@ -411,11 +411,6 @@ class chainset {
    */
   double mcse_mean(const int index) const {
     return analyze::mcse_mean(samples(index));
-    // if (num_samples() < 4
-    //     || !stan::analyze::is_finite_and_varies(samples(index)))
-    //   return std::numeric_limits<double>::quiet_NaN();
-    // double ess = analyze::ess(samples(index));
-    // return sd(index) / std::sqrt(ess);
   }
 
   /**
@@ -440,17 +435,6 @@ class chainset {
    */
   double mcse_sd(const int index) const {
     return analyze::mcse_sd(samples(index));
-    // if (num_samples() < 4
-    //     || !stan::analyze::is_finite_and_varies(samples(index)))
-    //   return std::numeric_limits<double>::quiet_NaN();
-    // Eigen::MatrixXd s = samples(index);
-    // Eigen::MatrixXd s2 = s.array().square();
-    // double ess_s = analyze::ess(s);
-    // double ess_s2 = analyze::ess(s2);
-    // double ess_sd = std::min(ess_s, ess_s2);
-    // return sd(index)
-    //        * std::sqrt(stan::math::e() * std::pow(1 - 1 / ess_sd, ess_sd - 1)
-    //                    - 1);
   }
 
   /**
diff --git a/src/test/unit/analyze/mcmc/split_rank_normalized_ess_test.cpp b/src/test/unit/analyze/mcmc/split_rank_normalized_ess_test.cpp
index 5632bc5908..34f34426b5 100644
--- a/src/test/unit/analyze/mcmc/split_rank_normalized_ess_test.cpp
+++ b/src/test/unit/analyze/mcmc/split_rank_normalized_ess_test.cpp
@@ -35,6 +35,7 @@ class RankNormalizedEss : public testing::Test {
 };
 
 TEST_F(RankNormalizedEss, test_bulk_tail_ess) {
+  // computed via R pkg posterior
   double ess_lp_bulk_expect = 1512.7684;
   double ess_lp_tail_expect = 1591.9707;
 
diff --git a/src/test/unit/analyze/mcmc/split_rank_normalized_rhat_test.cpp b/src/test/unit/analyze/mcmc/split_rank_normalized_rhat_test.cpp
index f043b9c96e..c9fd1ddb8d 100644
--- a/src/test/unit/analyze/mcmc/split_rank_normalized_rhat_test.cpp
+++ b/src/test/unit/analyze/mcmc/split_rank_normalized_rhat_test.cpp
@@ -37,15 +37,16 @@ class RankNormalizedRhat : public testing::Test {
 };
 
 TEST_F(RankNormalizedRhat, test_bulk_tail_rhat) {
-  double rhat_lp_expect = 1.0007301;
-  double rhat_theta_expect = 1.0067897;
+  // computed via R pkg posterior
+  double rhat_lp_expect = 1.00073;
+  double rhat_theta_expect = 1.006789;
 
   auto rhat_lp = stan::analyze::split_rank_normalized_rhat(chains_lp);
   auto rhat_theta = stan::analyze::split_rank_normalized_rhat(chains_theta);
 
-  EXPECT_NEAR(rhat_lp_expect, std::max(rhat_lp.first, rhat_lp.second), 0.00001);
+  EXPECT_NEAR(rhat_lp_expect, std::max(rhat_lp.first, rhat_lp.second), 1e-5);
   EXPECT_NEAR(rhat_theta_expect, std::max(rhat_theta.first, rhat_theta.second),
-              0.00001);
+              1e-5);
 }
 
 TEST_F(RankNormalizedRhat, const_fail) {
diff --git a/src/test/unit/mcmc/chainset_test.cpp b/src/test/unit/mcmc/chainset_test.cpp
index 0c59ec1a8e..6d5db4bd19 100644
--- a/src/test/unit/mcmc/chainset_test.cpp
+++ b/src/test/unit/mcmc/chainset_test.cpp
@@ -160,50 +160,50 @@ TEST_F(McmcChains, summary_stats) {
 
   Eigen::MatrixXd theta = bern_chains.samples("theta");
   // default summary statistics - via R pkg posterior
-  double theta_mean_expect = 0.2512974105;
-  EXPECT_NEAR(theta_mean_expect, bern_chains.mean("theta"), 0.00001);
+  double theta_mean_expect = 0.251297;
+  EXPECT_NEAR(theta_mean_expect, bern_chains.mean("theta"), 1e-5);
 
   double theta_median_expect = 0.237476;
-  EXPECT_NEAR(theta_median_expect, bern_chains.median("theta"), 0.00001);
+  EXPECT_NEAR(theta_median_expect, bern_chains.median("theta"), 1e-5);
 
-  double theta_sd_expect = 0.1215466867;
-  EXPECT_NEAR(theta_sd_expect, bern_chains.sd("theta"), 0.00001);
+  double theta_sd_expect = 0.121546;
+  EXPECT_NEAR(theta_sd_expect, bern_chains.sd("theta"), 1e-5);
 
-  double theta_mad_expect = 0.1230906411;
+  double theta_mad_expect = 0.12309;
   EXPECT_NEAR(theta_mad_expect, bern_chains.max_abs_deviation("theta"),
-              0.00001);
+              1e-5);
 
-  double theta_mcse_mean_expect = 0.0032339916;
-  EXPECT_NEAR(theta_mcse_mean_expect, bern_chains.mcse_mean("theta"), 0.0001);
+  double theta_mcse_mean_expect = 0.003234;
+  EXPECT_NEAR(theta_mcse_mean_expect, bern_chains.mcse_mean("theta"), 1e-4);
 
-  double theta_mcse_sd_expect = 0.0021642137;
-  EXPECT_NEAR(theta_mcse_sd_expect, bern_chains.mcse_sd("theta"), 0.0001);
+  double theta_mcse_sd_expect = 0.002164;
+  EXPECT_NEAR(theta_mcse_sd_expect, bern_chains.mcse_sd("theta"), 1e-4);
 
   Eigen::VectorXd probs(6);
   probs << 0.0, 0.01, 0.05, 0.95, 0.99, 1.0;
   Eigen::VectorXd quantiles_expect(6);
-  quantiles_expect << 0.004072430, 0.046281211, 0.07716935, 0.47388505,
-      0.574524110, 0.698401000;
+  quantiles_expect << 0.004072, 0.046281, 0.077169, 0.473885,
+      0.574524, 0.698401;
   Eigen::VectorXd theta_quantiles = bern_chains.quantiles("theta", probs);
   for (size_t i = 0; i < probs.size(); ++i) {
-    EXPECT_NEAR(quantiles_expect(i), theta_quantiles(i), 0.00001);
+    EXPECT_NEAR(quantiles_expect(i), theta_quantiles(i), 1e-5);
   }
 
-  double theta_rhat_expect = 1.0067897;
+  double theta_rhat_expect = 1.00679;
   auto rhat = bern_chains.split_rank_normalized_rhat("theta");
-  EXPECT_NEAR(theta_rhat_expect, std::max(rhat.first, rhat.second), 0.00001);
+  EXPECT_NEAR(theta_rhat_expect, std::max(rhat.first, rhat.second), 1e-5);
 
   double theta_ess_bulk_expect = 1407.5124;
   double theta_ess_tail_expect = 1291.7131;
   auto ess = bern_chains.split_rank_normalized_ess("theta");
-  EXPECT_NEAR(theta_ess_bulk_expect, ess.first, 0.0001);
-  EXPECT_NEAR(theta_ess_tail_expect, ess.second, 0.0001);
+  EXPECT_NEAR(theta_ess_bulk_expect, ess.first,  1e-4);
+  EXPECT_NEAR(theta_ess_tail_expect, ess.second, 1e-4);
 
   // autocorrelation - first 10 lags
   Eigen::VectorXd theta_ac_expect(10);
-  theta_ac_expect << 1.000000000000, 0.422042451075, 0.206832857945,
-      0.083833599168, 0.037326065784, 0.025076266911, 0.020038613922,
-      0.013467409681, 0.004762861453, 0.029494701819;
+  theta_ac_expect << 1.00000, 0.42204, 0.20683,
+      0.08383, 0.037326, 0.02507, 0.02003,
+      0.01347, 0.00476, 0.029495;
   auto theta_ac = bern_chains.autocorrelation(0, "theta");
   for (size_t i = 0; i < 10; ++i) {
     EXPECT_NEAR(theta_ac(i), theta_ac_expect(i), 0.0005);