Use testthat 3e and rely on posterior for ESS

R/effective_sample_sizes.R

@@ -202,72 +202,5 @@ psis_n_eff.matrix <- function(w, r_eff = NULL, ...) {
 #' @return MCMC effective sample size based on RStan's calculation.
 #'
 ess_rfun <- function(sims) {
-  if (is.vector(sims)) dim(sims) <- c(length(sims), 1)
-  chains <- ncol(sims)
-  n_samples <- nrow(sims)
-  acov <- lapply(1:chains, FUN = function(i) posterior::autocovariance(sims[,i]))
-  acov <- do.call(cbind, acov)
-  chain_mean <- colMeans(sims)
-  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
-}
-
-
-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
-  }
+  posterior::ess_basic(sims)
 }