diff --git a/R/centrality.R b/R/centrality.R
index deda22e2..2b47acb6 100644
--- a/R/centrality.R
+++ b/R/centrality.R
@@ -82,18 +82,22 @@ extract_centrality = function(bgms_object, measure = "strength", group = 1, ...)
 # Strength centrality of every node, on each row of a draws x pairs matrix of
 # edge weights.
 #
-# The incidence comes from the edge order rather than from the column names: a
-# variable name containing a hyphen would defeat splitting "A-B" back into its
-# nodes.
+# The incidence comes from the caller's edge index rather than from the column
+# names: a variable name containing a hyphen would defeat splitting "A-B" back
+# into its nodes. It cannot be assumed either -- a mixed fit lays its pairwise
+# draws out by block (discrete-discrete, continuous-continuous, cross), which is
+# not the row-major upper triangle whenever the discrete and continuous columns
+# interleave, and summing the wrong columns gives each node another node's
+# edges. indicator_pair_index() (R/verdicts.R) reads the layout off the fit.
 #
 # @param samples        draws x pairs matrix of edge weights.
+# @param pairs          pairs x 2 matrix of (row, column) node positions, in
+#   the column order of `samples`.
 # @param num_variables  Number of nodes.
 #
 # Returns: a draws x nodes matrix.
 # ------------------------------------------------------------------
