diff --git a/rstan/rstan/DESCRIPTION b/rstan/rstan/DESCRIPTION
index d5ca6f952..41b5324d0 100644
--- a/rstan/rstan/DESCRIPTION
+++ b/rstan/rstan/DESCRIPTION
@@ -43,7 +43,8 @@ Imports:
     loo (>= 2.4.1),
     pkgbuild (>= 1.2.0),
     QuickJSR,
-    ggplot2 (>= 3.3.5)
+    ggplot2 (>= 3.3.5),
+    posterior (>= 1.6.1)
 Depends:
     R (>= 3.4.0),
     StanHeaders (>= 2.36.0)
diff --git a/rstan/rstan/R/check_hmc_diagnostics.R b/rstan/rstan/R/check_hmc_diagnostics.R
index d306990fd..22896a64d 100644
--- a/rstan/rstan/R/check_hmc_diagnostics.R
+++ b/rstan/rstan/R/check_hmc_diagnostics.R
@@ -74,19 +74,19 @@ throw_sampler_warnings <- function(object) {
             call. = FALSE, noBreaks. = TRUE)
   
   sims <- as.array(object)
-  rhat <- apply(sims, MARGIN = 3, FUN = Rhat)
+  rhat <- apply(sims, MARGIN = 3, FUN = posterior::rhat)
   if (any(rhat > 1.05, na.rm = TRUE))
       warning("The largest R-hat is ", round(max(rhat), digits = 2),
             ", indicating chains have not mixed.\n",
             "Running the chains for more iterations may help. See\n",
             "https://mc-stan.org/misc/warnings.html#r-hat", call. = FALSE)
-  bulk_ess <- apply(sims, MARGIN = 3, FUN = ess_bulk)
+  bulk_ess <- apply(sims, MARGIN = 3, FUN = posterior::ess_bulk)
   if (any(bulk_ess < 100 * ncol(sims), na.rm = TRUE))
     warning("Bulk Effective Samples Size (ESS) is too low, ",
             "indicating posterior means and medians may be unreliable.\n",
             "Running the chains for more iterations may help. See\n",
             "https://mc-stan.org/misc/warnings.html#bulk-ess", call. = FALSE)
