fix(residual): precision diagonal, mixed anchor scale, configure fail-fast

R/anchor_curve.R

@@ -33,6 +33,11 @@ anchor_draws = function(fit) {
   if(identical(spec$model_type, "compare")) {
     return(compare_anchor_draws(fit, spec, raw))
   }
+  # GGM pairwise draws are K_ij; the slab prior is evaluated at -K_ij / 2
+  # (ggm_model.cpp:313). The mixed sampler stores every pairwise block on the
+  # slab's own scale: the discrete and cross blocks as is
+  # (mixed_mrf_gradient.cpp:503, 521), the continuous block as -Omega / 2
+  # (mixed_mrf_model.cpp:534, mixed_mrf_gradient.cpp:574).
   slab_factor = if(identical(spec$model_type, "ggm")) -0.5 else 1
   theta = lapply(raw$pairwise, function(m) slab_factor * m)
   gamma = raw$indicator
@@ -54,8 +59,14 @@ anchor_draws = function(fit) {
 # Under vary = "slab-and-diagonal" the refits set the diagonal rate to eta / s
 # (vary_diagonal_rate()), so the weight ratio carries the diagonal prior as
 # well as the slab. The sampler evaluates that prior at K_jj / 2
-# (ggm_model.cpp:320, mixed_mrf_gradient.cpp:559); the raw main draws store
-# K_jj.
+# (ggm_model.cpp:320, mixed_mrf_gradient.cpp:559). The two model types store
+# the diagonal on different scales:
+#
+#   ggm        raw main draws are K_jj (build_precision_from_draw());
+#              K_jj / 2 = 0.5 * draw.
+#   mixed_mrf  the "(precision diag)" columns are the diagonal of
+#              pairwise_effects_continuous_, stored as -Omega / 2 for the whole
+#              continuous block (mixed_mrf_model.cpp:534); K_jj / 2 = -draw.
 #
 # @param spec  The fit's spec.
 # @param raw   Its raw samples.
@@ -64,7 +75,8 @@ anchor_draws = function(fit) {
 #   diagonal elements, shape = the Gamma shape alpha), or NULL.
 # ------------------------------------------------------------------
 anchor_diagonal = function(spec, raw) {
-  cols = if(identical(spec$model_type, "ggm")) {
+  is_ggm = identical(spec$model_type, "ggm")
+  cols = if(is_ggm) {
     seq_len(ncol(raw$main[[1]]))
   } else if(identical(spec$model_type, "mixed_mrf")) {
     which(endsWith(raw$parameter_names[["main"]], "(precision diag)"))
@@ -74,8 +86,11 @@ anchor_diagonal = function(spec, raw) {
   if(length(cols) == 0L) {
     return(NULL)
   }
+  to_half_precision = if(is_ggm) 0.5 else -1
   list(
-    sum = lapply(raw$main, function(m) 0.5 * rowSums(m[, cols, drop = FALSE])),
+    sum = lapply(raw$main, function(m) {
+      to_half_precision * rowSums(m[, cols, drop = FALSE])
+    }),
     n = length(cols),
     shape = spec$prior$scale_shape %||% 1
   )

R/extract_prior_inclusion_probabilities.R

@@ -177,6 +177,7 @@ prior_only_chain_pips = function(spec, iter, warmup) {
     )
     conjugate_slab = p$interaction_prior_type %in% c("cauchy", "normal") &&
       p$scale_prior_type %in% c("gamma", "exponential")
+    chain_update_method = if(conjugate_slab) "gibbs" else "adaptive-metropolis"
     results = sample_ggm(
       inputFromR = list(
         n = 0L,
@@ -195,7 +196,10 @@ prior_only_chain_pips = function(spec, iter, warmup) {
       no_warmup = as.integer(warmup),
       no_chains = 1L,
       edge_selection = TRUE,
-      sampler_type = if(conjugate_slab) "gibbs" else "adaptive-metropolis",
+      sampler_type = chain_update_method,
+      # The C++ default is the NUTS target 0.80; use the target bgm() resolves
+      # for this update method.
+      target_acceptance = resolve_target_acceptance(chain_update_method),
       seed = as.integer(s$seed + 1L),
       no_threads = 1L,
       progress_type = 0L,
@@ -276,6 +280,7 @@ prior_only_chain_pips = function(spec, iter, warmup) {
     dirichlet_alpha = p$dirichlet_alpha,
     lambda = p$lambda,
     sampler_type = "adaptive-metropolis",
+    target_acceptance = resolve_target_acceptance("adaptive-metropolis"),
     delta = p$delta,
     edge_prior_correction = correction
   )

R/extractor_functions.R

@@ -1707,6 +1707,9 @@ extract_ess.bgmCompare = function(bgms_object, estimator = c("rb", "mixt")) {
 #' For mixed MRF models the precision matrix is reconstructed from the
 #' internal association-scale parameterization.
 #'
+#' Each entry, diagonal included, is the posterior mean of the corresponding
+#' precision draw. This is the matrix `predict()` and `simulate()` use.
+#'
 #' @param bgms_object A fitted model object of class `bgms` (from [bgm()]).
 #'
 #' @return A named numeric matrix containing the posterior mean precision
@@ -1745,23 +1748,33 @@ extract_precision.bgms = function(bgms_object) {
     return(invisible(NULL))
   }
 
-  rv = get_posterior_mean(bgms_object, "residual_variance")
   associations = get_posterior_mean(bgms_object, "pairwise")
 
+  # Diagonal: E[K_jj], the mean of the raw diagonal draws, as in predict() and
+  # simulate().
+
   if(isTRUE(arguments$is_mixed)) {
-    # Mixed MRF: extract the q x q continuous block, convert to precision
+    # Mixed MRF: the q x q continuous block, converted to precision. The sampler
+    # stores this block, diagonal included, as -K / 2.
     cont_idx = arguments$continuous_indices
     cont_names = arguments$data_columnnames_continuous
     cont_block = associations[cont_idx, cont_idx]
-    precision = -2 * cont_block
-    diag(precision) = 1 / rv
+    # Column means one at a time, as in build_mixed_params_mean(), for a
+    # bit-identical match with predict().
+    diag_draws = mixed_cont_diagonal_draws(bgms_object, arguments)
+    precision = reconstruct_precision(
+      cont_block,
+      -2 * apply(diag_draws, 2, mean)
+    )
     dimnames(precision) = list(cont_names, cont_names)
     return(precision)
   }
 
-  # GGM: associations are stored at half precision scale; convert to precision
-  precision = -2 * associations
-  diag(precision) = 1 / rv
+  # GGM: associations are stored as -K_ij / 2; the raw diagonal draws as K_jj.
+  precision = reconstruct_precision(
+    associations,
+    posterior_mean_precision_diagonal(bgms_object)
+  )
   return(precision)
 }
 

R/sample_ggm_prior.R

@@ -381,14 +381,8 @@ sample_ggm_prior = function(
     no_chains = 1L,
     edge_selection = TRUE,
     sampler_type = update_method,
-    # The target bgm() resolves for this update method (validate_sampler()):
-    # 0.44, the componentwise random-walk MH optimum, for the proposal-SD
-    # tuning of adaptive-metropolis; NA for gibbs, which tunes nothing. The
-    # C++ default of 0.80 is the NUTS target.
-    target_acceptance = switch(update_method,
-      "adaptive-metropolis" = 0.44,
-      "gibbs"               = NA_real_
-    ),
+    # The target bgm() resolves for this update method (validate_sampler()).
+    target_acceptance = resolve_target_acceptance(update_method),
     seed = as.integer(seed),
     no_threads = 1L,
     progress_type = if(verbose) 2L else 0L,

R/validate_sampler.R

@@ -21,6 +21,34 @@
 # check, and never exceed the machine. The sweep builders add a no-forking rule
 # on Windows on top of it (normalize_builder_cores()); the sampler's chain
 # workers are threads and need no such rule.
+# ------------------------------------------------------------------
+# resolve_target_acceptance (internal)
+# ------------------------------------------------------------------
+# The default acceptance target for an update method. Every chain the package
+# starts (the bgm() fit, the reference chain in sample_ggm_prior(), the
+# prior-only chain in prior_only_chain_pips()) takes its target from here, so
+# they all tune the same way for a given method. The C++ default is 0.80, the
+# NUTS target.
+#
+# @param update_method "adaptive-metropolis", "nuts", or "gibbs".
+#
+# Returns: numeric scalar.
+# ------------------------------------------------------------------
+resolve_target_acceptance = function(update_method) {
+  switch(update_method,
+    # The componentwise random-walk MH optimum the AM chain tunes its
+    # between-model proposal SDs to.
+    "adaptive-metropolis" = 0.44,
+    "nuts"                = 0.80,
+    # Exact draw: no acceptance target. NA_real_ records that honestly (not a
+    # fake 0.44); it stays numeric length-1 for the downstream contract and is
+    # unused by the Gibbs path (the C++ side ignores target_accept for gibbs).
+    "gibbs"               = NA_real_,
+    stop("Unknown update_method: ", update_method)
+  )
+}
+
+
 check_limit_cores = function() {
   check_limit = Sys.getenv("_R_CHECK_LIMIT_CORES_", "")
   nzchar(check_limit) && !identical(tolower(check_limit), "false")
@@ -154,14 +182,7 @@ validate_sampler = function(update_method,
     target_accept = min(target_accept, 1 - sqrt(.Machine$double.eps))
     target_accept = max(target_accept, 0 + sqrt(.Machine$double.eps))
   } else {
-    target_accept = switch(update_method,
-      "adaptive-metropolis" = 0.44,
-      "nuts"                = 0.80,
-      # Exact draw: no acceptance target. NA_real_ records that honestly (not a
-      # fake 0.44); it stays numeric length-1 for the downstream contract and is
-      # unused by the Gibbs path (the C++ side ignores target_accept for gibbs).
-      "gibbs"               = NA_real_
-    )
+    target_accept = resolve_target_acceptance(update_method)
   }
 
   # --- iter / warmup ----------------------------------------------------------

configure

@@ -1,13 +1,25 @@
 #!/bin/sh
 
+# Exit on the first failing command (for example an Rscript call with R_HOME
+# unset), instead of handing make an empty sources.mk.
+set -e
+
+RSCRIPT="${R_HOME}/bin/Rscript"
+
 # RcppParallel flags
-RCPP_PARALLEL_CPPFLAGS=`"${R_HOME}/bin/Rscript" -e "cat(RcppParallel::CxxFlags())"`
-RCPP_PARALLEL_LIBS=`"${R_HOME}/bin/Rscript" -e "cat(RcppParallel::LdFlags())"`
+RCPP_PARALLEL_CPPFLAGS=`"${RSCRIPT}" -e "cat(RcppParallel::CxxFlags())"`
+RCPP_PARALLEL_LIBS=`"${RSCRIPT}" -e "cat(RcppParallel::LdFlags())"`
 
 # Generate sources.mk using R (the script writes src/sources.mk itself)
-"${R_HOME}/bin/Rscript" inst/generate_makevars_sources.R
+"${RSCRIPT}" inst/generate_makevars_sources.R
 
 # Substitute into Makevars
 sed -e "s|@RCPP_PARALLEL_CPPFLAGS@|${RCPP_PARALLEL_CPPFLAGS}|" \
     -e "s|@RCPP_PARALLEL_LIBS@|${RCPP_PARALLEL_LIBS}|" \
 src/Makevars.in > src/Makevars
+
+# make needs a non-empty object list; check the file, not only the exit status.
+if [ ! -s src/sources.mk ]; then
+  echo "configure: src/sources.mk was not generated; cannot build bgms." >&2
+  exit 1
+fi

configure.win

@@ -1,13 +1,24 @@
 #!/bin/sh
 
+# Exit on the first failing command; see ./configure.
+set -e
+
+RSCRIPT="${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe"
+
 # RcppParallel flags
-RCPP_PARALLEL_CPPFLAGS=`"${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "cat(RcppParallel::CxxFlags())"`
-RCPP_PARALLEL_LIBS=`"${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" -e "cat(RcppParallel::LdFlags())"`
+RCPP_PARALLEL_CPPFLAGS=`"${RSCRIPT}" -e "cat(RcppParallel::CxxFlags())"`
+RCPP_PARALLEL_LIBS=`"${RSCRIPT}" -e "cat(RcppParallel::LdFlags())"`
 
 # Generate sources.mk using R (the script writes src/sources.mk itself)
-"${R_HOME}/bin${R_ARCH_BIN}/Rscript.exe" inst/generate_makevars_sources.R
+"${RSCRIPT}" inst/generate_makevars_sources.R
 
 # Substitute into Makevars.win
 sed -e "s|@RCPP_PARALLEL_CPPFLAGS@|${RCPP_PARALLEL_CPPFLAGS}|" \
     -e "s|@RCPP_PARALLEL_LIBS@|${RCPP_PARALLEL_LIBS}|" \
     src/Makevars.in > src/Makevars.win
+
+# make needs a non-empty object list; check the file, not only the exit status.
+if [ ! -s src/sources.mk ]; then
+  echo "configure.win: src/sources.mk was not generated; cannot build bgms." >&2
+  exit 1
+fi

tests/testthat/test-anchor-curve.R

@@ -199,17 +199,20 @@ test_that("anchor_diagonal reads the precision diagonal a model actually has", {
   expect_equal(d$sum[[1]], 0.5 * rowSums(main[[1]]))
   expect_equal(d$sum[[2]], 0.5 * rowSums(main[[2]]))
 
-  # Mixed MRF: only the columns that are precision diagonals.
+  # Mixed MRF: only the "(precision diag)" columns count, stored as -K_jj / 2.
+  kjj = matrix(seq(1, 10, length.out = 10), 5, 2)
+  stored = -kjj / 2
   raw_mixed = list(
-    main = list(matrix(1:20, 5, 4)),
+    main = list(cbind(matrix(0, 5, 2), stored)),
     parameter_names = list(main = c(
       "a (1)", "b (mean)", "b (precision diag)", "c (precision diag)"
     ))
   )
   spec_mixed = list(model_type = "mixed_mrf", prior = list(scale_shape = 1))
   dm = anchor_diagonal(spec_mixed, raw_mixed)
   expect_equal(dm$n, 2L)
-  expect_equal(dm$sum[[1]], 0.5 * rowSums(matrix(1:20, 5, 4)[, 3:4]))
+  expect_equal(dm$sum[[1]], rowSums(kjj) / 2)
+  expect_true(all(dm$sum[[1]] > 0))
 
   # A discrete model has no precision diagonal at all.
   expect_null(anchor_diagonal(
@@ -219,6 +222,57 @@ test_that("anchor_diagonal reads the precision diagonal a model actually has", {
 })
 
 
+test_that("the anchor draws of a real fit are in the frames the priors use", {
+  # Pin both scales to something the fit reports independently: the residual
+  # variance is 1 / K_jj averaged per draw, and the posterior mean pairwise
+  # matrix is on the association scale the slab is evaluated on.
+  for(spec in list(
+    list(label = "ggm", fit = get_bgms_fit_ggm()),
+    list(label = "mixed", fit = get_bgms_fit_mixed_mrf())
+  )) {
+    fit = spec$fit
+    ctx = spec$label
+    draws = anchor_draws(fit)
+
+    # Diagonal statistic: sum_j K_jj / 2, positive, and reproducing the
+    # posterior mean residual variance as the mean of 1 / K_jj.
+    kjj_half = unlist(draws$diagonal$sum)
+    expect_true(all(kjj_half > 0), info = ctx)
+
+    raw = get_raw_samples(fit)
+    cols = if(identical(ctx, "ggm")) {
+      seq_len(ncol(raw$main[[1]]))
+    } else {
+      which(endsWith(raw$parameter_names[["main"]], "(precision diag)"))
+    }
+    scale = if(identical(ctx, "ggm")) 1 else -2
+    kjj = do.call(rbind, lapply(raw$main, function(m) {
+      scale * m[, cols, drop = FALSE]
+    }))
+    expect_equal(
+      colMeans(1 / kjj),
+      unname(get_posterior_mean(fit, "residual_variance")),
+      info = ctx
+    )
+    expect_equal(
+      kjj_half, unname(rowSums(kjj) / 2),
+      info = ctx
+    )
+
+    # Slab frame: the association scale the posterior mean pairwise matrix is
+    # reported on. The GGM stores raw precision and needs the -0.5 factor; the
+    # mixed sampler already stores every block in its slab's own frame.
+    theta_means = colMeans(do.call(rbind, draws$theta))
+    pm = get_posterior_mean(fit, "pairwise")
+    expect_equal(
+      sort(round(theta_means, 10)),
+      sort(round(pm[upper.tri(pm)], 10)),
+      info = ctx
+    )
+  }
+})
+
+
 test_that("a slab family the reweighting identity is not written for stops", {
   # Only the normal and Cauchy identities are written; any other family stops.
   theta = matrix(rnorm(20), 10, 2)