-strength_from_pairwise = function(samples, num_variables) {
-  pairs = which(upper.tri(matrix(0, num_variables, num_variables)), arr.ind = TRUE)
-  pairs = pairs[order(pairs[, "row"], pairs[, "col"]), , drop = FALSE]
+strength_from_pairwise = function(samples, pairs, num_variables) {
   stopifnot(nrow(pairs) == ncol(samples))
 
   absolute = abs(samples)
@@ -114,9 +118,10 @@ extract_centrality.bgms = function(bgms_object, measure = "strength", group = 1,
 
   samples = extract_pairwise_interactions(bgms_object)
   nodes = extract_arguments(bgms_object)$data_columnnames
+  pairs = indicator_pair_index(bgms_object, length(nodes))
 
   structure(
-    strength_from_pairwise(samples, length(nodes)),
+    strength_from_pairwise(samples, pairs, length(nodes)),
     dimnames = list(NULL, nodes),
     class = c("bgms_centrality", "matrix", "array"),
     measure = measure
@@ -168,13 +173,15 @@ extract_centrality.bgmCompare = function(bgms_object, measure = "strength",
     out
   }
 
+  pairs = indicator_pair_index(bgms_object, num_variables)
   labels = compare_group_labels(arguments, num_groups)
   if(length(group) == 1L) {
-    centrality = strength_from_pairwise(weights_of(group), num_variables)
+    centrality = strength_from_pairwise(weights_of(group), pairs, num_variables)
     label = sprintf("%s %s centrality", group_tag(labels, group), measure)
   } else {
-    centrality = strength_from_pairwise(weights_of(group[1]), num_variables) -
-      strength_from_pairwise(weights_of(group[2]), num_variables)
+    centrality =
+      strength_from_pairwise(weights_of(group[1]), pairs, num_variables) -
+      strength_from_pairwise(weights_of(group[2]), pairs, num_variables)
     label = sprintf(
       "difference in %s centrality (%s - %s)",
       measure, group_tag(labels, group[1]), group_tag(labels, group[2])
diff --git a/R/predict_simulate_ggm.R b/R/predict_simulate_ggm.R
index 544e0a94..c00fd793 100644
--- a/R/predict_simulate_ggm.R
+++ b/R/predict_simulate_ggm.R
@@ -7,23 +7,41 @@
 #   GGM Prediction Helpers
 # ==============================================================================
 
-# Reconstruct the full precision matrix from association-scale values and
-# residual variances.
+# Reconstruct the full precision matrix from association-scale values and the
+# precision diagonal.
 #
 # @param associations p x p symmetric matrix of pairwise
 #   associations (zero diagonal).
-# @param residual_variance Named numeric vector of residual variances.
+# @param precision_diagonal Numeric vector of diagonal precision elements.
 #
 # @return p x p precision matrix.
-reconstruct_precision = function(associations, residual_variance) {
+reconstruct_precision = function(associations, precision_diagonal) {
   omega = -2 * associations
   # Excluded edges (NA) have zero precision
   omega[is.na(omega)] = 0
-  diag(omega) = 1 / residual_variance
+  diag(omega) = precision_diagonal
   return(omega)
 }
 
 
+# The posterior mean precision diagonal, E[K_jj], read off the raw diagonal
+# draws the sampler stores (the same values build_precision_from_draw() puts on
+# the diagonal per draw).
+#
+# posterior_mean_residual_variance is E[1 / K_jj] -- the right summary of a
+# residual variance, and deliberately not 1 / E[K_jj] -- but its reciprocal is
+# a harmonic mean of the diagonal, and pairing that with off-diagonals that are
+# arithmetic means E[K_ij] gives a matrix that is no posterior mean of anything
+# and need not stay positive definite. The mean of the draws is.
+#
+# @param object Fitted bgms object (GGM).
+#
+# @return Numeric vector of length p.
+posterior_mean_precision_diagonal = function(object) {
+  colMeans(do.call(rbind, get_raw_samples(object)$main))
+}
+
+
 # Reconstruct precision matrix from a single posterior draw.
 # GGM raw samples are already on precision scale.
 #
@@ -74,7 +92,7 @@ simulate_bgms_ggm = function(object, nsim, seed, method, ndraws,
     # Reconstruct precision matrix from off-diagonal + separate diagonal
     precision = reconstruct_precision(
       get_posterior_mean(object, "pairwise"),
-      get_posterior_mean(object, "residual_variance")
+      posterior_mean_precision_diagonal(object)
     )
 
     # Call simulate_mrf with variable_type = "continuous"
@@ -160,7 +178,7 @@ predict_bgms_ggm = function(object, newdata, predict_vars, data_columnnames,
     # Reconstruct precision matrix from posterior means
     omega = reconstruct_precision(
       get_posterior_mean(object, "pairwise"),
-      get_posterior_mean(object, "residual_variance")
+      posterior_mean_precision_diagonal(object)
     )
 
     result = compute_conditional_ggm(
@@ -246,12 +264,17 @@ predict_bgms_ggm = function(object, newdata, predict_vars, data_columnnames,
   }
 
   if(type == "response") {
-    # Return conditional means
-    pred_matrix = sapply(result, function(m) m[, "mean"])
-    if(is.vector(pred_matrix)) {
-      pred_matrix = matrix(pred_matrix, ncol = 1)
-    }
+    # Return conditional means. Preallocated rather than sapply()'d: with one
+    # row of newdata sapply() collapses to a vector, and the n x 1 matrix it
+    # was then reshaped into no longer takes one name per predicted variable.
+    pred_matrix = matrix(
+      NA_real_,
+      nrow = nrow(newdata_centered), ncol = length(predict_vars)
+    )
     colnames(pred_matrix) = data_columnnames[predict_vars]
+    for(v in seq_along(predict_vars)) {
+      pred_matrix[, v] = result[[v]][, "mean"]
+    }
     return(pred_matrix)
   }
 
diff --git a/R/predict_simulate_mixed.R b/R/predict_simulate_mixed.R
index 97c3aba3..6948c411 100644
--- a/R/predict_simulate_mixed.R
+++ b/R/predict_simulate_mixed.R
@@ -194,6 +194,16 @@ predict_bgms_mixed = function(object, newdata, predict_vars, arguments,
     }
   }
 
+  # Rows the kernel cannot condition on. A missing discrete value reaches
+  # compute_conditional_mixed() as the NA_integer_ sentinel and corrupts the
+  # rest score of every other variable in that row, discrete and continuous
+  # alike; the continuous block needs no such entry, since a missing continuous
+  # value propagates as NaN on its own.
+  na_mask = na_conditioning_mask(
+    cbind(is.na(x_data), matrix(FALSE, nrow = nrow(x_data), ncol = q)),
+    internal_predict_vars + 1L
+  )
+
   compute_one_draw = function(pairwise_disc, pairwise_cross, pairwise_cont, mux, muy) {
     compute_conditional_mixed(
       x_observations = x_data,
@@ -218,7 +228,8 @@ predict_bgms_mixed = function(object, newdata, predict_vars, arguments,
 
     probs = format_mixed_predictions(
       raw_result, predict_vars, internal_predict_vars,
-      p, num_categories, data_columnnames
+      p, num_categories, data_columnnames,
+      category_levels, blume_capel_shift
     )
   } else {
     sample_info = split_mixed_raw_samples(object, arguments)
@@ -239,11 +250,15 @@ predict_bgms_mixed = function(object, newdata, predict_vars, arguments,
     }
 
     if(isTRUE(return_draws)) {
-      return(lapply(all_results, format_mixed_predictions,
+      per_draw = lapply(all_results, format_mixed_predictions,
         predict_vars = predict_vars,
         internal_predict_vars = internal_predict_vars, p = p,
-        num_categories = num_categories, data_columnnames = data_columnnames
-      ))
+        num_categories = num_categories, data_columnnames = data_columnnames,
+        category_levels = category_levels,
+        blume_capel_shift = blume_capel_shift
+      )
+      warn_na_conditioning(na_mask)
+      return(lapply(per_draw, apply_na_conditioning_mask, mask = na_mask))
     }
 
     # Average predictions across draws
@@ -261,16 +276,25 @@ predict_bgms_mixed = function(object, newdata, predict_vars, arguments,
 
     probs = format_mixed_predictions(
       probs, predict_vars, internal_predict_vars,
-      p, num_categories, data_columnnames
+      p, num_categories, data_columnnames,
+      category_levels, blume_capel_shift
     )
     names(probs_sd) = names(probs)
     attr(probs, "sd") = probs_sd
   }
 
+  probs_sd = attr(probs, "sd")
+  probs = apply_na_conditioning_mask(probs, na_mask)
+  if(!is.null(probs_sd)) {
+    attr(probs, "sd") = apply_na_conditioning_mask(probs_sd, na_mask)
+  }
+  warn_na_conditioning(na_mask)
+
   if(type == "response") {
     return(format_mixed_response(
       probs, predict_vars, internal_predict_vars,
-      p, data_columnnames
+      p, data_columnnames, num_categories,
+      category_levels, blume_capel_shift
     ))
   }
 
@@ -321,10 +345,13 @@ build_mixed_params_mean = function(object, arguments) {
       pairwise_cont[i, j] = pmat[cont_idx[i], cont_idx[j]]
     }
   }
-  # Convert residual variance back to association-scale diagonal
-  rv = get_posterior_mean(object, "residual_variance")
+  # Association-scale diagonal, from the draws that carry it. Going through
+  # posterior_mean_residual_variance instead would put a harmonic mean of the
+  # diagonal next to arithmetic-mean off-diagonals; see
+  # posterior_mean_precision_diagonal() in predict_simulate_ggm.R.
+  diag_draws = mixed_cont_diagonal_draws(object, arguments)
   for(j in seq_len(q)) {
-    pairwise_cont[j, j] = -1 / (2 * rv[j])
+    pairwise_cont[j, j] = mean(diag_draws[, j])
   }
 
   pm_main = get_posterior_mean(object, "main")
@@ -337,6 +364,29 @@ build_mixed_params_mean = function(object, arguments) {
 }
 
 
+# ------------------------------------------------------------------
+# mixed_cont_diagonal_draws
+# ------------------------------------------------------------------
+# The raw continuous-diagonal draws of a mixed fit, on the association scale
+# the sampler stores them on (entry -K_jj/2). They sit at the tail of the main
+# draws, after the discrete thresholds and the continuous means -- the same
+# layout split_mixed_raw_samples() reads.
+#
+# @param object     Fitted bgms object (mixed MRF).
+# @param arguments  Output of extract_arguments().
+#
+# Returns: draws x q matrix.
+# ------------------------------------------------------------------
+mixed_cont_diagonal_draws = function(object, arguments) {
+  q = arguments$num_continuous
+  num_mux = sum(ifelse(
+    arguments$is_ordinal, arguments$num_categories, 2L
+  ))
+  main_all = do.call(rbind, get_raw_samples(object)$main)
+  main_all[, num_mux + q + seq_len(q), drop = FALSE]
+}
+
+
 # ------------------------------------------------------------------
 # split_mixed_raw_samples
 # ------------------------------------------------------------------
@@ -539,7 +589,9 @@ combine_mixed_result = function(result, disc_idx, cont_idx, colnames) {
 # ------------------------------------------------------------------
 format_mixed_predictions = function(raw_result, predict_vars,
                                     internal_predict_vars, p,
-                                    num_categories, data_columnnames) {
+                                    num_categories, data_columnnames,
+                                    category_levels = NULL,
+                                    blume_capel_shift = NULL) {
   probs = raw_result
   names(probs) = data_columnnames[predict_vars]
 
@@ -547,8 +599,9 @@ format_mixed_predictions = function(raw_result, predict_vars,
     int_idx = internal_predict_vars[k]
     if(int_idx < p) {
       s = int_idx + 1L
-      n_cats = num_categories[s] + 1
-      colnames(probs[[k]]) = paste0("cat_", 0:(n_cats - 1))
+      colnames(probs[[k]]) = probability_column_labels(original_category_values(
+        s, num_categories[s], category_levels, blume_capel_shift
+      ))
     } else {
       colnames(probs[[k]]) = c("mean", "sd")
     }
@@ -573,7 +626,9 @@ format_mixed_predictions = function(raw_result, predict_vars,
 # ------------------------------------------------------------------
 format_mixed_response = function(probs, predict_vars,
                                  internal_predict_vars, p,
-                                 data_columnnames) {
+                                 data_columnnames, num_categories,
+                                 category_levels = NULL,
+                                 blume_capel_shift = NULL) {
   n = nrow(probs[[1]])
   out = matrix(NA_real_, nrow = n, ncol = length(predict_vars))
   colnames(out) = data_columnnames[predict_vars]
@@ -581,7 +636,11 @@ format_mixed_response = function(probs, predict_vars,
   for(k in seq_along(predict_vars)) {
     int_idx = internal_predict_vars[k]
     if(int_idx < p) {
-      out[, k] = apply(probs[[k]], 1, which.max) - 1L
+      s = int_idx + 1L
+      values = original_category_values(
+        s, num_categories[s], category_levels, blume_capel_shift
+      )
+      out[, k] = values[row_modes(probs[[k]])]
     } else {
       out[, k] = probs[[k]][, 1]
     }
diff --git a/R/simulate_mrf.R b/R/simulate_mrf.R
index 110aa8c8..f846c254 100644
--- a/R/simulate_mrf.R
+++ b/R/simulate_mrf.R
@@ -481,7 +481,6 @@ simulate_mrf = function(num_states,
 }
 
 
-
 # ==============================================================================
 #   mrfSampler() - Deprecated Wrapper for simulate_mrf()
 # ==============================================================================
diff --git a/R/simulate_predict.R b/R/simulate_predict.R
index af07ba3f..b04a24b4 100644
--- a/R/simulate_predict.R
+++ b/R/simulate_predict.R
@@ -130,10 +130,22 @@ simulate.bgms = function(object,
   method = match.arg(method)
   progress_type = progress_type_from_display_progress(display_progress)
 
-  # Validate cores
-  check_positive_integer(cores, "cores")
+  # Validate cores. parallel::detectCores() -- the default -- is documented to
+  # return NA when it cannot tell, and check_positive_integer() turns that NA
+  # into "missing value where TRUE/FALSE needed" before
+  # normalize_parallel_cores() gets the chance to fall back to 1. The check runs
+  # on everything the user can actually supply; only that one NA is left to the
+  # normalizer.
+  if(!(length(cores) == 1L && is.na(cores))) {
+    check_positive_integer(cores, "cores")
+  }
   cores = normalize_parallel_cores(cores)
 
+  # nsim and iter reach four different simulation paths (OMRF mean, OMRF
+  # sample, GGM, mixed); validating them here covers all of them at once.
+  check_positive_integer(nsim, "nsim")
+  check_positive_integer(iter, "iter")
+
   # Setting the seed
   seed = check_seed(seed)
 
@@ -486,10 +498,12 @@ average_draws = function(per_draw_list, v) {
 #' predicted variable. Each element is a matrix with \code{n} rows and
 #' \code{num_categories + 1} columns containing
 #' \eqn{P(X_j = c | X_{-j})}{P(X_j = c | X_-j)}
-#' for each observation and category.
+#' for each observation and category. Columns are labelled with the original
+#' category values of the training data.
 #'
 #' For \code{type = "response"}: A matrix with \code{n} rows and
-#' \code{length(variables)} columns containing predicted categories.
+#' \code{length(variables)} columns containing predicted categories, on the
+#' original category scale that \code{\link{simulate.bgms}} returns.
 #'
 #' When \code{method = "posterior-sample"}, probabilities are averaged over
 #' posterior draws, and an attribute \code{"sd"} is included containing the
@@ -529,6 +543,18 @@ average_draws = function(per_draw_list, v) {
 #' and variance \eqn{\omega_{jj}^{-1}}{omega_jj^{-1}}, where \eqn{\Omega}{Omega}
 #' is the precision matrix.
 #'
+#' \code{newdata} is matched to the fitted model by position. When it carries
+#' column names they must be the model's variables in the model's order, or
+#' \code{predict()} stops; when it carries none it is read positionally, with a
+#' warning saying so.
+#'
+#' A discrete cell of \code{newdata} that is \code{NA}, or that holds a category
+#' value never observed in the training data, leaves the conditional
+#' distribution of every other variable in that row undefined. Those
+#' predictions are returned as \code{NA}, with one warning giving the number of
+#' rows affected; the variable whose own value is missing is unaffected, since
+#' its conditional distribution does not use it.
+#'
 #' @seealso \code{\link{simulate.bgms}} for generating new data from the model.
 #' @family prediction
 #'
@@ -600,6 +626,8 @@ predict.bgms = function(object,
     ))
   }
 
+  check_newdata_columns(newdata, data_columnnames)
+
   # Handle variable_type
   variable_type = expand_variable_type(variable_type, num_variables)
 
@@ -694,9 +722,10 @@ predict.bgms = function(object,
     # Add names
     names(probs) = data_columnnames[predict_vars]
     for(v in seq_along(probs)) {
-      var_idx = predict_vars[v]
-      n_cats = num_categories[var_idx] + 1
-      colnames(probs[[v]]) = paste0("cat_", 0:(n_cats - 1))
+      colnames(probs[[v]]) = probability_column_labels(original_category_values(
+        predict_vars[v], num_categories[predict_vars[v]],
+        arguments$category_levels, arguments$blume_capel_shift
+      ))
     }
   } else {
     # Use posterior samples
@@ -753,25 +782,32 @@ predict.bgms = function(object,
       probs[[v]] = avg$mean
       probs_sd[[v]] = avg$sd
 
-      var_idx = predict_vars[v]
-      n_cats = num_categories[var_idx] + 1
-      colnames(probs[[v]]) = paste0("cat_", 0:(n_cats - 1))
-      colnames(probs_sd[[v]]) = paste0("cat_", 0:(n_cats - 1))
+      labels = probability_column_labels(original_category_values(
+        predict_vars[v], num_categories[predict_vars[v]],
+        arguments$category_levels, arguments$blume_capel_shift
+      ))
+      colnames(probs[[v]]) = labels
+      colnames(probs_sd[[v]]) = labels
     }
 
     attr(probs, "sd") = probs_sd
   }
 
+  # Blank the rows the kernel could not condition on, before anything is read
+  # off them.
+  mask = na_conditioning_mask(is.na(newdata_recoded), predict_vars)
+  probs_sd = attr(probs, "sd")
+  probs = apply_na_conditioning_mask(probs, mask)
+  if(!is.null(probs_sd)) {
+    attr(probs, "sd") = apply_na_conditioning_mask(probs_sd, mask)
+  }
+  warn_na_conditioning(mask)
+
   if(type == "response") {
-    # Return predicted categories (mode)
-    pred_matrix = sapply(probs, function(p) {
-      apply(p, 1, which.max) - 1L # Convert to 0-based category
-    })
-    if(is.vector(pred_matrix)) {
-      pred_matrix = matrix(pred_matrix, ncol = 1)
-    }
-    colnames(pred_matrix) = data_columnnames[predict_vars]
-    return(pred_matrix)
+    return(format_discrete_response(
+      probs, predict_vars, data_columnnames, num_categories,
+      arguments$category_levels, arguments$blume_capel_shift
+    ))
   }
 
   return(probs)
@@ -819,10 +855,12 @@ predict.bgms = function(object,
 #' element per predicted variable. Each element is a matrix with
 #' \code{n} rows and \code{num_categories + 1} columns containing
 #' \eqn{P(X_j = c | X_{-j})}{P(X_j = c | X_-j)}
-#' for each observation and category.
+#' for each observation and category. Columns are labelled with the original
+#' category values of the training data.
 #'
 #' For \code{type = "response"}: A matrix with \code{n} rows and
-#' \code{length(variables)} columns containing predicted categories.
+#' \code{length(variables)} columns containing predicted categories, on the
+#' original category scale that \code{\link{simulate.bgmCompare}} returns.
 #'
 #' @details
 #' Group-specific parameters are obtained by applying the projection matrix
@@ -830,6 +868,9 @@ predict.bgms = function(object,
 #' The function then computes the conditional distribution of target variables
 #' given the observed values of all other variables.
 #'
+#' The \code{newdata} column-name and missing-value rules of
+#' \code{\link{predict.bgms}} apply here unchanged.
+#'
 #' @seealso \code{\link{predict.bgms}} for predicting
 #'   from single-group models,
 #'   \code{\link{simulate.bgmCompare}} for simulating
@@ -915,6 +956,8 @@ predict.bgmCompare = function(object,
     ))
   }
 
+  check_newdata_columns(newdata, data_columnnames)
+
   # Determine variable_type from is_ordinal
   variable_type = ifelse(is_ordinal, "ordinal", "blume-capel")
 
@@ -984,22 +1027,24 @@ predict.bgmCompare = function(object,
     # Add names
     names(probs) = data_columnnames[predict_vars]
     for(v in seq_along(probs)) {
-      var_idx = predict_vars[v]
-      n_cats = num_categories[var_idx] + 1
-      colnames(probs[[v]]) = paste0("cat_", 0:(n_cats - 1))
+      colnames(probs[[v]]) = probability_column_labels(original_category_values(
+        predict_vars[v], num_categories[predict_vars[v]],
+        arguments$category_levels, arguments$blume_capel_shift
+      ))
     }
   }
 
+  # Blank the rows the kernel could not condition on, before anything is read
+  # off them.
+  mask = na_conditioning_mask(is.na(newdata_recoded), predict_vars)
+  probs = apply_na_conditioning_mask(probs, mask)
+  warn_na_conditioning(mask)
+
   if(type == "response") {
-    # Return predicted categories (mode)
-    pred_matrix = sapply(probs, function(p) {
-      apply(p, 1, which.max) - 1L # Convert to 0-based category
-    })
-    if(is.vector(pred_matrix)) {
-      pred_matrix = matrix(pred_matrix, ncol = 1)
-    }
-    colnames(pred_matrix) = data_columnnames[predict_vars]
-    return(pred_matrix)
+    return(format_discrete_response(
+      probs, predict_vars, data_columnnames, num_categories,
+      arguments$category_levels, arguments$blume_capel_shift
+    ))
   }
 
   return(probs)
@@ -1073,18 +1118,24 @@ recode_data_for_prediction = function(x, is_ordinal,
       if(any(observed & is.na(recoded))) {
         warning(
           "newdata for variable ", v, " contains category values not ",
-          "observed in the training data; predictions for those cells are NA.",
+          "observed in the training data; those cells are treated as missing.",
           call. = FALSE
         )
       }
       x[, v] = recoded
     } else {
-      # Legacy fallback (fit has no recode map): shift to 0-based by the
-      # per-column minimum.
-      x[, v] = as.integer(x[, v])
-      if(min(x[, v], na.rm = TRUE) > 0) {
-        x[, v] = x[, v] - min(x[, v], na.rm = TRUE)
-      }
+      # No recode map for an ordinal variable. Every fit bgm() and bgmCompare()
+      # produce carries one, so this is only reachable for an object built by a
+      # bgms old enough to predate it. The shift-by-the-column-minimum fallback
+      # that used to stand here read the offset off newdata rather than off the
+      # training data, which is the wrong number whenever newdata does not
+      # happen to span the training range, and silently so.
+      stop(
+        "The fitted object carries no category recode map for variable ", v,
+        ", so 'newdata' cannot be put on the scale the model was fitted on. ",
+        "It predates the recode map; refit with the current bgms.",
+        call. = FALSE
+      )
     }
   }
 
@@ -1092,6 +1143,228 @@ recode_data_for_prediction = function(x, is_ordinal,
 }
 
 
+# ------------------------------------------------------------------------------
+# check_newdata_columns()
+# ------------------------------------------------------------------------------
+# newdata is matched to the fit by position, and a column count is not enough to
+# establish that the match is the intended one: a data frame whose columns were
+# reordered, or one built from a different subset of the same width, passes the
+# count check and then predicts every variable from the wrong neighbours,
+# silently. Named newdata therefore has to agree with the fit exactly. Unnamed
+# newdata is still accepted -- it carries nothing to check -- but says so.
+#
+# @param newdata           The (already coerced) newdata matrix.
+# @param data_columnnames  The fit's variable names.
+# ------------------------------------------------------------------------------
+check_newdata_columns = function(newdata, data_columnnames) {
+  if(is.null(data_columnnames)) {
+    return(invisible(NULL))
+  }
+
+  observed = colnames(newdata)
+  if(is.null(observed)) {
+    warning(
+      "'newdata' has no column names, so its columns are matched to the ",
+      "fitted model by position: ",
+      paste(data_columnnames, collapse = ", "), ".",
+      call. = FALSE
+    )
+    return(invisible(NULL))
+  }
+
+  if(identical(observed, data_columnnames)) {
+    return(invisible(NULL))
+  }
+
+  if(setequal(observed, data_columnnames)) {
+    stop(
+      "'newdata' holds the fitted model's variables in a different order. ",
+      "Columns are matched by position, so reorder 'newdata' to: ",
+      paste(data_columnnames, collapse = ", "), ".",
+      call. = FALSE
+    )
+  }
+
+  absent = setdiff(data_columnnames, observed)
+  unexpected = setdiff(observed, data_columnnames)
+  stop(
+    "'newdata' column names do not match the fitted model.",
+    if(length(absent) > 0) {
+      paste0(" Missing: ", paste(absent, collapse = ", "), ".")
+    },
+    if(length(unexpected) > 0) {
+      paste0(" Not in the model: ", paste(unexpected, collapse = ", "), ".")
+    },
+    call. = FALSE
+  )
+}
+
+
+# ------------------------------------------------------------------------------
+# original_category_values()
+# ------------------------------------------------------------------------------
+# The original category values behind the internal codes 0..num_categories of
+# one variable, in code order: the per-variable inverse of
+# recode_data_for_prediction(). simulate() returns data on this scale
+# (recode_simulated_to_original()), so predict() has to report responses and
+# label probability columns on it too, or the round trip the two promise each
+# other does not close for 1-based, non-contiguous or Blume-Capel variables.
+#
+# @param v                  Variable index.
+# @param num_categories     Number of categories on top of the baseline, for v.
+# @param category_levels    The fit's recode map (list, one entry per variable).
+# @param blume_capel_shift  The fit's Blume-Capel shifts (NA elsewhere).
+#
+# Returns: numeric vector of length num_categories + 1.
+# ------------------------------------------------------------------------------
+original_category_values = function(v, num_categories,
+                                    category_levels = NULL,
+                                    blume_capel_shift = NULL) {
+  codes = seq.int(0L, num_categories)
+
+  if(!is.null(blume_capel_shift) && !is.na(blume_capel_shift[v])) {
+    return(codes + blume_capel_shift[v])
+  }
+
+  levels_v = if(!is.null(category_levels)) category_levels[[v]] else NULL
+  if(is.null(levels_v)) {
+    return(as.numeric(codes))
+  }
+
+  if(!is.null(names(levels_v))) {
+    # Named lookup (bgmCompare), possibly many-to-one: invert to the smallest
+    # original value carrying each code, as recode_simulated_to_original() does,
+    # so the value reported recodes back to the code it came from.
+    inverse = tapply(
+      as.numeric(names(levels_v)), as.integer(unname(levels_v)), min
+    )
+    return(unname(inverse[as.character(codes)]))
+  }
+
+  as.numeric(levels_v[codes + 1L])
+}
+
+
+# ------------------------------------------------------------------------------
+# probability_column_labels()
+# ------------------------------------------------------------------------------
+# Column labels for a discrete variable's probability matrix, naming the
+# ORIGINAL category values rather than the internal codes.
+# ------------------------------------------------------------------------------
+probability_column_labels = function(values) {
+  paste0("cat_", values)
+}
+
+
+# ------------------------------------------------------------------------------
+# na_conditioning_mask()
+# ------------------------------------------------------------------------------
+# Which rows leave a target variable's conditional distribution undefined,
+# because some OTHER variable it conditions on is missing.
+#
+# The C++ kernels take the recoded observation matrix at face value: an NA cell
+# reaches them as R's NA_integer_ sentinel and enters the rest score of every
+# other variable in that row as a huge negative number, which comes back as a
+# confident one-hot distribution rather than as missingness. Both origins of the
+# NA -- a plain NA in newdata and a category value never seen in training -- are
+# caught here, so the masking is done once in R after the kernel returns.
+#
+# @param is_missing    n x p logical matrix, TRUE where the conditioning value
+#   is missing (in the kernel's own column order).
+# @param target_cols   Column of `is_missing` each prediction targets.
+#
+# Returns: list of logical vectors, one per target, TRUE for rows to blank.
+# ------------------------------------------------------------------------------
+na_conditioning_mask = function(is_missing, target_cols) {
+  total = rowSums(is_missing)
+  lapply(target_cols, function(v) total - is_missing[, v] > 0)
+}
+
+
+# ------------------------------------------------------------------------------
+# apply_na_conditioning_mask()
+# ------------------------------------------------------------------------------
+# Blank the masked rows of each prediction matrix and say how many rows of
+# newdata lost a prediction. Warning is the caller's, once per call.
+# ------------------------------------------------------------------------------
+apply_na_conditioning_mask = function(predictions, mask) {
+  for(k in seq_along(predictions)) {
+    rows = mask[[k]]
+    if(any(rows)) {
+      predictions[[k]][rows, ] = NA_real_
+    }
+  }
+  predictions
+}
+
+
+warn_na_conditioning = function(mask) {
+  affected = sum(Reduce(`|`, mask))
+  if(affected > 0) {
+    warning(
+      "newdata has ", affected, " row(s) in which a conditioning variable is ",
+      "missing or carries a category value not observed in the training data. ",
+      "The conditional distribution of the other variables is undefined in ",
+      "those rows, so their predictions are NA.",
+      call. = FALSE
+    )
+  }
+  invisible(affected)
+}
+
+
+# ------------------------------------------------------------------------------
+# row_modes()
+# ------------------------------------------------------------------------------
+# Index of the largest entry in each row (ties to the first), NA for rows that
+# carry any NA -- where which.max() would return integer(0) and collapse the
+# result to a list.
+# ------------------------------------------------------------------------------
+row_modes = function(probabilities) {
+  out = rep(NA_integer_, nrow(probabilities))
+  usable = !rowSums(is.na(probabilities)) > 0
+  if(any(usable)) {
+    out[usable] = max.col(
+      probabilities[usable, , drop = FALSE],
+      ties.method = "first"
+    )
+  }
+  out
+}
+
+
+# ------------------------------------------------------------------------------
+# format_discrete_response()
+# ------------------------------------------------------------------------------
+# Point predictions for discrete variables: the mode of each conditional
+# distribution, reported on the ORIGINAL category scale.
+#
+# The n x L result is preallocated rather than built with sapply(), which
+# collapses to a vector when n is 1 and then mis-shapes into an n x 1 matrix
+# that the column names no longer fit.
+# ------------------------------------------------------------------------------
+format_discrete_response = function(probabilities, predict_vars,
+                                    data_columnnames, num_categories,
+                                    category_levels = NULL,
+                                    blume_capel_shift = NULL) {
+  out = matrix(
+    NA_real_,
+    nrow = nrow(probabilities[[1]]), ncol = length(predict_vars)
+  )
+  colnames(out) = data_columnnames[predict_vars]
+
+  for(v in seq_along(predict_vars)) {
+    var_idx = predict_vars[v]
+    values = original_category_values(
+      var_idx, num_categories[var_idx], category_levels, blume_capel_shift
+    )
+    out[, v] = values[row_modes(probabilities[[v]])]
+  }
+
+  out
+}
+
+
 # ------------------------------------------------------------------------------
 # recode_simulated_to_original()
 # ------------------------------------------------------------------------------
diff --git a/man/predict.bgmCompare.Rd b/man/predict.bgmCompare.Rd
index 409bdc54..fb7aec03 100644
--- a/man/predict.bgmCompare.Rd
+++ b/man/predict.bgmCompare.Rd
@@ -51,10 +51,12 @@ For \code{type = "probabilities"}: A named list with one
 element per predicted variable. Each element is a matrix with
 \code{n} rows and \code{num_categories + 1} columns containing
 \eqn{P(X_j = c | X_{-j})}{P(X_j = c | X_-j)}
-for each observation and category.
+for each observation and category. Columns are labelled with the original
+category values of the training data.
 
 For \code{type = "response"}: A matrix with \code{n} rows and
-\code{length(variables)} columns containing predicted categories.
+\code{length(variables)} columns containing predicted categories, on the
+original category scale that \code{\link{simulate.bgmCompare}} returns.
 }
 \description{
 Computes conditional probability distributions for one or more variables
@@ -66,6 +68,9 @@ Group-specific parameters are obtained by applying the projection matrix
 to convert baseline parameters and differences into group-level estimates.
 The function then computes the conditional distribution of target variables
 given the observed values of all other variables.
+
+The \code{newdata} column-name and missing-value rules of
+\code{\link{predict.bgms}} apply here unchanged.
 }
 \examples{
 \donttest{
diff --git a/man/predict.bgms.Rd b/man/predict.bgms.Rd
index ec03a98b..8f65e221 100644
--- a/man/predict.bgms.Rd
+++ b/man/predict.bgms.Rd
@@ -60,10 +60,12 @@ For \code{type = "probabilities"}: A named list with one element per
 predicted variable. Each element is a matrix with \code{n} rows and
 \code{num_categories + 1} columns containing
 \eqn{P(X_j = c | X_{-j})}{P(X_j = c | X_-j)}
-for each observation and category.
+for each observation and category. Columns are labelled with the original
+category values of the training data.
 
 For \code{type = "response"}: A matrix with \code{n} rows and
-\code{length(variables)} columns containing predicted categories.
+\code{length(variables)} columns containing predicted categories, on the
+original category scale that \code{\link{simulate.bgms}} returns.
 
 When \code{method = "posterior-sample"}, probabilities are averaged over
 posterior draws, and an attribute \code{"sd"} is included containing the
@@ -107,6 +109,18 @@ For GGM (continuous) models, the conditional distribution of
 \omega_{jk} x_k}{-omega_jj^{-1} sum_{k != j} omega_jk x_k}
 and variance \eqn{\omega_{jj}^{-1}}{omega_jj^{-1}}, where \eqn{\Omega}{Omega}
 is the precision matrix.
+
+\code{newdata} is matched to the fitted model by position. When it carries
+column names they must be the model's variables in the model's order, or
+\code{predict()} stops; when it carries none it is read positionally, with a
+warning saying so.
+
+A discrete cell of \code{newdata} that is \code{NA}, or that holds a category
+value never observed in the training data, leaves the conditional
+distribution of every other variable in that row undefined. Those
+predictions are returned as \code{NA}, with one warning giving the number of
+rows affected; the variable whose own value is missing is unaffected, since
+its conditional distribution does not use it.
 }
 \examples{
 \donttest{
diff --git a/src/mrf_simulation.cpp b/src/mrf_simulation.cpp
index 75cccc47..2defc98f 100644
--- a/src/mrf_simulation.cpp
+++ b/src/mrf_simulation.cpp
@@ -54,6 +54,11 @@ arma::imat simulate_mrf(
   arma::imat observations(num_states, num_variables);
   int max_num_categories = arma::max(num_categories);
   arma::vec probabilities(max_num_categories + 1);
+  // Category exponents of the full conditional, held so the largest can be
+  // subtracted before exponentiating (as compute_probs_ordinal() and
+  // compute_probs_blume_capel() already do on the prediction side). Without it
+  // a large rest score overflows exp() to Inf and the draw degenerates.
+  arma::vec exponents(max_num_categories + 1);
   double exponent = 0.0;
   double rest_score = 0.0;
   double cumsum = 0.0;
@@ -97,8 +102,12 @@ arma::imat simulate_mrf(
           rest_score += 2.0 * (obs - ref) * pairwise_safe(vertex, variable);
         }
 
+        // Both branches collect the category exponents first, subtract their
+        // maximum, and only then exponentiate: the cumulative weights are
+        // sampled from up to a common factor, so the shift leaves the
+        // distribution untouched while keeping every term in range.
+        double max_exponent = 0.0;
         if(variable_type[variable] == "blume-capel") {
-          cumsum = 0.0;
           int ref = baseline_category[variable];
           for(int category = 0; category <= num_categories[variable]; category++) {
             const int s = category - ref;
@@ -108,18 +117,28 @@ arma::imat simulate_mrf(
             exponent += main(variable, 1) * s * s;
             // Pairwise effects
             exponent += rest_score * s;
-            cumsum += MY_EXP(exponent);
+            exponents[category] = exponent;
+            if(category == 0 || exponent > max_exponent) max_exponent = exponent;
+          }
+          cumsum = 0.0;
+          for(int category = 0; category <= num_categories[variable]; category++) {
+            cumsum += MY_EXP(exponents[category] - max_exponent);
             probabilities[category] = cumsum;
           }
         } else {
-          // Ordinal: baseline category 0 has probability 1 (unnormalized)
-          cumsum = 1.0;
-          probabilities[0] = cumsum;
+          // Ordinal: baseline category 0 carries exponent 0
+          exponents[0] = 0.0;
+          max_exponent = 0.0;
           for(int category = 0; category < num_categories[variable]; category++) {
             exponent = main(variable, category);
             exponent += (category + 1) * rest_score;
-            cumsum += MY_EXP(exponent);
-            probabilities[category + 1] = cumsum;
+            exponents[category + 1] = exponent;
+            if(exponent > max_exponent) max_exponent = exponent;
+          }
+          cumsum = 0.0;
+          for(int category = 0; category <= num_categories[variable]; category++) {
+            cumsum += MY_EXP(exponents[category] - max_exponent);
+            probabilities[category] = cumsum;
           }
         }
 
@@ -507,6 +526,14 @@ Rcpp::List run_simulation_parallel(
   }
   pm.finish();
 
+  // Every worker returns as soon as the interrupt flag is set, so the draws it
+  // had not reached hold default-constructed matrices. Those are not a shorter
+  // simulation, they are empty ones, so the interrupt is raised rather than
+  // handed back as data.
+  if (pm.shouldExit()) {
+    Rcpp::stop("Simulation interrupted by user.");
+  }
+
   // Convert results to R list
   Rcpp::List output(ndraws);
   for (int i = 0; i < ndraws; i++) {
@@ -664,6 +691,12 @@ Rcpp::List run_ggm_simulation_parallel(
   }
   pm.finish();
 
+  // See the OMRF runner: an interrupted worker leaves its remaining draws
+  // empty, which must not reach the caller as a result.
+  if (pm.shouldExit()) {
+    Rcpp::stop("Simulation interrupted by user.");
+  }
+
   Rcpp::List output(ndraws);
   for (int i = 0; i < ndraws; i++) {
     if (results[i].error) {
@@ -1132,6 +1165,12 @@ Rcpp::List run_mixed_simulation_parallel(
   }
   pm.finish();
 
+  // See the OMRF runner: an interrupted worker leaves its remaining draws
+  // empty, which must not reach the caller as a result.
+  if (pm.shouldExit()) {
+    Rcpp::stop("Simulation interrupted by user.");
+  }
+
   Rcpp::List output(ndraws);
   for (int i = 0; i < ndraws; i++) {
     if (results[i].error) {
diff --git a/tests/testthat/test-centrality.R b/tests/testthat/test-centrality.R
index 53b55b85..a3a49831 100644
--- a/tests/testthat/test-centrality.R
+++ b/tests/testthat/test-centrality.R
@@ -156,3 +156,85 @@ test_that("a centrality difference is summarized but not drawn", {
   on.exit(grDevices::dev.off(), add = TRUE)
   expect_invisible(plot(fit, type = "centrality", group = 2))
 })
+
+
+# ==============================================================================
+# Edge order. A mixed fit stores its pairwise draws by block --
+# discrete-discrete, then continuous-continuous, then cross -- which is not the
+# row-major upper triangle whenever the discrete and continuous columns
+# interleave. Summing the columns under the assumed order gives each node some
+# other node's edges, and the edge COUNT matches either way, so nothing shorter
+# than checking the incidence catches it.
+# ==============================================================================
+
+test_that("mixed-MRF strength sums each node's own edges", {
+  skip_on_cran()
+  fit = get_bgms_fit_mixed_mrf()
+  nodes = extract_arguments(fit)$data_columnnames
+  samples = extract_pairwise_interactions(fit)
+  strength = extract_centrality(fit)
+
+  expect_equal(colnames(strength), nodes)
+
+  # The stored column names carry the incidence; hand-sum from them.
+  pairs = do.call(rbind, strsplit(colnames(samples), "-", fixed = TRUE))
+  for(node in nodes) {
+    incident = pairs[, 1] == node | pairs[, 2] == node
+    expect_equal(
+      unname(strength[, node]),
+      unname(rowSums(abs(samples[, incident, drop = FALSE]))),
+      tolerance = 1e-12
+    )
+  }
+
+  # This fixture interleaves its discrete and continuous columns, so the block
+  # order really does differ from the row-major upper triangle here and the
+  # check above is not vacuous.
+  assumed = which(
+    upper.tri(matrix(0, length(nodes), length(nodes))),
+    arr.ind = TRUE
+  )
+  assumed = assumed[order(assumed[, "row"], assumed[, "col"]), , drop = FALSE]
+  expect_false(identical(
+    unname(bgms:::indicator_pair_index(fit, length(nodes))),
+    unname(assumed)
+  ))
+})
+
+test_that("ordinal, Blume-Capel and GGM centrality are bit-identical to before", {
+  skip_on_cran()
+  # Their draws ARE laid out as the row-major upper triangle, so reading the
+  # order off the fit has to reproduce the previous output exactly.
+  row_major_strength = function(samples, num_variables) {
+    pairs = which(
+      upper.tri(matrix(0, num_variables, num_variables)),
+      arr.ind = TRUE
+    )
+    pairs = pairs[order(pairs[, "row"], pairs[, "col"]), , drop = FALSE]
+    absolute = abs(samples)
+    vapply(
+      seq_len(num_variables),
+      function(v) {
+        rowSums(absolute[, pairs[, 1] == v | pairs[, 2] == v, drop = FALSE])
+      },
+      numeric(nrow(samples))
+    )
+  }
+
+  fits = list(
+    ordinal = get_bgms_fit_ordinal(),
+    blume_capel = get_bgms_fit_blumecapel(),
+    ggm = get_bgms_fit_ggm()
+  )
+  for(label in names(fits)) {
+    fit = fits[[label]]
+    nodes = extract_arguments(fit)$data_columnnames
+    expect_identical(
+      as.numeric(extract_centrality(fit)),
+      as.numeric(row_major_strength(
+        extract_pairwise_interactions(fit), length(nodes)
+      )),
+      info = label
+    )
+  }
+})
diff --git a/tests/testthat/test-input-validation.R b/tests/testthat/test-input-validation.R
index c140d463..495d3118 100644
--- a/tests/testthat/test-input-validation.R
+++ b/tests/testthat/test-input-validation.R
@@ -327,10 +327,11 @@ test_that("predict.bgms errors when newdata is missing", {
 
 test_that("predict.bgms errors on invalid variable names", {
   fit = get_bgms_fit()
-  args = extract_arguments(fit)
 
-  data("Wenchuan", package = "bgms")
-  newdata = Wenchuan[1:5, 1:args$num_variables]
+  # newdata has to be the fit's own variables: predict() now rejects a named
+  # newdata whose columns are not the model's, so borrowing another data set of
+  # the same width would stop on that instead of on the variable name.
+  newdata = get_prediction_data_binary(5)
 
   expect_error(
     predict(fit, newdata = newdata, variables = "NonexistentVar"),
@@ -340,10 +341,8 @@ test_that("predict.bgms errors on invalid variable names", {
 
 test_that("predict.bgms errors on out-of-range variable indices", {
   fit = get_bgms_fit()
-  args = extract_arguments(fit)
 
-  data("Wenchuan", package = "bgms")
-  newdata = Wenchuan[1:5, 1:args$num_variables]
+  newdata = get_prediction_data_binary(5)
 
   expect_error(
     predict(fit, newdata = newdata, variables = 999),
@@ -423,6 +422,8 @@ test_that("predict.bgmCompare errors on invalid variable names", {
   args = extract_arguments(fit)
 
   newdata = matrix(0L, nrow = 5, ncol = args$num_variables)
+  # Named, so predict() does not also warn that it is matching by position.
+  colnames(newdata) = args$data_columnnames
 
   expect_error(
     predict(fit, newdata = newdata, group = 1, variables = "NonexistentVar"),
diff --git a/tests/testthat/test-methods.R b/tests/testthat/test-methods.R
index de4e6dcc..45b03c08 100644
--- a/tests/testthat/test-methods.R
+++ b/tests/testthat/test-methods.R
@@ -636,8 +636,15 @@ test_that("predict.bgms GGM conditional mean matches analytic formula", {
 
   pred = predict(fit, newdata = newdata)
 
-  # Reconstruct the posterior mean precision matrix (precision = -2 * association)
-  omega_hat = extract_precision(fit)
+  # Reconstruct the posterior mean precision matrix the way predict() does:
+  # off-diagonals -2 * association, diagonal the mean of the raw diagonal draws.
+  # Not extract_precision(), which still takes 1 / posterior_mean_residual_
+  # variance for the diagonal -- a harmonic mean of the same draws, and so a
+  # different matrix.
+  omega_hat = bgms:::reconstruct_precision(
+    get_posterior_mean(fit, "pairwise"),
+    bgms:::posterior_mean_precision_diagonal(fit)
+  )
   p = args$num_variables
 
   # Center newdata on the training means (predict does the same internally)
diff --git a/tests/testthat/test-predict-contracts.R b/tests/testthat/test-predict-contracts.R
new file mode 100644
index 00000000..1eddfb25
--- /dev/null
+++ b/tests/testthat/test-predict-contracts.R
@@ -0,0 +1,349 @@
+# ==============================================================================
+# What the R layer around the prediction kernels owes its caller.
+#
+# The C++ kernels take the recoded observation matrix at face value and know
+# nothing about the fit's category scale, so four contracts live entirely in R:
+#
+#   - a row whose conditioning set is incomplete comes back as NA, not as a
+#     confident answer computed from R's NA_integer_ sentinel;
+#   - responses and probability labels are on the ORIGINAL category scale --
+#     the one simulate() returns and predict() reads back;
+#   - one row of newdata is a shape like any other;
+#   - the posterior-mean precision matrix is a posterior mean throughout,
+#     diagonal included;
+#   - newdata columns are matched to the fit deliberately, not by width alone.
+# ==============================================================================
+
+.contract_cache = new.env(parent = emptyenv())
+
+# Three recode maps, none of them the identity: an ordinal variable on 1..3, a
+# non-contiguous ordinal on {1, 5, 9}, and a Blume-Capel variable on 2..5.
+get_original_scale_fit = function() {
+  if(is.null(.contract_cache$fit)) {
+    set.seed(1)
+    n = 200
+    x = cbind(
+      a = sample(c(1, 2, 3), n, replace = TRUE),
+      b = sample(c(1, 5, 9), n, replace = TRUE),
+      bc = sample(2:5, n, replace = TRUE)
+    )
+    .contract_cache$fit = bgm(
+      x = x,
+      variable_type = c("ordinal", "ordinal", "blume-capel"),
+      baseline_category = c(1L, 1L, 3L),
+      iter = 200, warmup = 200, chains = 1, cores = 1,
+      update_method = "adaptive-metropolis",
+      display_progress = "none", seed = 7
+    )
+  }
+  .contract_cache$fit
+}
+
+collect_warnings = function(expr) {
+  seen = character()
+  withCallingHandlers(
+    force(expr),
+    warning = function(w) {
+      seen <<- c(seen, conditionMessage(w))
+      invokeRestart("muffleWarning")
+    }
+  )
+  seen
+}
+
+
+# ------------------------------------------------------------------------------
+# Original category scale (simulate <-> predict round trip)
+# ------------------------------------------------------------------------------
+
+test_that("simulate() and predict(type = 'response') meet on the original scale", {
+  skip_on_cran()
+  fit = get_original_scale_fit()
+
+  sim = simulate(fit, nsim = 30, method = "posterior-mean", seed = 3)
+  expect_true(all(sim[, "a"] %in% c(1, 2, 3)))
+  expect_true(all(sim[, "b"] %in% c(1, 5, 9)))
+  expect_true(all(sim[, "bc"] %in% 2:5))
+
+  resp = predict(fit, newdata = sim, type = "response")
+  expect_equal(colnames(resp), c("a", "b", "bc"))
+  expect_true(all(resp[, "a"] %in% c(1, 2, 3)))
+  expect_true(all(resp[, "b"] %in% c(1, 5, 9)))
+  expect_true(all(resp[, "bc"] %in% 2:5))
+
+  # The round trip closes: the responses are themselves valid newdata, with no
+  # unseen-category warning and no NA anywhere.
+  expect_no_warning(again <- predict(fit, newdata = resp, type = "response"))
+  expect_false(anyNA(again))
+})
+
+test_that("probability columns are labelled with the original categories", {
+  skip_on_cran()
+  fit = get_original_scale_fit()
+  sim = simulate(fit, nsim = 10, method = "posterior-mean", seed = 4)
+
+  probs = predict(fit, newdata = sim, type = "probabilities")
+  expect_equal(colnames(probs$a), paste0("cat_", c(1, 2, 3)))
+  expect_equal(colnames(probs$b), paste0("cat_", c(1, 5, 9)))
+  expect_equal(colnames(probs$bc), paste0("cat_", 2:5))
+
+  # The reported response is the value its own label names.
+  resp = predict(fit, newdata = sim, type = "response")
+  for(v in names(probs)) {
+    labels = as.numeric(sub("^cat_", "", colnames(probs[[v]])))
+    expect_equal(
+      unname(labels[max.col(probs[[v]], ties.method = "first")]),
+      unname(resp[, v])
+    )
+  }
+})
+
+test_that("posterior-sample prediction carries the original labels too", {
+  skip_on_cran()
+  fit = get_original_scale_fit()
+  sim = simulate(fit, nsim = 8, method = "posterior-mean", seed = 11)
+
+  probs = predict(
+    fit,
+    newdata = sim, type = "probabilities",
+    method = "posterior-sample", ndraws = 5, seed = 2
+  )
+  expect_equal(colnames(probs$b), paste0("cat_", c(1, 5, 9)))
+  expect_equal(colnames(attr(probs, "sd")$b), paste0("cat_", c(1, 5, 9)))
+})
+
+
+# ------------------------------------------------------------------------------
+# Incomplete conditioning
+# ------------------------------------------------------------------------------
+
+test_that("an unseen category or a plain NA blanks the row's other variables", {
+  skip_on_cran()
+  fit = get_original_scale_fit()
+  sim = simulate(fit, nsim = 6, method = "posterior-mean", seed = 5)
+
+  newdata = sim
+  newdata[1, "a"] = 99 # a value never observed in training
+  newdata[2, "b"] = NA # a plain NA
+
+  probs = suppressWarnings(
+    predict(fit, newdata = newdata, type = "probabilities")
+  )
+  # The other variables in a corrupted row are undefined ...
+  expect_true(all(is.na(probs$b[1, ])))
+  expect_true(all(is.na(probs$bc[1, ])))
+  expect_true(all(is.na(probs$a[2, ])))
+  expect_true(all(is.na(probs$bc[2, ])))
+  # ... while the variable whose own value is missing is not: its conditional
+  # never reads that value. Clean rows are untouched.
+  expect_false(anyNA(probs$a[1, ]))
+  expect_false(anyNA(probs$b[2, ]))
+  expect_false(anyNA(probs$a[3:6, ]))
+  expect_false(anyNA(probs$bc[3:6, ]))
+
+  resp = suppressWarnings(predict(fit, newdata = newdata, type = "response"))
+  expect_true(all(is.na(resp[1, c("b", "bc")])))
+  expect_true(all(is.na(resp[2, c("a", "bc")])))
+  expect_false(anyNA(resp[3:6, ]))
+})
+
+test_that("incomplete conditioning warns once, with the affected row count", {
+  skip_on_cran()
+  fit = get_original_scale_fit()
+  sim = simulate(fit, nsim = 6, method = "posterior-mean", seed = 5)
+
+  newdata = sim
+  newdata[1, "a"] = 99
+  newdata[2, "b"] = NA
+
+  seen = collect_warnings(predict(fit, newdata = newdata))
+  hits = grep("conditioning variable is missing", seen, fixed = TRUE)
+  expect_length(hits, 1L)
+  expect_match(seen[hits], "2 row\\(s\\)")
+
+  # A plain NA on its own used to pass in silence.
+  only_na = sim
+  only_na[4, "bc"] = NA
+  seen_na = collect_warnings(predict(fit, newdata = only_na))
+  expect_length(grep("conditioning variable is missing", seen_na, fixed = TRUE), 1L)
+  expect_match(seen_na, "1 row\\(s\\)", all = FALSE)
+})
+
+test_that("a mixed fit blanks continuous predictions on a missing discrete cell", {
+  skip_on_cran()
+  fit = get_bgms_fit_mixed_mrf()
+  newdata = get_prediction_data_mixed(n = 5)
+  newdata[1, "d1"] = NA
+
+  probs = suppressWarnings(
+    predict(fit, newdata = newdata, type = "probabilities")
+  )
+  # The discrete rest score is what the sentinel corrupts, and it feeds the
+  # continuous conditional as well as the discrete ones.
+  expect_true(all(is.na(probs$c1[1, ])))
+  expect_true(all(is.na(probs$d2[1, ])))
+  expect_false(anyNA(probs$d1[1, ]))
+  expect_false(anyNA(probs$c1[2:5, ]))
+
+  resp = suppressWarnings(predict(fit, newdata = newdata, type = "response"))
+  expect_true(all(is.na(resp[1, c("c1", "d2", "c2", "d3")])))
+  expect_false(anyNA(resp[2:5, ]))
+})
+
+
+# ------------------------------------------------------------------------------
+# One row of newdata
+# ------------------------------------------------------------------------------
+
+test_that("a single row of newdata predicts for every model family", {
+  skip_on_cran()
+
+  single_row_ok = function(fit, newdata, ...) {
+    arguments = extract_arguments(fit)
+    resp = predict(fit, newdata = newdata, type = "response", ...)
+    expect_true(is.matrix(resp))
+    expect_equal(dim(resp), c(1L, arguments$num_variables))
+    expect_equal(colnames(resp), arguments$data_columnnames)
+    expect_false(anyNA(resp))
+  }
+
+  single_row_ok(get_bgms_fit(), get_prediction_data_binary(n = 1))
+  single_row_ok(get_bgms_fit_ggm(), get_prediction_data_ggm(n = 1))
+  single_row_ok(get_bgms_fit_mixed_mrf(), get_prediction_data_mixed(n = 1))
+  single_row_ok(
+    get_bgmcompare_fit(), get_prediction_data_bgmcompare_binary(n = 1),
+    group = 1
+  )
+})
+
+
+# ------------------------------------------------------------------------------
+# The posterior-mean precision diagonal
+# ------------------------------------------------------------------------------
+
+test_that("the GGM precision diagonal is the pooled mean of the raw draws", {
+  skip_on_cran()
+  fit = get_bgms_fit_ggm()
+
+  raw_diagonal = do.call(rbind, get_raw_samples(fit)$main)
+  omega = bgms:::reconstruct_precision(
+    get_posterior_mean(fit, "pairwise"),
+    bgms:::posterior_mean_precision_diagonal(fit)
+  )
+  expect_equal(
+    unname(diag(omega)), unname(colMeans(raw_diagonal)),
+    tolerance = 1e-12
+  )
+
+  # Not the reciprocal of the residual variance, which is a harmonic mean of
+  # the same draws and so a different number.
+  expect_false(isTRUE(all.equal(
+    unname(diag(omega)),
+    unname(1 / get_posterior_mean(fit, "residual_variance")),
+    tolerance = 1e-8
+  )))
+
+  # Omega is now exactly the mean of the per-draw precision matrices, and a
+  # mean of positive-definite matrices is positive definite.
+  raw_pairwise = do.call(rbind, get_raw_samples(fit)$pairwise)
+  p = extract_arguments(fit)$num_variables
+  pooled = matrix(0, p, p)
+  for(i in seq_len(nrow(raw_pairwise))) {
+    pooled = pooled + bgms:::build_precision_from_draw(
+      raw_pairwise[i, ], raw_diagonal[i, ], p
+    )
+  }
+  pooled = pooled / nrow(raw_pairwise)
+  expect_equal(unname(omega), unname(pooled), tolerance = 1e-12)
+  expect_true(all(eigen(omega, symmetric = TRUE, only.values = TRUE)$values > 0))
+})
+
+test_that("the mixed continuous diagonal comes from its own draws", {
+  skip_on_cran()
+  fit = get_bgms_fit_mixed_mrf()
+  arguments = extract_arguments(fit)
+
+  params = bgms:::build_mixed_params_mean(fit, arguments)
+  draws = bgms:::mixed_cont_diagonal_draws(fit, arguments)
+  expect_equal(
+    unname(diag(params$pairwise_cont)), unname(colMeans(draws)),
+    tolerance = 1e-12
+  )
+
+  precision = -2 * params$pairwise_cont
+  expect_true(
+    all(eigen(precision, symmetric = TRUE, only.values = TRUE)$values > 0)
+  )
+})
+
+
+# ------------------------------------------------------------------------------
+# newdata columns
+# ------------------------------------------------------------------------------
+
+test_that("named newdata must carry the fit's variables, in its order", {
+  fit = get_bgms_fit()
+  nodes = extract_arguments(fit)$data_columnnames
+  newdata = as.matrix(get_prediction_data_binary(n = 4))
+
+  renamed = newdata
+  colnames(renamed) = paste0("x", seq_along(nodes))
+  expect_error(
+    predict(fit, newdata = renamed),
+    "column names do not match the fitted model"
+  )
+
+  reordered = newdata[, rev(seq_along(nodes)), drop = FALSE]
+  expect_error(predict(fit, newdata = reordered), "in a different order")
+
+  expect_no_error(predict(fit, newdata = newdata))
+})
+
+test_that("unnamed newdata is accepted by position, and says so", {
+  fit = get_bgms_fit()
+  newdata = as.matrix(get_prediction_data_binary(n = 4))
+  dimnames(newdata) = NULL
+
+  expect_warning(
+    predict(fit, newdata = newdata),
+    "matched to the fitted model by position"
+  )
+})
+
+test_that("bgmCompare applies the same newdata column rules", {
+  fit = get_bgmcompare_fit()
+  nodes = extract_arguments(fit)$data_columnnames
+  newdata = as.matrix(get_prediction_data_bgmcompare_binary(n = 4))
+
+  renamed = newdata
+  colnames(renamed) = paste0("x", seq_along(nodes))
+  expect_error(
+    predict(fit, newdata = renamed, group = 1),
+    "column names do not match the fitted model"
+  )
+
+  dimnames(newdata) = NULL
+  expect_warning(
+    predict(fit, newdata = newdata, group = 1),
+    "matched to the fitted model by position"
+  )
+})
+
+
+# ------------------------------------------------------------------------------
+# simulate() entry validation
+# ------------------------------------------------------------------------------
+
+test_that("simulate.bgms validates nsim and iter for every path", {
+  fit = get_bgms_fit()
+  expect_error(simulate(fit, nsim = 0), "nsim")
+  expect_error(simulate(fit, nsim = 10.5), "nsim")
+  expect_error(simulate(fit, nsim = 10, iter = -1), "iter")
+})
+
+test_that("simulate.bgms survives a machine that cannot report its cores", {
+  # parallel::detectCores() is documented to return NA when it cannot tell, and
+  # that is the argument's default.
+  fit = get_bgms_fit()
+  expect_no_error(simulate(fit, nsim = 5, iter = 10, cores = NA_integer_))
+})
diff --git a/tests/testthat/test-predict-recode-map.R b/tests/testthat/test-predict-recode-map.R
index e86cd3a7..34f8871d 100644
--- a/tests/testthat/test-predict-recode-map.R
+++ b/tests/testthat/test-predict-recode-map.R
@@ -34,16 +34,21 @@ test_that("values not observed in training give NA with a warning", {
   expect_equal(out[1, 1], 0)
 })
 
-test_that("Blume-Capel columns (NULL level entry) fall back to subtract-min", {
-  # category_levels present but NULL for this variable -> legacy behaviour.
+test_that("an ordinal column with a NULL level entry is an error, not a guess", {
+  # The old fallback shifted by the newdata column minimum, which is the
+  # training offset only by coincidence. Every current fit carries a map.
   levels = list(NULL)
-  out = recode(matrix(c(2, 3, 4), ncol = 1), TRUE, category_levels = levels)
-  expect_equal(out[, 1], c(0, 1, 2))
+  expect_error(
+    recode(matrix(c(2, 3, 4), ncol = 1), TRUE, category_levels = levels),
+    "predates the recode map"
+  )
 })
 
-test_that("no map at all (old fit) preserves legacy subtract-min behaviour", {
-  out = recode(matrix(c(2, 3, 4), ncol = 1), TRUE, category_levels = NULL)
-  expect_equal(out[, 1], c(0, 1, 2))
+test_that("no map at all (old fit) is an error naming the refit", {
+  expect_error(
+    recode(matrix(c(2, 3, 4), ncol = 1), TRUE, category_levels = NULL),
+    "predates the recode map"
+  )
 })
 
 test_that("continuous/non-ordinal columns are left unchanged", {
diff --git a/tests/testthat/test-simulate-mrf.R b/tests/testthat/test-simulate-mrf.R
index 1557f1a5..f0563baf 100644
--- a/tests/testthat/test-simulate-mrf.R
+++ b/tests/testthat/test-simulate-mrf.R
@@ -460,3 +460,74 @@ test_that("simulate_mrf rejects mixed continuous and ordinal", {
     "all variables must be of type"
   )
 })
+
+
+# ------------------------------------------------------------------------------
+# Numerical stability of the full conditional
+# ------------------------------------------------------------------------------
+# The full conditional subtracts the largest category exponent before
+# exponentiating. Without that, a large rest score sends exp() to Inf, the
+# cumulative weights become Inf or NaN, and the category search degenerates.
+
+test_that("large parameters do not overflow the full conditional", {
+  num_categories = c(2L, 2L, 3L)
+  pairwise = matrix(0, 3, 3)
+  pairwise[lower.tri(pairwise)] = c(80, -60, 50)
+  pairwise = pairwise + t(pairwise)
+  main = matrix(NA_real_, 3, max(num_categories))
+  main[1, 1:2] = c(60, -40)
+  main[2, 1:2] = c(-100, 80)
+  main[3, 1:3] = c(20, 40, -80)
+
+  result = simulate_mrf(
+    num_states = 40, num_variables = 3, num_categories = num_categories,
+    pairwise = pairwise, main = main,
+    variable_type = rep("ordinal", 3), baseline_category = rep(0L, 3),
+    iter = 50, seed = 3
+  )
+
+  expect_true(all(is.finite(result)))
+  expect_true(all(result >= 0))
+  expect_true(all(result <= rep(num_categories, each = nrow(result))))
+})
+
+test_that("the stabilized ordinal conditional samples the intended law", {
+  skip_on_cran()
+  # Three variables small enough to enumerate the joint exactly, so the
+  # simulated marginals have something other than themselves to agree with.
+  num_categories = c(2L, 2L, 3L)
+  pairwise = matrix(0, 3, 3)
+  pairwise[lower.tri(pairwise)] = c(0.4, -0.3, 0.25)
+  pairwise = pairwise + t(pairwise)
+  main = matrix(NA_real_, 3, max(num_categories))
+  main[1, 1:2] = c(0.3, -0.2)
+  main[2, 1:2] = c(-0.5, 0.4)
+  main[3, 1:3] = c(0.1, 0.2, -0.4)
+
+  simulated = simulate_mrf(
+    num_states = 50000, num_variables = 3, num_categories = num_categories,
+    pairwise = pairwise, main = main,
+    variable_type = rep("ordinal", 3), baseline_category = rep(0L, 3),
+    iter = 100, seed = 99
+  )
+
+  grid = as.matrix(expand.grid(lapply(num_categories, function(k) 0:k)))
+  log_weight = apply(grid, 1, function(z) {
+    total = 0
+    for(v in seq_len(3)) if(z[v] > 0) total = total + main[v, z[v]]
+    for(i in 1:2) {
+      for(j in (i + 1):3) total = total + 2 * z[i] * z[j] * pairwise[i, j]
+    }
+    total
+  })
+  weight = exp(log_weight - max(log_weight))
+  weight = weight / sum(weight)
+
+  for(v in seq_len(3)) {
+    exact = as.numeric(tapply(weight, grid[, v], sum))
+    empirical = as.numeric(
+      table(factor(simulated[, v], levels = 0:num_categories[v]))
+    ) / nrow(simulated)
+    expect_lt(max(abs(empirical - exact)), 0.01)
+  }
+})
diff --git a/tests/testthat/test-simulate-predict-regression.R b/tests/testthat/test-simulate-predict-regression.R
index 423e08fa..57f523d5 100644
--- a/tests/testthat/test-simulate-predict-regression.R
+++ b/tests/testthat/test-simulate-predict-regression.R
@@ -697,10 +697,14 @@ test_that("mixed predict recodes discrete newdata to the original scale", {
   probs_shift = predict(fit_shift, newdata = x_shift[1:20, ], type = "probabilities")
   probs_base = predict(fit_base, newdata = x_base[1:20, ], type = "probabilities")
   for(nm in names(probs_shift)) {
-    expect_equal(probs_shift[[nm]], probs_base[[nm]],
+    # Values, not labels: the labels name the ORIGINAL categories, so the two
+    # relabelings are supposed to differ there and only there.
+    expect_equal(unname(probs_shift[[nm]]), unname(probs_base[[nm]]),
       info = sprintf("relabel invariance %s", nm)
     )
   }
+  expect_equal(colnames(probs_shift$d), paste0("cat_", 1:3))
+  expect_equal(colnames(probs_base$d), paste0("cat_", 0:2))
 
   # simulate() returns the original discrete scale and the round trip is clean.
   sim = simulate(fit_shift, nsim = 40, method = "posterior-mean", seed = 5)
@@ -771,9 +775,13 @@ test_that("sparse category codings recode to the fitted categories", {
   probs_sparse = predict(fit_sparse, newdata = x_sparse[1:20, ], type = "probabilities")
   probs_dense = predict(fit_dense, newdata = x_dense[1:20, ], type = "probabilities")
   for(nm in names(probs_sparse)) {
-    expect_equal(probs_sparse[[nm]], probs_dense[[nm]],
+    # Values, not labels: the labels name the ORIGINAL categories, so the
+    # sparse fit's read cat_1/2/4/5 against the dense fit's cat_0..3.
+    expect_equal(unname(probs_sparse[[nm]]), unname(probs_dense[[nm]]),
       info = sprintf("sparse/dense relabel invariance %s", nm)
     )
+    expect_equal(colnames(probs_sparse[[nm]]), paste0("cat_", sparse))
+    expect_equal(colnames(probs_dense[[nm]]), paste0("cat_", 0:3))
   }
 
   # --- simulate() -> predict() round trip stays on the sparse original scale ---