fix(prior-sampler): acceptance targets, honest family reporting, dead knobs

NEWS.md

@@ -41,6 +41,12 @@ released from this line yet.
   follows it, as the measurement that says whether the flag is noise rather than
   as the first thing to try. `vignette("diagnostics")` gives the same order.
 
+## Deprecations
+
+* The `step_size` argument of `sample_ggm_prior()` is deprecated and ignored:
+  the NUTS warmup sets its own initial step size. The returned `step_size`
+  field has been removed.
+
 ## Documentation
 
 * New vignette, `vignette("fast-ggm")`, on the Gaussian graphical model route

R/RcppExports.R

@@ -45,8 +45,8 @@ ggm_test_forward_map <- function(theta, edge_indicators) {
     .Call(`_bgms_ggm_test_forward_map`, theta, edge_indicators)
 }
 
-sample_ggm_prior_cpp <- function(p, n_samples, n_warmup = 1000L, pairwise_scale = 2.5, interaction_prior_type = "cauchy", scale_prior_type = "gamma", gamma_shape = 1.0, gamma_rate = 1.0, step_size = 0.1, max_depth = 10L, seed = 1L, verbose = TRUE, edge_indicators_nullable = NULL, delta = 0.0) {
-    .Call(`_bgms_sample_ggm_prior`, p, n_samples, n_warmup, pairwise_scale, interaction_prior_type, scale_prior_type, gamma_shape, gamma_rate, step_size, max_depth, seed, verbose, edge_indicators_nullable, delta)
+sample_ggm_prior_cpp <- function(p, n_samples, n_warmup = 1000L, pairwise_scale = 2.5, interaction_prior_type = "cauchy", scale_prior_type = "gamma", gamma_shape = 1.0, gamma_rate = 1.0, max_depth = 10L, seed = 1L, verbose = TRUE, edge_indicators_nullable = NULL, delta = 0.0) {
+    .Call(`_bgms_sample_ggm_prior`, p, n_samples, n_warmup, pairwise_scale, interaction_prior_type, scale_prior_type, gamma_shape, gamma_rate, max_depth, seed, verbose, edge_indicators_nullable, delta)
 }
 
 .compute_ess_cpp <- function(array3d) {

R/sample_ggm_prior.R

@@ -67,10 +67,8 @@
 #'   \code{s}). Default: \code{exponential_prior(eta = 1)}; with the
 #'   default \code{normal_prior(scale = 1)} interaction prior this
 #'   resolves to \eqn{K_{ii}/2 \sim \textrm{Exponential}(1)}.
-#' @param step_size Positive numeric. Initial NUTS step size used to seed
-#'   dual-averaging adaptation. Default \code{0.1}. Used only for
-#'   \code{spec = "conditional"} (NUTS path); ignored for the
-#'   \code{"joint"} MH path.
+#' @param step_size `r lifecycle::badge("deprecated")` Deprecated and ignored.
+#'   The NUTS warmup determines its own initial step size.
 #' @param max_depth Integer. Maximum NUTS tree depth. Default \code{10}.
 #'   Used only for \code{spec = "conditional"}.
 #' @param seed Integer. RNG seed for the chain. Default \code{1L}.
@@ -85,9 +83,13 @@
 #'   \eqn{(K, \Gamma)} jointly from the un-normalised joint prior), or
 #'   \code{"hierarchical"} (sample \eqn{(K, \Gamma)} from the per-graph
 #'   normalized specification via the Z-ratio approximation).
-#' @param edge_inclusion_prob Probability in \eqn{(0, 1)} for the
-#'   Bernoulli edge prior used when \code{spec = "joint"}. Default
-#'   \code{0.5}. Ignored when \code{spec = "conditional"}.
+#' @param edge_inclusion_prob Probability in \eqn{(0, 1)} for the Bernoulli
+#'   edge prior the chain falls back to when \code{edge_prior} is
+#'   \code{NULL}. Default \code{0.5}. Read only when \code{spec = "joint"}
+#'   or \code{"hierarchical"} \emph{and} \code{edge_prior = NULL}: it is
+#'   ignored for \code{spec = "conditional"} (the graph is fixed) and
+#'   whenever an \code{edge_prior} object is supplied, including
+#'   \code{bernoulli_prior()}, which carries its own probability.
 #' @param update_method One of \code{"adaptive-metropolis"} (default) or
 #'   \code{"gibbs"}. Sampler driving the \code{spec = "joint"} chain; the
 #'   Gibbs chain uses the conjugate row and edge updates and needs no
@@ -97,13 +99,18 @@
 #'   or \code{\link{sbm_prior}()}, or \code{NULL} (default) for a Bernoulli
 #'   prior with probability \code{edge_inclusion_prob}. Only for
 #'   \code{spec = "joint"}.
-#' @param apply_correction Logical. For the hierarchical edge priors
-#'   (\code{beta_bernoulli_prior()}, \code{sbm_prior()}), apply the
-#'   normalizing-constant correction to the hyperparameter updates (default
-#'   \code{TRUE}; the correction table is built from the tilted prior
-#'   sampler and cached across calls). With \code{FALSE} the plain conjugate
-#'   updates are used, whose hyperparameter marginals do not match the
-#'   hyperpriors under the determinant tilt.
+#' @param apply_correction Logical. Apply the normalizing-constant correction
+#'   to the hyperparameter updates (default \code{TRUE}; the correction table
+#'   is built from the tilted prior sampler and cached across calls). With
+#'   \code{FALSE} the plain conjugate updates are used, whose hyperparameter
+#'   marginals do not match the hyperpriors under the determinant tilt. Read
+#'   only when \code{spec = "joint"} \emph{and} the edge prior is a
+#'   hierarchical one (\code{beta_bernoulli_prior()}, \code{sbm_prior()});
+#'   it is ignored for \code{spec = "conditional"} (the graph is fixed), for
+#'   \code{spec = "hierarchical"} (the per-edge Z-ratio carries the
+#'   normalizer instead, and the hyperparameter updates stay clean
+#'   conjugate), and for a Bernoulli edge prior (no hyperparameters to
+#'   correct).
 #' @param zratio_diagnostics Logical (default \code{TRUE}). Only for
 #'   \code{spec = "hierarchical"}: run the trust gauge
 #'   (\code{\link{summarize_zratio_gauge}}) on the returned chain and attach
@@ -196,7 +203,7 @@ sample_ggm_prior = function(
   n_warmup = 2e3,
   interaction_prior = normal_prior(scale = 1),
   precision_scale_prior = exponential_prior(eta = 1),
-  step_size = 0.1,
+  step_size = lifecycle::deprecated(),
   max_depth = 10L,
   seed = 1L,
   verbose = TRUE,
@@ -211,6 +218,15 @@ sample_ggm_prior = function(
 ) {
   spec = match.arg(spec)
   update_method = match.arg(update_method)
+  if(lifecycle::is_present(step_size)) {
+    lifecycle::deprecate_warn(
+      "0.2.1", "sample_ggm_prior(step_size = )",
+      details = paste0(
+        "The NUTS warmup determines its own initial step size; the supplied ",
+        "value is ignored."
+      )
+    )
+  }
   ep = if(is.null(edge_prior)) {
     NULL
   } else {
@@ -235,7 +251,6 @@ sample_ggm_prior = function(
   validate_integer(n_samples, "n_samples", min_value = 1L)
   validate_integer(n_warmup, "n_warmup", min_value = 0L)
   validate_integer(max_depth, "max_depth", min_value = 1L)
-  validate_finite_scalar(step_size, "step_size", positive = TRUE)
   validate_integer(seed, "seed", min_value = 0L)
   if(is.null(delta)) {
     delta = 0.5 * log(p)
@@ -277,7 +292,6 @@ sample_ggm_prior = function(
       scale_prior_type         = sp$scale_prior_type,
       gamma_shape              = sp$scale_shape,
       gamma_rate               = sp$scale_rate,
-      step_size                = step_size,
       max_depth                = as.integer(max_depth),
       seed                     = as.integer(seed),
       verbose                  = verbose,
@@ -367,6 +381,14 @@ 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_
+    ),
     seed = as.integer(seed),
     no_threads = 1L,
     progress_type = if(verbose) 2L else 0L,
@@ -545,16 +567,6 @@ validate_integer = function(x, name, min_value = 1L) {
   invisible(as.integer(x))
 }
 
-validate_finite_scalar = function(x, name, positive = FALSE) {
-  if(!is.numeric(x) || length(x) != 1L || is.na(x) || !is.finite(x)) {
-    stop(sprintf("'%s' must be a single finite numeric.", name))
-  }
-  if(positive && x <= 0) {
-    stop(sprintf("'%s' must be positive.", name))
-  }
-  invisible(x)
-}
-
 validate_ggm_prior_edge_indicators = function(edge_indicators, p) {
   if(is.null(edge_indicators)) {
     return(NULL)

R/sample_graph_prior.R

@@ -142,6 +142,10 @@ graph_prior_conditioning = function(ep, p, theta, allocations, block_probs) {
 #' a Bernoulli or Beta-Bernoulli prior, and \code{allocations} plus
 #' \code{block_probs} fix the block structure of a Stochastic-Block prior
 #' (reducing it to independent pair flips at the given block probabilities).
+#' Either way the graph law that is sampled is a Bernoulli prior at the fixed
+#' pair probabilities, under both specifications, so the returned
+#' \code{edge_prior} reports \code{"Bernoulli"} and the requested family is
+#' returned separately in \code{requested_edge_prior}.
 #'
 #' Under \code{spec = "joint"} the tilted graph law depends on the precision
 #' prior through its normalizer, so \code{interaction_prior},
@@ -192,8 +196,16 @@ graph_prior_conditioning = function(ep, p, theta, allocations, block_probs) {
 #'     \item{\code{allocations}}{Only with an unconditioned
 #'       \code{sbm_prior()}: integer matrix (\code{n_samples x p}) of sampled
 #'       block allocations.}
-#'     \item{\code{spec}, \code{edge_prior}, \code{p}}{The specification,
-#'       edge-prior family, and node count of the draw.}
+#'     \item{\code{edge_prior}}{The edge-prior family of the draw:
+#'       \code{"Bernoulli"} whenever \code{theta} or (\code{allocations},
+#'       \code{block_probs}) is supplied, since conditioning fixes the
+#'       hyperparameters and leaves independent pair flips; otherwise the
+#'       family of the \code{edge_prior} argument.}
+#'     \item{\code{requested_edge_prior}}{The family of the \code{edge_prior}
+#'       argument. Equal to \code{edge_prior} unless a conditioning argument
+#'       was supplied.}
+#'     \item{\code{spec}, \code{p}}{The specification and node count of the
+#'       draw.}
 #'   }
 #'
 #' @examples
@@ -256,6 +268,12 @@ sample_graph_prior = function(
   ep = unpack_indicator_prior(edge_prior, num_variables = p)
   cond_prob = graph_prior_conditioning(ep, p, theta, allocations, block_probs)
 
+  # Conditioning fixes the hyperparameters, which collapses either spec to
+  # independent pair flips at `cond_prob`: a Bernoulli prior, whichever family
+  # was requested. `edge_prior` reports the family that ran; the request is
+  # kept in `requested_edge_prior`.
+  run_edge_prior = if(is.null(cond_prob)) ep$edge_prior else "Bernoulli"
+
   pairs = graph_pair_indices(p)
   pair_names = paste0(pairs[, 1L], "-", pairs[, 2L])
 
@@ -282,7 +300,8 @@ sample_graph_prior = function(
       theta = draws$theta,
       allocations = draws$allocations,
       spec = spec,
-      edge_prior = ep$edge_prior,
+      edge_prior = run_edge_prior,
+      requested_edge_prior = ep$edge_prior,
       p = p
     )
     return(out)
@@ -330,7 +349,8 @@ sample_graph_prior = function(
     theta = res$theta,
     allocations = res$allocations,
     spec = spec,
-    edge_prior = ep$edge_prior,
+    edge_prior = run_edge_prior,
+    requested_edge_prior = ep$edge_prior,
     p = p
   )
 }

src/RcppExports.cpp

@@ -219,8 +219,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // sample_ggm_prior
-Rcpp::List sample_ggm_prior(int p, int n_samples, int n_warmup, double pairwise_scale, const std::string& interaction_prior_type, const std::string& scale_prior_type, double gamma_shape, double gamma_rate, double step_size, int max_depth, int seed, bool verbose, Rcpp::Nullable<Rcpp::IntegerMatrix> edge_indicators_nullable, double delta);
-RcppExport SEXP _bgms_sample_ggm_prior(SEXP pSEXP, SEXP n_samplesSEXP, SEXP n_warmupSEXP, SEXP pairwise_scaleSEXP, SEXP interaction_prior_typeSEXP, SEXP scale_prior_typeSEXP, SEXP gamma_shapeSEXP, SEXP gamma_rateSEXP, SEXP step_sizeSEXP, SEXP max_depthSEXP, SEXP seedSEXP, SEXP verboseSEXP, SEXP edge_indicators_nullableSEXP, SEXP deltaSEXP) {
+Rcpp::List sample_ggm_prior(int p, int n_samples, int n_warmup, double pairwise_scale, const std::string& interaction_prior_type, const std::string& scale_prior_type, double gamma_shape, double gamma_rate, int max_depth, int seed, bool verbose, Rcpp::Nullable<Rcpp::IntegerMatrix> edge_indicators_nullable, double delta);
+RcppExport SEXP _bgms_sample_ggm_prior(SEXP pSEXP, SEXP n_samplesSEXP, SEXP n_warmupSEXP, SEXP pairwise_scaleSEXP, SEXP interaction_prior_typeSEXP, SEXP scale_prior_typeSEXP, SEXP gamma_shapeSEXP, SEXP gamma_rateSEXP, SEXP max_depthSEXP, SEXP seedSEXP, SEXP verboseSEXP, SEXP edge_indicators_nullableSEXP, SEXP deltaSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -232,13 +232,12 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< const std::string& >::type scale_prior_type(scale_prior_typeSEXP);
     Rcpp::traits::input_parameter< double >::type gamma_shape(gamma_shapeSEXP);
     Rcpp::traits::input_parameter< double >::type gamma_rate(gamma_rateSEXP);
-    Rcpp::traits::input_parameter< double >::type step_size(step_sizeSEXP);
     Rcpp::traits::input_parameter< int >::type max_depth(max_depthSEXP);
     Rcpp::traits::input_parameter< int >::type seed(seedSEXP);
     Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP);
     Rcpp::traits::input_parameter< Rcpp::Nullable<Rcpp::IntegerMatrix> >::type edge_indicators_nullable(edge_indicators_nullableSEXP);
     Rcpp::traits::input_parameter< double >::type delta(deltaSEXP);
-    rcpp_result_gen = Rcpp::wrap(sample_ggm_prior(p, n_samples, n_warmup, pairwise_scale, interaction_prior_type, scale_prior_type, gamma_shape, gamma_rate, step_size, max_depth, seed, verbose, edge_indicators_nullable, delta));
+    rcpp_result_gen = Rcpp::wrap(sample_ggm_prior(p, n_samples, n_warmup, pairwise_scale, interaction_prior_type, scale_prior_type, gamma_shape, gamma_rate, max_depth, seed, verbose, edge_indicators_nullable, delta));
     return rcpp_result_gen;
 END_RCPP
 }
@@ -999,7 +998,7 @@ static const R_CallMethodDef CallEntries[] = {
     {"_bgms_ggm_test_gibbs_sweep", (DL_FUNC) &_bgms_ggm_test_gibbs_sweep, 8},
     {"_bgms_ggm_test_logp_and_gradient", (DL_FUNC) &_bgms_ggm_test_logp_and_gradient, 5},
     {"_bgms_ggm_test_forward_map", (DL_FUNC) &_bgms_ggm_test_forward_map, 2},
-    {"_bgms_sample_ggm_prior", (DL_FUNC) &_bgms_sample_ggm_prior, 14},
+    {"_bgms_sample_ggm_prior", (DL_FUNC) &_bgms_sample_ggm_prior, 13},
     {"_bgms_compute_ess_cpp", (DL_FUNC) &_bgms_compute_ess_cpp, 1},
     {"_bgms_compute_rhat_cpp", (DL_FUNC) &_bgms_compute_rhat_cpp, 1},
     {"_bgms_compute_indicator_ess_cpp", (DL_FUNC) &_bgms_compute_indicator_ess_cpp, 1},