-  tail_ess <- apply(sims, MARGIN = 3, FUN = ess_tail)
+  tail_ess <- apply(sims, MARGIN = 3, FUN = posterior::ess_tail)
   if (any(tail_ess < 100 * ncol(sims), na.rm = TRUE))
     warning("Tail Effective Samples Size (ESS) is too low, indicating ",
             "posterior variances and tail quantiles may be unreliable.\n",
diff --git a/rstan/rstan/R/monitor.R b/rstan/rstan/R/monitor.R
index 8af411bdb..547748cea 100644
--- a/rstan/rstan/R/monitor.R
+++ b/rstan/rstan/R/monitor.R
@@ -16,260 +16,6 @@
 # along with this program; if not, write to the Free Software
 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 
-fft_next_good_size <- function(N) {
-  # Find the optimal next size for the FFT so that
-  # a minimum number of zeros are padded.
-  if (N <= 2)
-    return(2)
-  while (TRUE) {
-    m <- N
-    while ((m %% 2) == 0) m <- m / 2
-    while ((m %% 3) == 0) m <- m / 3
-    while ((m %% 5) == 0) m <- m / 5
-    if (m <= 1)
-      return(N)
-    N <- N + 1
-  }
-}
-
-#' Autocovariance estimates
-#'
-#' Compute autocovariance estimates for every lag for the specified
-#' input sequence using a fast Fourier transform approach. Estimate
-#' for lag t is scaled by N-t.
-#'
-#' @param y A numeric vector forming a sequence of values.
-#'
-#' @return A numeric vector of autocovariances at every lag (scaled by N-lag).
-autocovariance <- function(y) {
-  N <- length(y)
-  M <- fft_next_good_size(N)
-  Mt2 <- 2 * M
-  yc <- y - mean(y)
-  yc <- c(yc, rep.int(0, Mt2 - N))
-  transform <- fft(yc)
-  ac <- fft(Conj(transform) * transform, inverse = TRUE)
-  # use "biased" estimate as recommended by Geyer (1992)
-  ac <- Re(ac)[1:N] / (N^2 * 2)
-  ac
-}
-
-#' Autocorrelation estimates
-#'
-#' Compute autocorrelation estimates for every lag for the specified
-#' input sequence using a fast Fourier transform approach. Estimate
-#' for lag t is scaled by N-t.
-#'
-#' @param y A numeric vector forming a sequence of values.
-#'
-#' @return A numeric vector of autocorrelations at every lag (scaled by N-lag).
-autocorrelation <- function(y) {
-  ac <- autocovariance(y)
-  ac <- ac / ac[1]
-}
-
-#' Rank normalization
-#'
-#' Compute rank normalization for a numeric array. First replace each
-#' value by its rank. Average rank for ties are used to conserve the
-#' number of unique values of discrete quantities. Second, normalize
-#' ranks via the inverse normal transformation.
-#'
-#' @param x A numeric array of values.
-#'
-#' @return A numeric array of rank normalized values with the same
-#'     size as input.
-z_scale <- function(x) {
-  S <- length(x)
-  r <- rank(x, ties.method = 'average')
-  z <- qnorm((r - 1 / 2) / S)
-  z[is.na(x)] <- NA
-  if (!is.null(dim(x))) {
-    # output should have the input dimension
-    z <- array(z, dim = dim(x), dimnames = dimnames(x))
-  }
-  z
-}
-
-#' Rank uniformization
-#'
-#' Compute rank uniformization for a numeric array. First replace each
-#' value by its rank. Average rank for ties are used to conserve the
-#' number of unique values of discrete quantities. Second, uniformize
-#' ranks to scale \code{[1/(2S), 1-1/(2S)]}, where \code{S} is the the number 
-#' of values.
-#'
-#' @param x A numeric array of values.
-#'
-#' @return A numeric array of rank uniformized values with the same
-#'     size as input.
-#'     
-u_scale <- function(x) {
-  S <- length(x)
-  r <- rank(x, ties.method = 'average')
-  u <- (r - 1 / 2) / S
-  u[is.na(x)] <- NA
-  if (!is.null(dim(x))) {
-    # output should have the input dimension
-    u <- array(u, dim = dim(x), dimnames = dimnames(x))
-  }
-  u
-}
-
-#' Rank values
-#'
-#' Compute ranks for a numeric array. First replace each
-#' value by its rank. Average rank for ties are used to conserve the
-#' number of unique values of discrete quantities. Second, normalize
-#' ranks via the inverse normal transformation.
-#'
-#' @param x A numeric array of values.
-#'
-#' @return A numeric array of ranked values with the same
-#'     size as input.
-#'     
-r_scale <- function(x) {
-  S <- length(x)
-  r <- rank(x, ties.method = 'average')
-  r[is.na(x)] <- NA
-  if (!is.null(dim(x))) {
-    # output should have the input dimension
-    r <- array(r, dim = dim(x), dimnames = dimnames(x))
-  }
-  r
-}
-
-split_chains <- function(sims) {
-  # split Markov chains
-  # Args:
-  #   sims: a 2D array of samples (# iter * # chains)
-  if (is.vector(sims)) {
-    dim(sims) <- c(length(sims), 1)
-  }
-  niter <- dim(sims)[1]
-  if (niter == 1L) return(sims)
-  half <- niter / 2
-  cbind(sims[1:floor(half), ], sims[ceiling(half + 1):niter, ])
-}
-
-is_constant <- function(x, tol = .Machine$double.eps) {
-  abs(max(x) - min(x)) < tol
-}
-
-#' Traditional Rhat convergence diagnostic
-#'
-#' Compute the Rhat convergence diagnostic for a single parameter
-#' For split-Rhat, call this with split chains.
-#'
-#' @param sims A 2D array _without_ warmup samples (# iter * # chains).
-#'
-#' @return A single numeric value for Rhat.
-#' 
-#' @references
-#' Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
-#' Paul-Christian Bürkner (2019). Rank-normalization, folding, and
-#' localization: An improved R-hat for assessing convergence of
-#' MCMC. \emph{arXiv preprint} \code{arXiv:1903.08008}.
-rhat_rfun <- function(sims) {
-  if (anyNA(sims)) {
-    return(NA)
-  }
-  if (any(!is.finite(sims))) {
-    return(NaN)
-  }
-  if (is_constant(sims)) {
-    return(NA)
-  }
-  if (is.vector(sims)) {
-    dim(sims) <- c(length(sims), 1)
-  }
-  chains <- ncol(sims)
-  n_samples <- nrow(sims)
-  chain_mean <- numeric(chains)
-  chain_var <- numeric(chains)
-  for (i in seq_len(chains)) {
-    chain_mean[i] <- mean(sims[, i])
-    chain_var[i] <- var(sims[, i])
-  }
-  var_between <- n_samples * var(chain_mean)
-  var_within <- mean(chain_var)
-  sqrt((var_between / var_within + n_samples - 1) / n_samples)
-}
-
-#' Effective sample size
-#'
-#' Compute the effective sample size estimate for a sample of several chains
-#' for one parameter. For split-ESS, call this with split chains.
-#'
-#' @param sims A 2D array _without_ warmup samples (# iter * # chains).
-#'
-#' @return A single numeric value for the effective sample size.
-#'
-#' @references
-#' Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
-#' Paul-Christian Bürkner (2019). Rank-normalization, folding, and
-#' localization: An improved R-hat for assessing convergence of
-#' MCMC. \emph{arXiv preprint} \code{arXiv:1903.08008}.
-ess_rfun <- function(sims) {
-  if (is.vector(sims)) {
-    dim(sims) <- c(length(sims), 1)
-  }
-  chains <- ncol(sims)
-  n_samples <- nrow(sims)
-  if (n_samples < 3L || should_return_NA(sims)) {
-    return(NA_real_)
-  }
-  acov <- lapply(seq_len(chains), function(i) autocovariance(sims[, i]))
-  acov <- do.call(cbind, acov)
-  chain_mean <- apply(sims, 2, mean)
-  mean_var <- mean(acov[1, ]) * n_samples / (n_samples - 1)
-  var_plus <- mean_var * (n_samples - 1) / n_samples
-  if (chains > 1)
-    var_plus <- var_plus + var(chain_mean)
-
-  # Geyer's initial positive sequence
-  rho_hat_t <- rep.int(0, n_samples)
-  t <- 0
-  rho_hat_even <- 1
-  rho_hat_t[t + 1] <- rho_hat_even
-  rho_hat_odd <- 1 - (mean_var - mean(acov[t + 2, ])) / var_plus
-  rho_hat_t[t + 2] <- rho_hat_odd
-  while (t < nrow(acov) - 5 && !is.nan(rho_hat_even + rho_hat_odd) &&
-         (rho_hat_even + rho_hat_odd > 0)) {
-    t <- t + 2
-    rho_hat_even = 1 - (mean_var - mean(acov[t + 1, ])) / var_plus
-    rho_hat_odd = 1 - (mean_var - mean(acov[t + 2, ])) / var_plus
-    if ((rho_hat_even + rho_hat_odd) >= 0) {
-      rho_hat_t[t + 1] <- rho_hat_even
-      rho_hat_t[t + 2] <- rho_hat_odd
-    }
-  }
-  max_t <- t
-  # this is used in the improved estimate
-  if (rho_hat_even>0)
-      rho_hat_t[max_t + 1] <- rho_hat_even
-  
-  # Geyer's initial monotone sequence
-  t <- 0
-  while (t <= max_t - 4) {
-    t <- t + 2
-    if (rho_hat_t[t + 1] + rho_hat_t[t + 2] >
-        rho_hat_t[t - 1] + rho_hat_t[t]) {
-      rho_hat_t[t + 1] = (rho_hat_t[t - 1] + rho_hat_t[t]) / 2;
-      rho_hat_t[t + 2] = rho_hat_t[t + 1];
-    }
-  }
-  ess <- chains * n_samples
-  # Geyer's truncated estimate
-  # tau_hat <- -1 + 2 * sum(rho_hat_t[1:max_t])
-  # Improved estimate reduces variance in antithetic case
-  tau_hat <- -1 + 2 * sum(rho_hat_t[1:max_t]) + rho_hat_t[max_t+1]
-  # Safety check for negative values and with max ess equal to ess*log10(ess)
-  tau_hat <- max(tau_hat, 1/log10(ess))
-  ess <- ess / tau_hat
-  ess
-}
-
 #' Rhat convergence diagnostic
 #'
 #' Compute Rhat convergence diagnostic as the maximum of rank normalized
@@ -287,10 +33,7 @@ ess_rfun <- function(sims) {
 #' 
 #' @export
 Rhat <- function(sims) {
-  bulk_rhat <- rhat_rfun(z_scale(split_chains(sims)))
-  sims_folded <- abs(sims - median(sims))
-  tail_rhat <- rhat_rfun(z_scale(split_chains(sims_folded)))
-  max(bulk_rhat, tail_rhat)
+  posterior::rhat(sims)
 }
 
 #' Bulk effective sample size (bulk-ESS)
@@ -312,7 +55,7 @@ Rhat <- function(sims) {
 #' 
 #' @export
 ess_bulk <- function(sims) {
-  ess_rfun(z_scale(split_chains(sims)))
+  posterior::ess_bulk(sims)
 }
 
 #' Tail effective sample size (tail-ESS)
@@ -334,11 +77,7 @@ ess_bulk <- function(sims) {
 #' 
 #' @export
 ess_tail <- function(sims) {
-  I05 <- sims <= quantile(sims, 0.05, na.rm = TRUE)
-  q05_ess <- ess_rfun(split_chains(I05))
-  I95 <- sims <= quantile(sims, 0.95, na.rm = TRUE)
-  q95_ess <- ess_rfun(split_chains(I95))
-  min(q05_ess, q95_ess)
+  posterior::ess_tail(sims)
 }
 
 #' Quantile effective sample size
@@ -360,11 +99,7 @@ ess_tail <- function(sims) {
 #' 
 #' @export
 ess_quantile <- function(sims, prob) {
-  if (should_return_NA(sims)) {
-    return(NA_real_)
-  }
-  I <- sims <= quantile(sims, prob, na.rm = TRUE)
-  ess_rfun(split_chains(I))
+  posterior::ess_quantile(sims, prob)
 }
 
 #' Effective sample size
@@ -385,7 +120,7 @@ ess_quantile <- function(sims, prob) {
 #' 
 #' @export
 ess_mean <- function(sims) {
-  ess_rfun(split_chains(sims))
+  posterior::ess_mean(sims)
 }
 
 #' Effective sample size
@@ -407,43 +142,11 @@ ess_mean <- function(sims) {
 #' 
 #' @export
 ess_sd <- function(sims) {
-  min(ess_rfun(split_chains(sims)), ess_rfun(split_chains(sims^2)))
+  # TODO: are these two equivalent/ok to change to posterior's implementation?
+  # min(ess_rfun(split_chains(sims)), ess_rfun(split_chains(sims^2)))
+  posterior::ess_sd(sims)
 }
 
-#' Monte Carlo diagnostics for a quantile
-#'
-#' Compute Monte Carlo standard error, 5%-quantile, 95%-quantile, and
-#' effective sample size estimate for a quantile estimate of a single
-#' parameter.
-#'
-#' @param sims A 2D array _without_ warmup samples (# iter * # chains).
-#' @param prob A single numeric value of probability.
-#'
-#' @return A data frame with Monte Carlo standard error (mcse),
-#'     5%-quantile (Q05), 95%-quantile (Q95), and effective sample
-#'     size estimate (ess).
-#'
-#' @references
-#' Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and
-#' Paul-Christian Bürkner (2019). Rank-normalization, folding, and
-#' localization: An improved R-hat for assessing convergence of
-#' MCMC. \emph{arXiv preprint} \code{arXiv:1903.08008}.
-conv_quantile <- function(sims, prob) {
-  if (is.vector(sims)) {
-    dim(sims) <- c(length(sims), 1)
-  }
-  ess <- ess_quantile(sims, prob)
-  p <- c(0.1586553, 0.8413447, 0.05, 0.95)
-  a <- qbeta(p, ess * prob + 1, ess * (1 - prob) + 1)
-  ssims <- sort(sims)
-  S <- length(ssims)
-  th1 <- ssims[max(round(a[1] * S), 1)]
-  th2 <- ssims[min(round(a[2] * S), S)]
-  mcse <- (th2 - th1) / 2
-  th1 <- ssims[max(round(a[3] * S), 1)]
-  th2 <- ssims[min(round(a[4] * S), S)]
-  data.frame(mcse = mcse, Q05 = th1, Q95 = th2, ess = ess)
-}
 
 #' Monte Carlo standard error for a quantile
 #'
@@ -464,7 +167,7 @@ conv_quantile <- function(sims, prob) {
 #' 
 #' @export
 mcse_quantile <- function(sims, prob) {
-  conv_quantile(sims, prob)$mcse
+  posterior::mcse_quantile(sims, prob)
 }
 
 #' Monte Carlo standard error for mean
@@ -513,9 +216,7 @@ mcse_mean <- function(sims) {
 #' 
 #' @export
 mcse_sd <- function(sims) {
-  # assumes normality of sims and uses Stirling's approximation
-  ess_sd <- ess_sd(sims)
-  sd(sims) * sqrt(exp(1) * (1 - 1 / ess_sd)^(ess_sd - 1) - 1)
+  posterior::mcse_sd(sims)
 }
 
 #' Summary of General Simulation Results
@@ -548,9 +249,9 @@ mcse_sd <- function(sims) {
 #' MCMC. \emph{arXiv preprint} \code{arXiv:1903.08008}.
 #' 
 #' @export
-monitor <- function(sims, warmup = floor(dim(sims)[1] / 2), 
-                    probs = c(0.025, 0.25, 0.50, 0.75, 0.975), 
-                    digits_summary = 1, print = TRUE, ...) { 
+monitor <- function(sims, warmup = floor(dim(sims)[1] / 2),
+                    probs = c(0.025, 0.25, 0.50, 0.75, 0.975),
+                    digits_summary = 1, print = TRUE, ...) {
   if (inherits(sims, "stanfit")) {
     chains <- sims@sim$chains
     iter <- sims@sim$iter
@@ -586,30 +287,32 @@ monitor <- function(sims, warmup = floor(dim(sims)[1] / 2),
   for (i in seq_along(out)) {
     sims_i <- sims[, , i]
     valid <- all(is.finite(sims_i))
-    quan <- unname(quantile(sims_i, probs = probs, na.rm = TRUE))
-    quan2 <- quantile(sims_i, probs = c(0.05, 0.5, 0.95), na.rm = TRUE)
+    quan <- unname(posterior::quantile2(sims_i, probs = probs))
+    #quan2 <- posterior::quantile2(sims_i, probs = c(0.05, 0.5, 0.95))
     mean <- mean(sims_i)
     sd <- sd(sims_i)
-    mcse_quan <- sapply(probs, mcse_quantile, sims = sims_i)
-    mcse_mean <- mcse_mean(sims_i)
-    mcse_sd <- mcse_sd(sims_i)
-    rhat <- Rhat(sims_i)
-    ess_bulk <- round(ess_bulk(sims_i))
-    ess_tail <- round(ess_tail(sims_i))
-    ess <- round(ess_rfun(sims_i))
+    mcse_quan <- sapply(probs, function(p) posterior::mcse_quantile(sims_i, probs = p))
+    mcse_mean <- posterior::mcse_mean(sims_i)
+    mcse_sd <- posterior::mcse_sd(sims_i)
+    rhat <- posterior::rhat(sims_i)
+    ess_bulk <- round(posterior::ess_bulk(sims_i))
+    ess_tail <- round(posterior::ess_tail(sims_i))
+    ess <- round(posterior::ess_bulk(sims_i))
     out[[i]] <- c(
       mean, mcse_mean, sd, quan, ess, rhat,
-      valid, quan2, mcse_quan, mcse_sd, ess_bulk, ess_tail
+      valid, #quan2,
+      mcse_quan, mcse_sd, ess_bulk, ess_tail
     )
   }
   
   out <- as.data.frame(do.call(rbind, out))
-  probs_str <- names(quantile(sims_i, probs = probs, na.rm = TRUE))
+  #probs_str <- names(quantile(sims_i, probs = probs, na.rm = TRUE))
   str_quan <- paste0("Q", probs * 100)
-  str_quan2 <- paste0("Q", c(0.05, 0.5, 0.95) * 100)
+  #str_quan2 <- paste0("Q", c(0.05, 0.5, 0.95) * 100)
   str_mcse_quan <- paste0("MCSE_", str_quan)
-  colnames(out) <- c("mean", "se_mean", "sd", probs_str, "n_eff", "Rhat",
-                     "valid", str_quan2, str_mcse_quan, "MCSE_SD", "Bulk_ESS", "Tail_ESS")
+  colnames(out) <- c("mean", "se_mean", "sd", str_quan, "n_eff", "Rhat",
+                     "valid", #str_quan2,
+                     str_mcse_quan, "MCSE_SD", "Bulk_ESS", "Tail_ESS")
   rownames(out) <- parnames
 
   # replace NAs with appropriate values if draws are valid
diff --git a/rstan/rstan/tests/testthat/_snaps/monitor.md b/rstan/rstan/tests/testthat/_snaps/monitor.md
new file mode 100644
index 000000000..a44a93222
--- /dev/null
+++ b/rstan/rstan/tests/testthat/_snaps/monitor.md
@@ -0,0 +1,240 @@
+# rstan monitor snapshot is stable
+
+    {
+      "type": "list",
+      "attributes": {
+        "names": {
+          "type": "character",
+          "attributes": {},
+          "value": ["mean", "se_mean", "sd", "Q2.5", "Q25", "Q50", "Q75", "Q97.5", "n_eff", "Rhat", "valid", "MCSE_Q2.5", "MCSE_Q25", "MCSE_Q50", "MCSE_Q75", "MCSE_Q97.5", "MCSE_SD", "Bulk_ESS", "Tail_ESS"]
+        },
+        "row.names": {
+          "type": "character",
+          "attributes": {},
+          "value": ["y[1]", "y[2]", "lp__"]
+        },
+        "chains": {
+          "type": "integer",
+          "attributes": {},
+          "value": [4]
+        },
+        "iter": {
+          "type": "integer",
+          "attributes": {},
+          "value": [5]
+        },
+        "warmup": {
+          "type": "double",
+          "attributes": {},
+          "value": [2]
+        },
+        "class": {
+          "type": "character",
+          "attributes": {},
+          "value": ["data.frame"]
+        }
+      },
+      "value": [
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [-0.053848, 0.092345, -0.667653]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.397088, 0.503423, 0.227732]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.794177, 1.006846, 0.455463]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [-1.203275, -1.125077, -1.571521]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [-0.623548, -0.473779, -0.880547]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [-0.049542, -0.239292, -0.615925]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.488907, 0.478004, -0.454857]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [1.135881, 2.017441, -0.086398]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [4, 4, 4]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.992545, 0.939105, 0.992545]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [1, 1, 1]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.32821, "NA", 0.469387]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.590813, 0.491098, 0.610703]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.590063, 0.503575, 0.197583]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.571268, 1.050291, 0.28849]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.537518, 0.982866, 0.235722]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.184646, 0.321772, 0.206564]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [4, 4, 4]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [4, 12, 4]
+        }
+      ]
+    }
+
+---
+
+    {
+      "type": "list",
+      "attributes": {
+        "names": {
+          "type": "character",
+          "attributes": {},
+          "value": ["mean", "se_mean", "sd", "Q10", "Q90", "n_eff", "Rhat", "valid", "MCSE_Q10", "MCSE_Q90", "MCSE_SD", "Bulk_ESS", "Tail_ESS"]
+        },
+        "row.names": {
+          "type": "character",
+          "attributes": {},
+          "value": ["y[1]", "y[2]", "lp__"]
+        },
+        "chains": {
+          "type": "integer",
+          "attributes": {},
+          "value": [4]
+        },
+        "iter": {
+          "type": "integer",
+          "attributes": {},
+          "value": [5]
+        },
+        "warmup": {
+          "type": "double",
+          "attributes": {},
+          "value": [2]
+        },
+        "class": {
+          "type": "character",
+          "attributes": {},
+          "value": ["data.frame"]
+        }
+      },
+      "value": [
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [-0.053848, 0.092345, -0.667653]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.397088, 0.503423, 0.227732]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.794177, 1.006846, 0.455463]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [-1.11838, -0.911257, -0.914781]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.867839, 1.651193, -0.142213]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [4, 4, 4]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.992545, 0.939105, 0.992545]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [1, 1, 1]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.476953, 0.428398, 0.562226]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.63768, 1.161523, 0.261648]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [0.184646, 0.321772, 0.206564]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [4, 4, 4]
+        },
+        {
+          "type": "double",
+          "attributes": {},
+          "value": [4, 12, 4]
+        }
+      ]
+    }
+
diff --git a/rstan/rstan/tests/testthat/test-monitor.R b/rstan/rstan/tests/testthat/test-monitor.R
new file mode 100644
index 000000000..800872117
--- /dev/null
+++ b/rstan/rstan/tests/testthat/test-monitor.R
@@ -0,0 +1,30 @@
+test_that("rstan monitor snapshot is stable", {
+  x <- structure(c(-0.665824869609518, 1.58830022066106, -0.565535490555016,
+              -0.797584842767636, 1.23502996632411, 0.772693257183606, -1.42729906097828,
+              -0.268049350062425, 0.874487281622664, 0.38254123399846, 0.245743405430849,
+              -0.416200361337001, 0.808003109795885, -0.0587539995594358, 0.159993102962675,
+              1.00945624491801, -0.157315599905492, -0.0403299275412378, -1.15402427067623,
+              -1.22195574738444, 1.41275401125985, 0.568417727166425, -0.948278118922285,
+              -0.33534143563398, 2.11310307328816, 1.42064239670371, 0.222846831903693,
+              1.76523990347241, 0.624771161221411, -0.578067144520362, -1.1061183309293,
+              -1.13675871572812, -0.439016877278619, -1.19213836974204, -0.209942138956473,
+              -0.22163939207362, 0.0688433795044292, 0.147371145111933, 0.429082133577958,
+              -0.268640866689185, -0.928038384125191, -1.54555765905919, -0.634054254999794,
+              -0.485741508523327, -1.81920104550335, -1.00884863320036, -1.13001472068659,
+              -0.918545178770649, -0.694749583470604, -0.362202470114714, -0.583254076121026,
+              -0.654990728252587, -0.545942951359219, -0.597795201103136, -0.117769965976049,
+              -0.620320651238796, -0.046795788739027, -0.0744988240837071,
+              -0.880427075443884, -0.880908357627523), dim = 5:3, dimnames = list(
+                iterations = NULL, chains = c("chain:1", "chain:2", "chain:3",
+                                              "chain:4"), parameters = c("y[1]", "y[2]", "lp__")))
+
+  expect_snapshot_value(
+    monitor(x) |> as.data.frame() |> round(6),
+    style = "json2"
+  )
+
+  expect_snapshot_value(
+    monitor(x, probs = c(0.1, 0.9)) |> as.data.frame() |> round(6),
+    style = "json2"
+  )
+})