diff --git a/DESCRIPTION b/DESCRIPTION
index 916a81f..93526b8 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: poco
 Type: Package
 Title: poco (POsterior COmpression): Compress Bayesian Posteriors Using Mixture and Density Approximations
-Version: 0.1.0
+Version: 0.2.0
 Authors@R: c(
     person("Chen", "Zhan", , "chen.zhan@adelaide.edu.au",
            role = c("aut", "cre")),
@@ -26,7 +26,8 @@ Imports:
     utils,
     methods,
     cli,
-    rlang
+    rlang,
+    ranger
 Suggests:
     posterior,
     cmdstanr,
diff --git a/NAMESPACE b/NAMESPACE
index e2be97a..5f930ce 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+S3method(print,compression_fidelity)
 S3method(print,posterior_compressed)
 S3method(print,summary.posterior_compressed)
 S3method(summary,posterior_compressed)
@@ -11,6 +12,7 @@ export(compress_posterior)
 export(compress_sccomp)
 export(compression_methods)
 export(density_posterior)
+export(evaluate_compression)
 export(get_n_draws_from_fit)
 export(has_cmdstan)
 export(poco_stan_models_cache_dir)
@@ -24,6 +26,7 @@ importFrom(methods,is)
 importFrom(mvtnorm,dmvnorm)
 importFrom(mvtnorm,rmvnorm)
 importFrom(stats,cov)
+importFrom(stats,predict)
 importFrom(stats,sd)
 importFrom(utils,head)
 importFrom(utils,object.size)
diff --git a/R/evaluate.R b/R/evaluate.R
new file mode 100644
index 0000000..f481750
--- /dev/null
+++ b/R/evaluate.R
@@ -0,0 +1,386 @@
+#' Evaluate the fidelity of a compressed posterior
+#'
+#' Compares draws regenerated from a compressed posterior against a
+#' reference set of draws (typically the original MCMC output) using
+#' two distribution-free, backend-independent two-sample diagnostics:
+#'
+#' \itemize{
+#'   \item \strong{Energy distance} (Szekely & Rizzo, 2013): a
+#'         scale-invariant distance between two samples that is zero
+#'         iff their distributions match. The raw distance is anchored
+#'         against a \emph{self-baseline} obtained by bootstrapping
+#'         `reference_draws` to two independent samples matched in
+#'         size to the reconstructed sample. The 90th percentile of
+#'         this baseline distribution defines the Monte Carlo noise
+#'         envelope; the reproduction score is
+#'         `100 * min(1, noise_envelope / max(distance, noise_envelope))`,
+#'         so 100\% means the reconstruction is within the bootstrap
+#'         envelope and the score decays as 1/ratio for larger
+#'         distances.
+#'   \item \strong{Classifier two-sample test (C2ST)} (Lopez-Paz &
+#'         Oquab, 2017): trains a classifier (random forest via
+#'         `ranger::ranger()` under cross-validation
+#'         to discriminate reference draws from reconstructed draws and
+#'         reports the out-of-fold ROC AUC. AUC = 0.5 means the two
+#'         samples are indistinguishable; the score is mapped to
+#'         `100 * (1 - 2 * |AUC - 0.5|)`.
+#' }
+#'
+#' Both metrics are computed only from samples and are independent of
+#' the compression backend (`mclust`, `mvdens_gmm`, `mvdens_kde`), so
+#' the score is not "circular": it does not depend on the same density
+#' model that produced the compression.
+#'
+#' For a strict out-of-sample evaluation, hold out a fraction of the
+#' draws \emph{before} fitting:
+#'
+#' \preformatted{
+#'   idx  <- sample.int(nrow(draws), size = 0.8 * nrow(draws))
+#'   comp <- compress_posterior(draws[idx, ], method = "mclust")
+#'   evaluate_compression(comp, reference_draws = draws[-idx, ])
+#' }
+#'
+#' @param comp A `posterior_compressed` object (or a path to an `.rds`
+#'   file containing one).
+#' @param reference_draws A draws matrix, data.frame, or
+#'   `posterior::draws_*` object whose columns include all parameters
+#'   in `comp`. Treated as the ground-truth distribution.
+#' @param metric Character vector of metrics to compute. One or more of
+#'   `"energy"` and `"c2st"`. Default: both.
+#' @param n_eval Integer number of draws to regenerate from `comp`.
+#'   Defaults to `min(max_n, nrow(reference_draws))`.
+#' @param max_n Cap on the number of points used in any pairwise
+#'   distance / classifier computation. Both samples are subsampled to
+#'   at most this many rows. Default `2000`.
+#' @param n_self_reps Number of self-baseline replicates for the energy
+#'   metric. Default `20`.
+#' @param classifier Classifier for C2ST. `"ranger"` (default) uses a
+#'   random forest (`ranger::ranger()`). `"knn"` uses a k-NN probability
+#'   estimate instead (no random forest).
+#' @param cv_folds Cross-validation folds for C2ST. Default `5`.
+#' @param seed Optional integer seed for reproducibility.
+#' @param verbose Logical; print progress messages.
+#'
+#' @return An S3 object of class `compression_fidelity` containing:
+#' \describe{
+#'   \item{`reproduction_pct`}{Headline reproduction score in `[0, 100]`,
+#'     averaged across requested metrics.}
+#'   \item{`metrics`}{Named list with detailed per-metric results.}
+#'   \item{`n_reference`, `n_eval`, `n_params`}{Sample sizes used.}
+#' }
+#'
+#' @references
+#' Szekely, G. J. & Rizzo, M. L. (2013). Energy statistics: A class of
+#' statistics based on distances. *Journal of Statistical Planning and
+#' Inference*, 143(8), 1249-1272.
+#'
+#' Lopez-Paz, D. & Oquab, M. (2017). Revisiting Classifier Two-Sample
+#' Tests. *ICLR*.
+#'
+#' @examples
+#' set.seed(1)
+#' draws <- matrix(rnorm(2000 * 3), ncol = 3,
+#'                 dimnames = list(NULL, c("alpha", "beta", "sigma")))
+#' comp <- compress_posterior(draws, method = "mclust", n_components = 2)
+#' fidelity <- evaluate_compression(comp, reference_draws = draws, seed = 1L)
+#' fidelity
+#'
+#' @importFrom stats predict
+#'
+#' @export
+evaluate_compression <- function(
+    comp,
+    reference_draws,
+    metric       = c("energy", "c2st"),
+    n_eval       = NULL,
+    max_n        = 2000L,
+    n_self_reps  = 20L,
+    classifier   = c("ranger", "knn"),
+    cv_folds     = 5L,
+    seed         = NULL,
+    verbose      = FALSE) {
+  comp <- .resolve_compressed(comp)
+  if (!inherits(comp, "posterior_compressed")) {
+    stop(
+      "`comp` must be a posterior_compressed object (or a path to one).",
+      call. = FALSE
+    )
+  }
+
+  metric     <- match.arg(metric, choices = c("energy", "c2st"),
+                          several.ok = TRUE)
+  classifier <- match.arg(classifier)
+
+  ref_mat <- .as_draws_matrix(reference_draws, variables = comp$param_names)
+  if (ncol(ref_mat) == 0L) {
+    stop("`reference_draws` has no columns matching the compressed posterior.",
+         call. = FALSE)
+  }
+  if (nrow(ref_mat) < 4L) {
+    stop("`reference_draws` must contain at least 4 rows for evaluation.",
+         call. = FALSE)
+  }
+
+  if (!is.null(seed)) set.seed(seed)
+
+  if (is.null(n_eval)) n_eval <- min(max_n, nrow(ref_mat))
+  if (verbose) {
+    .message("Generating ", n_eval, " draws from compressed posterior ...")
+  }
+  recon_mat <- sample_posterior(comp, n_draws = n_eval)
+  recon_mat <- recon_mat[, comp$param_names, drop = FALSE]
+
+  ref_used   <- .subsample_rows(ref_mat,   max_n)
+  recon_used <- .subsample_rows(recon_mat, max_n)
+
+  scaling      <- .scaling_from(ref_used)
+  ref_scaled   <- .apply_scaling(ref_used,   scaling)
+  recon_scaled <- .apply_scaling(recon_used, scaling)
+
+  metrics_out <- list()
+
+  if ("energy" %in% metric) {
+    if (verbose) .message("Computing energy distance ...")
+    metrics_out$energy <- .energy_metric(
+      ref_scaled, recon_scaled, n_self_reps = n_self_reps
+    )
+  }
+  if ("c2st" %in% metric) {
+    if (verbose) .message("Computing classifier two-sample test ...")
+    metrics_out$c2st <- .c2st_metric(
+      ref_scaled, recon_scaled,
+      classifier = classifier,
+      cv_folds   = cv_folds
+    )
+  }
+
+  pcts <- vapply(metrics_out, function(m) m$reproduction_pct, numeric(1))
+  reproduction_pct <- if (length(pcts) > 0L) mean(pcts) else NA_real_
+
+  out <- list(
+    reproduction_pct = reproduction_pct,
+    metrics          = metrics_out,
+    n_reference      = nrow(ref_used),
+    n_eval           = nrow(recon_used),
+    n_params         = ncol(ref_used),
+    method           = comp$method,
+    param_names      = comp$param_names
+  )
+  class(out) <- c("compression_fidelity", "list")
+  out
+}
+
+
+#' @rdname evaluate_compression
+#' @param x A `compression_fidelity` object.
+#' @param ... Currently unused.
+#' @export
+print.compression_fidelity <- function(x, ...) {
+  cat("<compression_fidelity>\n")
+  cat(sprintf(
+    "  method        : %s\n  parameters    : %d\n  reference n   : %d\n  eval n        : %d\n",
+    x$method, x$n_params, x$n_reference, x$n_eval
+  ))
+  cat("  ----------------------------------------\n")
+  if (!is.null(x$metrics$energy)) {
+    e <- x$metrics$energy
+    cat(sprintf(
+      "  energy        : %5.1f%% reproduction\n",
+      e$reproduction_pct
+    ))
+    cat(sprintf(
+      "    distance    : %.4f   noise envelope (q90): %.4f   ratio: %.2fx\n",
+      e$energy, e$noise_envelope, e$ratio
+    ))
+  }
+  if (!is.null(x$metrics$c2st)) {
+    c2 <- x$metrics$c2st
+    cat(sprintf(
+      "  C2ST          : %5.1f%% reproduction\n",
+      c2$reproduction_pct
+    ))
+    cat(sprintf(
+      "    AUC         : %.3f   classifier: %s   cv_folds: %d\n",
+      c2$auc, c2$classifier, c2$cv_folds
+    ))
+  }
+  cat("  ----------------------------------------\n")
+  cat(sprintf("  reproduction  : %.1f%%\n", x$reproduction_pct))
+  invisible(x)
+}
+
+
+# -- internals ---------------------------------------------------------------
+
+#' @keywords internal
+#' @noRd
+.subsample_rows <- function(x, max_n) {
+  if (nrow(x) <= max_n) return(x)
+  idx <- sample.int(nrow(x), size = max_n, replace = FALSE)
+  x[idx, , drop = FALSE]
+}
+
+#' @keywords internal
+#' @noRd
+.scaling_from <- function(x) {
+  mu <- colMeans(x)
+  sd <- apply(x, 2, stats::sd)
+  sd[!is.finite(sd) | sd == 0] <- 1
+  list(center = mu, scale = sd)
+}
+
+#' @keywords internal
+#' @noRd
+.apply_scaling <- function(x, scaling) {
+  sweep(sweep(x, 2, scaling$center, "-"), 2, scaling$scale, "/")
+}
+
+
+# -- energy distance ---------------------------------------------------------
+
+#' @keywords internal
+#' @noRd
+.energy_metric <- function(ref, recon, n_self_reps = 20L) {
+  ed_test <- .energy_distance(ref, recon)
+
+  m <- nrow(recon)
+  n <- nrow(ref)
+  ed_self_vals <- replicate(n_self_reps, {
+    idx_a <- sample.int(n, size = m, replace = TRUE)
+    idx_b <- sample.int(n, size = m, replace = TRUE)
+    .energy_distance(
+      ref[idx_a, , drop = FALSE],
+      ref[idx_b, , drop = FALSE]
+    )
+  })
+  ed_self_mean <- mean(ed_self_vals)
+  ed_self_q90  <- as.numeric(stats::quantile(ed_self_vals, 0.9, names = FALSE))
+
+  envelope <- max(ed_self_q90, .Machine$double.eps)
+  ratio    <- ed_test / envelope
+  pct      <- 100 * min(1, 1 / max(1, ratio))
+
+  list(
+    reproduction_pct = pct,
+    energy           = ed_test,
+    self_baseline    = ed_self_mean,
+    noise_envelope   = ed_self_q90,
+    ratio            = ratio,
+    n_self_reps      = n_self_reps
+  )
+}
+
+#' @keywords internal
+#' @noRd
+.energy_distance <- function(x, y) {
+  d_xy <- .mean_pdist_cross(x, y)
+  d_xx <- .mean_pdist_within(x)
+  d_yy <- .mean_pdist_within(y)
+  ed2 <- 2 * d_xy - d_xx - d_yy
+  sqrt(max(0, ed2))
+}
+
+#' @keywords internal
+#' @noRd
+.mean_pdist_within <- function(x) {
+  n <- nrow(x)
+  if (n < 2L) return(0)
+  d <- as.numeric(stats::dist(x))
+  (2 / (n^2)) * sum(d)
+}
+
+#' @keywords internal
+#' @noRd
+.mean_pdist_cross <- function(x, y) {
+  xx <- rowSums(x^2)
+  yy <- rowSums(y^2)
+  d2 <- outer(xx, yy, "+") - 2 * tcrossprod(x, y)
+  d2[d2 < 0] <- 0
+  mean(sqrt(d2))
+}
+
+# -- C2ST --------------------------------------------------------------------
+
+#' @keywords internal
+#' @noRd
+.c2st_metric <- function(ref, recon,
+                         classifier = "ranger",
+                         cv_folds   = 5L) {
+  Z <- rbind(ref, recon)
+  labels <- c(rep(0L, nrow(ref)), rep(1L, nrow(recon)))
+
+  n <- nrow(Z)
+  cv_folds <- max(2L, min(as.integer(cv_folds), n - 1L))
+  folds <- sample(rep(seq_len(cv_folds), length.out = n))
+
+  preds <- numeric(n)
+  for (k in seq_len(cv_folds)) {
+    test_idx  <- which(folds == k)
+    train_idx <- which(folds != k)
+    if (length(test_idx) == 0L || length(train_idx) == 0L) next
+    preds[test_idx] <- .classifier_predict(
+      classifier,
+      Z[train_idx, , drop = FALSE],
+      labels[train_idx],
+      Z[test_idx, , drop = FALSE]
+    )
+  }
+
+  auc <- .auc(labels, preds)
+  pct <- 100 * max(0, 1 - 2 * abs(auc - 0.5))
+
+  list(
+    reproduction_pct = pct,
+    auc              = auc,
+    classifier       = classifier,
+    cv_folds         = cv_folds
+  )
+}
+
+#' @keywords internal
+#' @noRd
+.classifier_predict <- function(classifier, x_train, y_train, x_test) {
+  if (classifier == "ranger") {
+    df <- as.data.frame(x_train)
+    df$.label <- factor(y_train, levels = c(0L, 1L))
+    mod <- ranger::ranger(
+      .label ~ .,
+      data        = df,
+      probability = TRUE,
+      num.trees   = 200L,
+      verbose     = FALSE
+    )
+    df_test <- as.data.frame(x_test)
+    pr <- predict(mod, df_test)$predictions
+    cn <- colnames(pr)
+    pr[, if ("1" %in% cn) "1" else cn[length(cn)]]
+  } else {
+    k <- max(3L, min(50L, floor(sqrt(nrow(x_train)))))
+    .knn_proba(x_train, y_train, x_test, k = k)
+  }
+}
+
+#' @keywords internal
+#' @noRd
+.knn_proba <- function(x_train, y_train, x_test, k) {
+  xx <- rowSums(x_test^2)
+  yy <- rowSums(x_train^2)
+  d2 <- outer(xx, yy, "+") - 2 * tcrossprod(x_test, x_train)
+  d2[d2 < 0] <- 0
+  apply(d2, 1, function(d_row) {
+    nn <- order(d_row)[seq_len(min(k, length(d_row)))]
+    mean(y_train[nn])
+  })
+}
+
+#' @keywords internal
+#' @noRd
+.auc <- function(labels, scores) {
+  pos <- scores[labels == 1L]
+  neg <- scores[labels == 0L]
+  if (length(pos) == 0L || length(neg) == 0L) return(0.5)
+  r <- rank(c(pos, neg), ties.method = "average")
+  r_pos <- r[seq_along(pos)]
+  (sum(r_pos) - length(pos) * (length(pos) + 1) / 2) /
+    (length(pos) * length(neg))
+}
diff --git a/README.md b/README.md
index 30c1fad..bad00e7 100644
--- a/README.md
+++ b/README.md
@@ -51,6 +51,8 @@ consumed by:
 
 - `sample_posterior()` - regenerate draws (any number),
 - `density_posterior()` - evaluate p.d.f. (or log p.d.f.) at given points,
+- `evaluate_compression()` - quantify how much information was lost
+  (a "98.3% reproduction"-style score, see below),
 - `reconstruct_brmsfit()` / `reconstruct_sccomp()` - rebuild a usable
   fit object so the standard post-processing keeps working.
 
@@ -117,6 +119,48 @@ new_draws <- sample_posterior(comp, n_draws = 4000)
 
 `compression_methods()` returns the current set.
 
+## How much information was lost?
+
+Like JPEG quality scores, `evaluate_compression()` reports a
+`reproduction %` for the compressed posterior. It uses two
+*backend-independent* two-sample diagnostics on samples regenerated
+from the compressed object versus the reference draws — so the score
+is not circular (it doesn't reuse the same density model that did the
+compression):
+
+- **Energy distance** (Szekely & Rizzo) anchored against a bootstrap
+  Monte Carlo noise envelope from the reference itself.
+- **Classifier two-sample test (C2ST)** — a random forest via
+  `ranger::ranger()` under cross-validation tries to tell reference
+  draws from reconstructed draws; AUC of 0.5 means indistinguishable.
+  Pass `classifier = "knn"` for a k-NN alternative instead of the forest.
+
+```r
+fidelity <- evaluate_compression(comp, reference_draws = draws)
+fidelity
+#> <compression_fidelity>
+#>   method        : mclust
+#>   parameters    : 3
+#>   reference n   : 1500
+#>   eval n        : 1500
+#>   ----------------------------------------
+#>   energy        : 100.0% reproduction
+#>     distance    : 0.054   noise envelope (q90): 0.062   ratio: 0.87x
+#>   C2ST          :  97.7% reproduction
+#>     AUC         : 0.511   classifier: ranger   cv_folds: 5
+#>   ----------------------------------------
+#>   reproduction  : 98.9%
+```
+
+For a strict out-of-sample evaluation, hold the reference draws out
+*before* fitting:
+
+```r
+idx  <- sample.int(nrow(draws), 0.8 * nrow(draws))
+comp <- compress_posterior(draws[idx, ], method = "mclust")
+evaluate_compression(comp, reference_draws = draws[-idx, ])
+```
+
 ## Vignettes
 
 - `vignette("introduction", package = "poco")` - hello world,
@@ -129,6 +173,9 @@ new_draws <- sample_posterior(comp, n_draws = 4000)
   `sccomp` workflow.
 - `vignette("methods-benchmark", package = "poco")` - side-by-side
   comparison of the supported compression methods.
+- `vignette("evaluate-compression", package = "poco")` - measure
+  reproduction quality with `evaluate_compression()` (sanity checks
+  on lost correlation, mode collapse, and a sliding-difficulty sweep).
 
 ## License
 
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 5e5d27f..976305f 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -4,5 +4,6 @@ articles:
     contents:
       - introduction
       - methods-benchmark
+      - evaluate-compression
       - compress-funnel
       - compress-brms
diff --git a/man/evaluate_compression.Rd b/man/evaluate_compression.Rd
new file mode 100644
index 0000000..173db20
--- /dev/null
+++ b/man/evaluate_compression.Rd
@@ -0,0 +1,125 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/evaluate.R
+\name{evaluate_compression}
+\alias{evaluate_compression}
+\alias{print.compression_fidelity}
+\title{Evaluate the fidelity of a compressed posterior}
+\usage{
+evaluate_compression(
+  comp,
+  reference_draws,
+  metric = c("energy", "c2st"),
+  n_eval = NULL,
+  max_n = 2000L,
+  n_self_reps = 20L,
+  classifier = c("ranger", "knn"),
+  cv_folds = 5L,
+  seed = NULL,
+  verbose = FALSE
+)
+
+\method{print}{compression_fidelity}(x, ...)
+}
+\arguments{
+\item{comp}{A \code{posterior_compressed} object (or a path to an \code{.rds}
+file containing one).}
+
+\item{reference_draws}{A draws matrix, data.frame, or
+\verb{posterior::draws_*} object whose columns include all parameters
+in \code{comp}. Treated as the ground-truth distribution.}
+
+\item{metric}{Character vector of metrics to compute. One or more of
+\code{"energy"} and \code{"c2st"}. Default: both.}
+
+\item{n_eval}{Integer number of draws to regenerate from \code{comp}.
+Defaults to \code{min(max_n, nrow(reference_draws))}.}
+
+\item{max_n}{Cap on the number of points used in any pairwise
+distance / classifier computation. Both samples are subsampled to
+at most this many rows. Default \code{2000}.}
+
+\item{n_self_reps}{Number of self-baseline replicates for the energy
+metric. Default \code{20}.}
+
+\item{classifier}{Classifier for C2ST. \code{"ranger"} (default) uses a
+random forest (\code{ranger::ranger()}). \code{"knn"} uses a k-NN probability
+estimate instead (no random forest).}
+
+\item{cv_folds}{Cross-validation folds for C2ST. Default \code{5}.}
+
+\item{seed}{Optional integer seed for reproducibility.}
+
+\item{verbose}{Logical; print progress messages.}
+
+\item{x}{A \code{compression_fidelity} object.}
+
+\item{...}{Currently unused.}
+}
+\value{
+An S3 object of class \code{compression_fidelity} containing:
+\describe{
+\item{\code{reproduction_pct}}{Headline reproduction score in \verb{[0, 100]},
+averaged across requested metrics.}
+\item{\code{metrics}}{Named list with detailed per-metric results.}
+\item{\code{n_reference}, \code{n_eval}, \code{n_params}}{Sample sizes used.}
+}
+}
+\description{
+Compares draws regenerated from a compressed posterior against a
+reference set of draws (typically the original MCMC output) using
+two distribution-free, backend-independent two-sample diagnostics:
+}
+\details{
+\itemize{
+\item \strong{Energy distance} (Szekely & Rizzo, 2013): a
+scale-invariant distance between two samples that is zero
+iff their distributions match. The raw distance is anchored
+against a \emph{self-baseline} obtained by bootstrapping
+\code{reference_draws} to two independent samples matched in
+size to the reconstructed sample. The 90th percentile of
+this baseline distribution defines the Monte Carlo noise
+envelope; the reproduction score is
+\code{100 * min(1, noise_envelope / max(distance, noise_envelope))},
+so 100\\% means the reconstruction is within the bootstrap
+envelope and the score decays as 1/ratio for larger
+distances.
+\item \strong{Classifier two-sample test (C2ST)} (Lopez-Paz &
+Oquab, 2017): trains a classifier (random forest via
+\code{ranger::ranger()} under cross-validation
+to discriminate reference draws from reconstructed draws and
+reports the out-of-fold ROC AUC. AUC = 0.5 means the two
+samples are indistinguishable; the score is mapped to
+\verb{100 * (1 - 2 * |AUC - 0.5|)}.
+}
+
+Both metrics are computed only from samples and are independent of
+the compression backend (\code{mclust}, \code{mvdens_gmm}, \code{mvdens_kde}), so
+the score is not "circular": it does not depend on the same density
+model that produced the compression.
+
+For a strict out-of-sample evaluation, hold out a fraction of the
+draws \emph{before} fitting:
+
+\preformatted{
+  idx  <- sample.int(nrow(draws), size = 0.8 * nrow(draws))
+  comp <- compress_posterior(draws[idx, ], method = "mclust")
+  evaluate_compression(comp, reference_draws = draws[-idx, ])
+}
+}
+\examples{
+set.seed(1)
+draws <- matrix(rnorm(2000 * 3), ncol = 3,
+                dimnames = list(NULL, c("alpha", "beta", "sigma")))
+comp <- compress_posterior(draws, method = "mclust", n_components = 2)
+fidelity <- evaluate_compression(comp, reference_draws = draws, seed = 1L)
+fidelity
+
+}
+\references{
+Szekely, G. J. & Rizzo, M. L. (2013). Energy statistics: A class of
+statistics based on distances. \emph{Journal of Statistical Planning and
+Inference}, 143(8), 1249-1272.
+
+Lopez-Paz, D. & Oquab, M. (2017). Revisiting Classifier Two-Sample
+Tests. \emph{ICLR}.
+}
diff --git a/tests/testthat/test-evaluate.R b/tests/testthat/test-evaluate.R
new file mode 100644
index 0000000..5ea2a3c
--- /dev/null
+++ b/tests/testthat/test-evaluate.R
@@ -0,0 +1,175 @@
+test_that("evaluate_compression returns a compression_fidelity object", {
+  draws <- make_two_blob_draws()
+  comp <- compress_posterior(draws, method = "mclust", n_components = 2)
+
+  fid <- evaluate_compression(
+    comp,
+    reference_draws = draws,
+    seed            = 1L,
+    n_self_reps     = 5L,
+    max_n           = 600L
+  )
+
+  expect_s3_class(fid, "compression_fidelity")
+  expect_named(fid$metrics, c("energy", "c2st"))
+  expect_true(is.finite(fid$reproduction_pct))
+  expect_gte(fid$reproduction_pct, 0)
+  expect_lte(fid$reproduction_pct, 100)
+  expect_equal(fid$n_params, 2L)
+  expect_equal(fid$method, "mclust")
+})
+
+test_that("a faithful compression scores high on both metrics", {
+  draws <- make_two_blob_draws(n_per = 1500L)
+  comp <- compress_posterior(draws, method = "mclust", n_components = 2)
+
+  fid <- evaluate_compression(
+    comp,
+    reference_draws = draws,
+    seed            = 1L,
+    n_self_reps     = 5L,
+    max_n           = 600L
+  )
+
+  expect_gt(fid$metrics$energy$reproduction_pct, 80)
+  expect_gt(fid$metrics$c2st$reproduction_pct,   60)
+  expect_lt(abs(fid$metrics$c2st$auc - 0.5),     0.2)
+})
+
+test_that("a deliberately bad compression scores lower than a good one", {
+  draws <- make_two_blob_draws(n_per = 1500L)
+
+  good <- compress_posterior(draws, method = "mclust", n_components = 2)
+  bad  <- compress_posterior(draws, method = "mclust", n_components = 1L,
+                             model_name = "EII")
+
+  fid_good <- evaluate_compression(
+    good, reference_draws = draws,
+    metric = "energy", seed = 1L, n_self_reps = 5L, max_n = 600L
+  )
+  fid_bad  <- evaluate_compression(
+    bad,  reference_draws = draws,
+    metric = "energy", seed = 1L, n_self_reps = 5L, max_n = 600L
+  )
+
+  expect_gt(fid_good$reproduction_pct, fid_bad$reproduction_pct)
+})
+
+test_that("diagonal-only mclust scores lower than full covariance on a correlated posterior", {
+  set.seed(2)
+  Sigma <- matrix(c(1.0, 0.85, 0.85, 1.0), 2, 2)
+  draws <- mvtnorm::rmvnorm(2000, mean = c(0, 0), sigma = Sigma)
+  colnames(draws) <- c("alpha", "beta")
+
+  comp_diag <- compress_posterior(
+    draws, method = "mclust", n_components = 1L,
+    model_name = c("EII", "VII", "EEI", "VEI", "EVI", "VVI")
+  )
+  comp_full <- compress_posterior(
+    draws, method = "mclust", n_components = 1L,
+    model_name = c("EEE", "VVV")
+  )
+
+  fid_diag <- evaluate_compression(
+    comp_diag, draws,
+    seed = 1L, n_self_reps = 10L, max_n = 1000L
+  )
+  fid_full <- evaluate_compression(
+    comp_full, draws,
+    seed = 1L, n_self_reps = 10L, max_n = 1000L
+  )
+
+  expect_lt(fid_diag$reproduction_pct, fid_full$reproduction_pct - 20)
+  expect_lt(
+    fid_diag$metrics$energy$reproduction_pct,
+    fid_full$metrics$energy$reproduction_pct
+  )
+  expect_lt(
+    fid_diag$metrics$c2st$reproduction_pct,
+    fid_full$metrics$c2st$reproduction_pct
+  )
+
+  expect_gt(fid_full$reproduction_pct, 90)
+  expect_lt(fid_diag$reproduction_pct, 70)
+})
+
+test_that("evaluate_compression accepts a single metric", {
+  draws <- make_two_blob_draws()
+  comp <- compress_posterior(draws, method = "mclust", n_components = 2)
+
+  fid_e <- evaluate_compression(
+    comp, draws, metric = "energy",
+    seed = 1L, n_self_reps = 5L, max_n = 400L
+  )
+  expect_named(fid_e$metrics, "energy")
+
+  fid_c <- evaluate_compression(
+    comp, draws, metric = "c2st",
+    seed = 1L, max_n = 400L
+  )
+  expect_named(fid_c$metrics, "c2st")
+})
+
+test_that("evaluate_compression errors on mismatched parameters", {
+  draws <- make_two_blob_draws()
+  comp <- compress_posterior(draws, method = "mclust", n_components = 2)
+
+  bad_ref <- matrix(rnorm(20), ncol = 2,
+                    dimnames = list(NULL, c("foo", "bar")))
+  expect_error(evaluate_compression(comp, bad_ref), regexp = "matching|found")
+})
+
+test_that("evaluate_compression supports a held-out reference set", {
+  draws <- make_two_blob_draws(n_per = 1500L)
+  set.seed(42)
+  idx <- sample.int(nrow(draws), size = floor(0.8 * nrow(draws)))
+  comp <- compress_posterior(draws[idx, ], method = "mclust", n_components = 2)
+
+  fid <- evaluate_compression(
+    comp,
+    reference_draws = draws[-idx, ],
+    metric          = "energy",
+    seed            = 1L,
+    n_self_reps     = 5L,
+    max_n           = 400L
+  )
+  expect_s3_class(fid, "compression_fidelity")
+  expect_gt(fid$reproduction_pct, 50)
+})
+
+test_that("print.compression_fidelity prints a reproduction line", {
+  draws <- make_two_blob_draws()
+  comp <- compress_posterior(draws, method = "mclust", n_components = 2)
+
+  fid <- evaluate_compression(
+    comp, draws,
+    metric = "energy", seed = 1L, n_self_reps = 5L, max_n = 400L
+  )
+  expect_output(print(fid), "reproduction")
+})
+
+test_that("default C2ST classifier is ranger", {
+  draws <- make_two_blob_draws()
+  comp <- compress_posterior(draws, method = "mclust", n_components = 2)
+
+  fid <- evaluate_compression(
+    comp, draws,
+    metric = "c2st",
+    seed = 1L, max_n = 400L
+  )
+  expect_equal(fid$metrics$c2st$classifier, "ranger")
+  expect_true(is.finite(fid$metrics$c2st$auc))
+})
+
+test_that("classifier = knn uses k-NN path", {
+  draws <- make_two_blob_draws()
+  comp <- compress_posterior(draws, method = "mclust", n_components = 2)
+
+  fid <- evaluate_compression(
+    comp, draws,
+    metric = "c2st", classifier = "knn",
+    seed = 1L, max_n = 300L
+  )
+  expect_equal(fid$metrics$c2st$classifier, "knn")
+  expect_true(is.finite(fid$metrics$c2st$auc))
+})
diff --git a/vignettes/evaluate-compression.Rmd b/vignettes/evaluate-compression.Rmd
new file mode 100644
index 0000000..e6b2d63
--- /dev/null
+++ b/vignettes/evaluate-compression.Rmd
@@ -0,0 +1,435 @@
+---
+title: "Evaluating compression fidelity"
+output:
+  rmarkdown::html_vignette:
+    toc: true
+vignette: >
+  %\VignetteIndexEntry{Evaluating compression fidelity}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include = FALSE}
+knitr::opts_chunk$set(
+  collapse   = TRUE,
+  comment    = "#>",
+  fig.width  = 7,
+  fig.height = 4,
+  message    = FALSE,
+  warning    = FALSE
+)
+```
+
+Posterior compression is lossy by construction - just like JPEG for
+images. `evaluate_compression()` reports a **reproduction percentage**
+that quantifies how much information was lost, and it does so without
+being *circular*: the score does not reuse the same density model
+(`mclust`, `mvdens_gmm`, `mvdens_kde`) that performed the compression.
+
+This vignette walks through the metric on three sanity checks where the
+visual answer is obvious, then a sweep that pins down the metric's
+behaviour as we increase the difficulty of the problem.
+
+## Two non-circular yardsticks
+
+`evaluate_compression()` combines two distribution-free, sample-based,
+two-sample diagnostics:
+
+- **Energy distance** (Szekely & Rizzo, 2013): a multivariate distance
+  between two empirical samples that is zero iff the two distributions
+  match. It is anchored against a *bootstrap noise envelope* obtained
+  by resampling the reference draws, so a reproduction of 100% means
+  the regenerated sample is statistically indistinguishable from a
+  fresh draw of the reference.
+- **Classifier two-sample test (C2ST)** (Lopez-Paz & Oquab, 2017): a
+  random forest (`ranger::ranger()`, default) is trained under
+  cross-validation to discriminate reference from reconstructed draws.
+  Use `classifier = "knn"` for a k-NN classifier instead. AUC = 0.5 means
+  indistinguishable; the score is `100 * (1 - 2 * |AUC - 0.5|)`.
+
+Both metrics work purely from samples returned by `sample_posterior()`;
+they never inspect the GMM weights or KDE bandwidths.
+
+## Setup
+
+```{r libs}
+library(poco)
+library(ggplot2)
+library(patchwork)
+library(mvtnorm)
+
+set.seed(1)
+```
+
+```{r helpers, include = FALSE}
+plot_overlay <- function(reference, recon, title) {
+  ref_df   <- data.frame(reference, source = "reference")
+  recon_df <- data.frame(recon,     source = "reconstructed")
+  df <- rbind(ref_df, recon_df)
+  cols <- colnames(reference)
+  ggplot(df, aes(x = .data[[cols[1]]], y = .data[[cols[2]]],
+                 colour = source)) +
+    geom_point(alpha = 0.3, size = 0.6) +
+    stat_ellipse(
+      aes(group = source),
+      type = "norm",
+      level = 0.68,
+      linewidth = 0.7,
+      linetype = 2
+    ) +
+    scale_colour_manual(values = c(reference = "gray40",
+                                   reconstructed = "firebrick")) +
+    labs(title = title, x = cols[1], y = cols[2]) +
+    theme_bw() +
+    theme(legend.position = "bottom")
+}
+```
+
+# Sanity check 1: a faithful compression
+
+A bivariate Gaussian with strong correlation, compressed with a single
+**full-covariance** Gaussian (the textbook right answer).
+
+```{r demo1-data}
+Sigma <- matrix(c(1.0, 0.7, 0.7, 1.0), 2, 2)
+draws <- rmvnorm(2000, mean = c(2, -1), sigma = Sigma)
+colnames(draws) <- c("alpha", "beta")
+```
+
+```{r demo1-fit}
+comp_full <- compress_posterior(
+  draws, method = "mclust",
+  n_components = 1L, model_name = "VVV"
+)
+recon_full <- sample_posterior(comp_full, n_draws = 2000)
+
+fid_full <- evaluate_compression(comp_full, draws, seed = 1L)
+fid_full
+```
+
+```{r demo1-plot, fig.width = 6, fig.height = 5, echo = FALSE}
+plot_overlay(
+  draws, recon_full,
+  sprintf("Faithful (full covariance) - %.1f%% reproduction",
+          fid_full$reproduction_pct)
+)
+```
+
+The reconstructed cloud is visually indistinguishable from the
+reference. Both metrics agree it is at the noise floor.
+
+# Sanity check 2: the killer test - lost correlation
+
+Now compress the **same** correlated posterior with a *diagonal-only*
+mclust model. By construction, every component has zero off-diagonal
+covariance, so any joint correlation between `alpha` and `beta` simply
+cannot be represented.
+
+```{r demo2-fit}
+comp_diag <- compress_posterior(
+  draws, method = "mclust",
+  n_components = 1L,
+  model_name = c("EII", "VII", "EEI", "VEI", "EVI", "VVI")
+)
+recon_diag <- sample_posterior(comp_diag, n_draws = 2000)
+
+fid_diag <- evaluate_compression(comp_diag, draws, seed = 1L)
+fid_diag
+```
+
+```{r demo2-plot, fig.width = 11, fig.height = 5, echo = FALSE}
+(plot_overlay(draws, recon_full,
+              sprintf("Full covariance: %.1f%%",
+                      fid_full$reproduction_pct)) +
+ plot_overlay(draws, recon_diag,
+              sprintf("Diagonal only: %.1f%%",
+                      fid_diag$reproduction_pct)))
+```
+
+Visually, the orange cloud on the right is *circular* where the gray
+reference is *elongated*: the reconstruction has lost the correlation.
+The empirical Pearson correlations confirm what the eye sees:
+
+```{r demo2-corr}
+data.frame(
+  source           = c("reference", "diagonal recon", "full recon"),
+  cor_alpha_beta   = c(cor(draws)[1, 2],
+                       cor(recon_diag)[1, 2],
+                       cor(recon_full)[1, 2])
+) |>
+  transform(cor_alpha_beta = round(cor_alpha_beta, 3))
+```
+
+The reproduction score drops from `r round(fid_full$reproduction_pct)`%
+to `r round(fid_diag$reproduction_pct)`% - and crucially, **both**
+metrics independently catch it (energy
+`r round(fid_diag$metrics$energy$reproduction_pct)`%, C2ST
+`r round(fid_diag$metrics$c2st$reproduction_pct)`%):
+
+```{r demo2-bothmetrics, echo = FALSE}
+data.frame(
+  metric                 = c("energy", "c2st", "headline"),
+  full_covariance        = round(c(
+    fid_full$metrics$energy$reproduction_pct,
+    fid_full$metrics$c2st$reproduction_pct,
+    fid_full$reproduction_pct
+  ), 1),
+  diagonal_only          = round(c(
+    fid_diag$metrics$energy$reproduction_pct,
+    fid_diag$metrics$c2st$reproduction_pct,
+    fid_diag$reproduction_pct
+  ), 1)
+)
+```
+
+# Sanity check 3: a collapsed mode
+
+A bimodal posterior, compressed with two or one components.
+
+```{r demo3-data}
+m1 <- rmvnorm(1500, mean = c(-2, -1),
+              sigma = matrix(c(0.4, 0.0, 0.0, 0.3), 2, 2))
+m2 <- rmvnorm(1500, mean = c( 2,  1),
+              sigma = matrix(c(0.3, 0.0, 0.0, 0.5), 2, 2))
+bi  <- rbind(m1, m2)
+colnames(bi) <- c("alpha", "beta")
+```
+
+```{r demo3-fits}
+comp_two <- compress_posterior(bi, method = "mclust", n_components = 2L)
+comp_one <- compress_posterior(bi, method = "mclust", n_components = 1L)
+
+recon_two <- sample_posterior(comp_two, n_draws = 3000)
+recon_one <- sample_posterior(comp_one, n_draws = 3000)
+
+fid_two <- evaluate_compression(comp_two, bi, seed = 1L)
+fid_one <- evaluate_compression(comp_one, bi, seed = 1L)
+```
+
+```{r demo3-plot, fig.width = 11, fig.height = 5, echo = FALSE}
+(plot_overlay(bi, recon_two,
+              sprintf("2 components: %.1f%%",
+                      fid_two$reproduction_pct)) +
+ plot_overlay(bi, recon_one,
+              sprintf("1 component (mode collapse): %.1f%%",
+                      fid_one$reproduction_pct)))
+```
+
+The single-component reconstruction visibly *bridges* the two real
+modes - and the score drops from
+`r round(fid_two$reproduction_pct)`% to
+`r round(fid_one$reproduction_pct)`%. Once again, both metrics agree
+something is wrong.
+
+# Sanity check 4: Neal's funnel (hard geometry)
+
+Neal's funnel is a classic stress test: one dimension controls the scale
+of the others, creating a *narrow neck* and strong non-linear geometry.
+Even in two dimensions, a density approximation that looks reasonable in
+the bulk can struggle in the neck.
+
+We generate draws from the canonical funnel:
+
+- \(y \sim \mathcal{N}(0, 3)\)
+- \(x \mid y \sim \mathcal{N}(0, \exp(y/2))\)
+
+```{r funnel-data}
+set.seed(1)
+n_funnel <- 4000
+y <- rnorm(n_funnel, mean = 0, sd = 3)
+x <- rnorm(n_funnel, mean = 0, sd = exp(y / 2))
+funnel <- cbind(y = y, x = x)
+```
+
+```{r funnel-fit}
+# Higher-capacity model (near-optimal for this example)
+comp_funnel_rich <- compress_posterior(
+  funnel, method = "mclust", n_components = 8L
+)
+recon_funnel_rich <- sample_posterior(comp_funnel_rich, n_draws = 4000)
+fid_funnel_rich <- evaluate_compression(
+  comp_funnel_rich,
+  reference_draws = funnel,
+  seed        = 1L,
+  max_n       = 1200L,
+  n_self_reps = 10L
+)
+
+# Oversimplified model 1: single full-covariance Gaussian
+comp_funnel_one <- compress_posterior(
+  funnel, method = "mclust", n_components = 1L, model_name = "VVV"
+)
+recon_funnel_one <- sample_posterior(comp_funnel_one, n_draws = 4000)
+fid_funnel_one <- evaluate_compression(
+  comp_funnel_one,
+  reference_draws = funnel,
+  seed        = 1L,
+  max_n       = 1200L,
+  n_self_reps = 10L
+)
+
+# Oversimplified model 2: single diagonal Gaussian
+comp_funnel_diag <- compress_posterior(
+  funnel, method = "mclust", n_components = 1L,
+  model_name = c("EII", "VII", "EEI", "VEI", "EVI", "VVI")
+)
+recon_funnel_diag <- sample_posterior(comp_funnel_diag, n_draws = 4000)
+fid_funnel_diag <- evaluate_compression(
+  comp_funnel_diag,
+  reference_draws = funnel,
+  seed        = 1L,
+  max_n       = 1200L,
+  n_self_reps = 10L
+)
+
+data.frame(
+  model = c("8 components (richer)", "1 component full-cov", "1 component diagonal"),
+  reproduction = round(c(
+    fid_funnel_rich$reproduction_pct,
+    fid_funnel_one$reproduction_pct,
+    fid_funnel_diag$reproduction_pct
+  ), 1),
+  energy = round(c(
+    fid_funnel_rich$metrics$energy$reproduction_pct,
+    fid_funnel_one$metrics$energy$reproduction_pct,
+    fid_funnel_diag$metrics$energy$reproduction_pct
+  ), 1),
+  c2st = round(c(
+    fid_funnel_rich$metrics$c2st$reproduction_pct,
+    fid_funnel_one$metrics$c2st$reproduction_pct,
+    fid_funnel_diag$metrics$c2st$reproduction_pct
+  ), 1)
+)
+```
+
+```{r funnel-plot, fig.width = 15, fig.height = 5, echo = FALSE}
+plot_overlay(
+  funnel, recon_funnel_rich,
+  sprintf("Funnel: 8 components - %.1f%%", fid_funnel_rich$reproduction_pct)
+) +
+plot_overlay(
+  funnel, recon_funnel_one,
+  sprintf("Funnel: 1 full-cov - %.1f%%", fid_funnel_one$reproduction_pct)
+) +
+plot_overlay(
+  funnel, recon_funnel_diag,
+  sprintf("Funnel: 1 diagonal - %.1f%%", fid_funnel_diag$reproduction_pct)
+)
+```
+
+As expected, the richer model tracks the funnel geometry better, while
+oversimplified one-component fits smooth out the neck and lose shape.
+
+# Sanity sweep: reproduction vs the difficulty of the problem
+
+If the metric is doing its job, increasing the *true* correlation
+should make a diagonal-only compression progressively worse. We sweep
+the correlation from 0 to 0.95 and watch the reproduction score.
+
+```{r sweep, fig.width = 8, fig.height = 4}
+sweep_one_rho <- function(rho, n = 1500, seed = 7L) {
+  set.seed(seed)
+  S <- matrix(c(1, rho, rho, 1), 2, 2)
+  d <- rmvnorm(n, sigma = S)
+  colnames(d) <- c("alpha", "beta")
+  comp <- compress_posterior(
+    d, method = "mclust", n_components = 1L,
+    model_name = c("EII", "VII", "EEI", "VEI", "EVI", "VVI")
+  )
+  fid <- evaluate_compression(
+    comp, d,
+    seed        = 1L,
+    n_self_reps = 10L,
+    max_n       = 800L
+  )
+  data.frame(
+    rho      = rho,
+    energy   = fid$metrics$energy$reproduction_pct,
+    c2st     = fid$metrics$c2st$reproduction_pct,
+    headline = fid$reproduction_pct
+  )
+}
+
+rhos <- seq(0, 0.95, by = 0.05)
+sweep_df <- do.call(rbind, lapply(rhos, sweep_one_rho))
+
+sweep_long <- rbind(
+  data.frame(rho = sweep_df$rho, metric = "energy",
+             pct = sweep_df$energy),
+  data.frame(rho = sweep_df$rho, metric = "c2st",
+             pct = sweep_df$c2st),
+  data.frame(rho = sweep_df$rho, metric = "headline (mean)",
+             pct = sweep_df$headline)
+)
+
+ggplot(sweep_long, aes(x = rho, y = pct, colour = metric)) +
+  geom_hline(yintercept = 100, linetype = 3, colour = "gray70") +
+  geom_line(linewidth = 0.9) +
+  geom_point(size = 1.5) +
+  scale_colour_manual(values = c(
+    "energy"           = "steelblue",
+    "c2st"             = "firebrick",
+    "headline (mean)"  = "gray30"
+  )) +
+  scale_y_continuous(limits = c(0, 105),
+                     breaks = seq(0, 100, by = 20)) +
+  labs(
+    title    = "Diagonal-only compression: fidelity drops as correlation grows",
+    subtitle = "Both metrics independently track the lost joint structure",
+    x        = expression("True correlation " * rho),
+    y        = "reproduction (%)",
+    colour   = "metric"
+  ) +
+  theme_bw()
+```
+
+Both energy and C2ST decrease monotonically with the strength of the
+joint structure that the diagonal model cannot capture. At `rho = 0`
+the diagonal model *is* the truth and reproduction sits at 100%; by
+`rho = 0.95` the metric correctly reports a substantial loss.
+
+# Strict out-of-sample evaluation
+
+If you want to rule out any residual concern that the GMM has seen the
+reference draws (it has, because we fit on the same `draws`), split
+*before* fitting:
+
+```{r holdout}
+n   <- nrow(bi)
+idx <- sample.int(n, size = floor(0.8 * n))
+comp_train <- compress_posterior(
+  bi[idx, ], method = "mclust", n_components = 2L
+)
+
+fid_holdout <- evaluate_compression(
+  comp_train,
+  reference_draws = bi[-idx, ],
+  seed            = 1L
+)
+fid_holdout
+```
+
+The interpretation is identical, but the reference is data the GMM
+never saw. Useful as a final sign-off that the score is not
+optimistically biased.
+
+# Wrap-up
+
+`evaluate_compression()` gives you:
+
+- a **single, interpretable number** (`reproduction_pct`) suitable for
+  reporting alongside compressed posteriors;
+- two **independent** corroborating metrics (energy distance and
+  C2ST), so a single high score is unlikely to be a fluke;
+- a workflow that is **backend-agnostic**: the same function works
+  identically for `mclust`, `mvdens_gmm`, and `mvdens_kde` outputs.
+
+The sanity checks above are also a useful template for authors of new
+compression backends to verify that a method has not silently lost
+correlation, multimodality, or any other joint structure.
+
+## Session info
+
+```{r si}
+sessionInfo()
+```