From 6be09842e36c6ddfd6dd535978f910fc99490c34 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Fri, 15 May 2026 06:44:54 +0200 Subject: [PATCH 01/10] feat(ggm): determinant-tilt prior for the NUTS gradient Adds a non-negative tilt parameter delta to the GGM gradient engine: multiplies the prior by |K|^delta, which in Cholesky coordinates is + 2*delta * sum(psi). One line on the log-posterior in each gradient path (logp_and_gradient, logp_and_gradient_full) and +2*delta on each psi-bar entry in the adjoint. Cost O(p) per leapfrog step inside an O(p^3) gradient; effectively free. delta is plumbed through GGMGradientEngine::rebuild and stored on the engine. GGMModel exposes a set_determinant_tilt(delta) setter that forwards on the next ensure_constraint_structure() rebuild; the value is also carried in the copy constructor so per-chain clones target the intended posterior. Currently consumed only by the NUTS path. The MH ratios in ggm_model.cpp do not yet include the corresponding +delta * (log|K_prop| - log|K_curr|) term, so this is safe only when no MH moves fire (i.e. sample_ggm_prior with edge_selection = false). The MH-path tilt will land before delta is exposed on bgm(). --- src/ggm_gradient_interface.cpp | 4 +++- src/models/ggm/ggm_gradient.cpp | 20 +++++++++++++++++--- src/models/ggm/ggm_gradient.h | 6 +++++- src/models/ggm/ggm_model.cpp | 2 +- src/models/ggm/ggm_model.h | 19 ++++++++++++++++++- 5 files changed, 44 insertions(+), 7 deletions(-) diff --git a/src/ggm_gradient_interface.cpp b/src/ggm_gradient_interface.cpp index cff375da..bd58fc53 100644 --- a/src/ggm_gradient_interface.cpp +++ b/src/ggm_gradient_interface.cpp @@ -385,7 +385,8 @@ Rcpp::List sample_precision_prior( int max_depth = 10, int seed = 1, bool verbose = true, - Rcpp::Nullable edge_indicators_nullable = R_NilValue) + Rcpp::Nullable edge_indicators_nullable = R_NilValue, + double delta = 0.0) { // Build edge indicators (default: full graph, no constraints) arma::imat edge_indicators; @@ -407,6 +408,7 @@ Rcpp::List sample_precision_prior( arma::mat inc_prob(p, p, arma::fill::value(0.5)); GGMModel model(0, suf_stat, inc_prob, edge_indicators, false, std::move(ip), std::move(dp)); + model.set_determinant_tilt(delta); // Configure NUTS with the standard windowed warmup (dual averaging + // diagonal Welford mass-matrix adaptation). This is the same path used diff --git a/src/models/ggm/ggm_gradient.cpp b/src/models/ggm/ggm_gradient.cpp index 09dec889..229cad20 100644 --- a/src/models/ggm/ggm_gradient.cpp +++ b/src/models/ggm/ggm_gradient.cpp @@ -13,7 +13,8 @@ void GGMGradientEngine::rebuild( size_t n, const arma::mat& suf_stat, const BaseParameterPrior& interaction_prior, - const BaseParameterPrior& diagonal_prior) + const BaseParameterPrior& diagonal_prior, + double determinant_tilt) { structure_ = &structure; n_ = n; @@ -21,6 +22,7 @@ void GGMGradientEngine::rebuild( suf_stat_ = &suf_stat; interaction_prior_ = &interaction_prior; diagonal_prior_ = &diagonal_prior; + delta_ = determinant_tilt; } // ===================================================================== @@ -285,7 +287,11 @@ std::pair GGMGradientEngine::logp_and_gradient( log_diag_prior += diagonal_prior_->logp(0.5 * K(i, i)); } - double lp = log_lik + log_slab + log_diag_prior + fm.log_det_jacobian; + // Determinant tilt: adds delta_ * log|K| = 2*delta_ * sum(psi) to the + // log-prior. Pushes the chain away from the PD-cone boundary. + double log_tilt = delta_ * log_det_K; + + double lp = log_lik + log_slab + log_diag_prior + fm.log_det_jacobian + log_tilt; if (!std::isfinite(lp)) { return {lp, arma::vec(theta.n_elem, arma::fill::zeros)}; @@ -335,6 +341,8 @@ std::pair GGMGradientEngine::logp_and_gradient( if (q + 1 < p_) { psi_bar += static_cast(p_ - 1 - q); } + // Determinant tilt: d/dpsi_q [delta_ * 2 * sum(psi)] = 2 * delta_ + psi_bar += 2.0 * delta_; gradient(offset + d_q) = psi_bar; if (q == 0) continue; @@ -591,7 +599,11 @@ std::pair GGMGradientEngine::logp_and_gradient_full( pfaffian += arma::accu(arma::log(arma::diagvec(L_q))); } - double lp = log_lik + log_slab + log_diag_prior + ldj - pfaffian; + // Determinant tilt: adds delta_ * log|K| = 2*delta_ * sum(psi) to the + // log-prior. Pushes the chain away from the PD-cone boundary. + double log_tilt = delta_ * log_det_K; + + double lp = log_lik + log_slab + log_diag_prior + ldj - pfaffian + log_tilt; if (!std::isfinite(lp)) { return {lp, arma::vec(x.n_elem, arma::fill::zeros)}; @@ -676,6 +688,8 @@ std::pair GGMGradientEngine::logp_and_gradient_full( if (q + 1 < p_) { psi_bar += static_cast(p_ - 1 - q); } + // Determinant tilt: d/dpsi_q [delta_ * 2 * sum(psi)] = 2 * delta_ + psi_bar += 2.0 * delta_; gradient(offset + q) = psi_bar; } diff --git a/src/models/ggm/ggm_gradient.h b/src/models/ggm/ggm_gradient.h index d95d6b65..2b944114 100644 --- a/src/models/ggm/ggm_gradient.h +++ b/src/models/ggm/ggm_gradient.h @@ -79,7 +79,8 @@ class GGMGradientEngine { size_t n, const arma::mat& suf_stat, const BaseParameterPrior& interaction_prior, - const BaseParameterPrior& diagonal_prior); + const BaseParameterPrior& diagonal_prior, + double determinant_tilt = 0.0); /** * Forward map: theta -> (Phi, K, log|det J|). @@ -159,4 +160,7 @@ class GGMGradientEngine { const arma::mat* suf_stat_ = nullptr; const BaseParameterPrior* interaction_prior_ = nullptr; const BaseParameterPrior* diagonal_prior_ = nullptr; + // Determinant-tilt exponent: adds delta_ * log|K| to the (unnormalised) + // log-prior. delta_ = 0 recovers the untilted target. + double delta_ = 0.0; }; diff --git a/src/models/ggm/ggm_model.cpp b/src/models/ggm/ggm_model.cpp index 75f124e1..8b9496f7 100644 --- a/src/models/ggm/ggm_model.cpp +++ b/src/models/ggm/ggm_model.cpp @@ -12,7 +12,7 @@ void GGMModel::ensure_constraint_structure() { if (!constraint_dirty_) return; constraint_structure_.build(edge_indicators_); - gradient_engine_.rebuild(constraint_structure_, n_, suf_stat_, *interaction_prior_, *diagonal_prior_); + gradient_engine_.rebuild(constraint_structure_, n_, suf_stat_, *interaction_prior_, *diagonal_prior_, determinant_tilt_); constraint_dirty_ = false; theta_valid_ = false; } diff --git a/src/models/ggm/ggm_model.h b/src/models/ggm/ggm_model.h index 438f7704..f0cc5c5d 100644 --- a/src/models/ggm/ggm_model.h +++ b/src/models/ggm/ggm_model.h @@ -152,7 +152,8 @@ class GGMModel : public BaseModel { constraint_dirty_(other.constraint_dirty_), theta_valid_(other.theta_valid_), theta_(other.theta_), - pcg_lambda_cache_(other.pcg_lambda_cache_) + pcg_lambda_cache_(other.pcg_lambda_cache_), + determinant_tilt_(other.determinant_tilt_) {} /** @return true (GGM supports NUTS via free-element Cholesky gradient). */ @@ -207,6 +208,18 @@ class GGMModel : public BaseModel { */ void init_metropolis_adaptation(const WarmupSchedule& schedule) override; + /** + * Set the determinant-tilt exponent delta. Adds delta * log|K| to the + * log-prior, pushing the chain away from the PD-cone boundary. delta = 0 + * (default) recovers the untilted target. Currently consumed only by + * the NUTS gradient engine; the MH path is unchanged. Triggers an engine + * rebuild on the next gradient call. + */ + void set_determinant_tilt(double delta) { + determinant_tilt_ = delta; + constraint_dirty_ = true; + } + /** Shuffle edge visit order (random scan). */ void prepare_iteration() override; @@ -473,6 +486,10 @@ class GGMModel : public BaseModel { /// stays null and the stage-3b path in tune_proposal_sd is used instead). std::unique_ptr metropolis_adapter_; + // Determinant-tilt exponent (see set_determinant_tilt). Forwarded to + // GGMGradientEngine on every rebuild. + double determinant_tilt_ = 0.0; + /** Extract upper triangle of the precision matrix into a vector. */ arma::vec extract_upper_triangle() const { arma::vec result(dim_); From 987cc8bc4dd9a0fcb5c755549b46c51b163e9536 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Fri, 15 May 2026 06:45:47 +0200 Subject: [PATCH 02/10] refactor: rename sample_precision_prior to sample_ggm_prior The function samples from a prior whose density is specified on the partial-association scale K_yy = -K/2, not on the raw precision K, so "precision" was misleading. Output samples are still entries of K. Renames the R wrapper, the Rcpp export (sample_ggm_prior_cpp), the internal C++ function, the test file, and the docstring; regenerates RcppExports and roxygen docs to match. --- NAMESPACE | 2 +- R/RcppExports.R | 4 +- ...e_precision_prior.R => sample_ggm_prior.R} | 70 ++++++++++++------- ...precision_prior.Rd => sample_ggm_prior.Rd} | 66 ++++++++++------- src/RcppExports.cpp | 11 +-- src/ggm_gradient_interface.cpp | 11 +-- ...cision-prior.R => test-sample-ggm-prior.R} | 9 +-- 7 files changed, 102 insertions(+), 71 deletions(-) rename R/{sample_precision_prior.R => sample_ggm_prior.R} (71%) rename man/{sample_precision_prior.Rd => sample_ggm_prior.Rd} (58%) rename tests/testthat/{test-sample-precision-prior.R => test-sample-ggm-prior.R} (95%) diff --git a/NAMESPACE b/NAMESPACE index 2765cf8a..2404bd92 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -70,7 +70,7 @@ export(extract_sbm) export(gamma_prior) export(mrfSampler) export(normal_prior) -export(sample_precision_prior) +export(sample_ggm_prior) export(sbm_prior) export(simulate_mrf) import(RcppParallel) diff --git a/R/RcppExports.R b/R/RcppExports.R index 47938d7b..4080c26d 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -49,8 +49,8 @@ ggm_test_leapfrog_constrained_checked <- function(x0, r0, step_size, n_steps, su .Call(`_bgms_ggm_test_leapfrog_constrained_checked`, x0, r0, step_size, n_steps, suf_stat, n, edge_indicators, pairwise_scale, reverse_check_tol, inv_mass_in) } -sample_precision_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) { - .Call(`_bgms_sample_precision_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) +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) } .compute_ess_cpp <- function(array3d) { diff --git a/R/sample_precision_prior.R b/R/sample_ggm_prior.R similarity index 71% rename from R/sample_precision_prior.R rename to R/sample_ggm_prior.R index a68c23db..04a65b28 100644 --- a/R/sample_precision_prior.R +++ b/R/sample_ggm_prior.R @@ -1,30 +1,36 @@ -#' Sample Precision Matrices from the GGM Prior +#' Sample from the GGM (Partial-Association) Prior #' -#' Draws precision matrices \eqn{K} from the prior of a Gaussian graphical -#' model using the same constrained NUTS sampler that drives \code{\link{bgm}} -#' for continuous data. The likelihood is omitted (\eqn{n = 0}, -#' \eqn{S = 0}), so the chain targets the prior alone. +#' Draws from the prior of a Gaussian graphical model using the same +#' constrained NUTS sampler that drives \code{\link{bgm}} for continuous +#' data. The likelihood is omitted (\eqn{n = 0}, \eqn{S = 0}), so the chain +#' targets the prior alone. #' -#' Off-diagonals are placed on the association scale -#' \eqn{K_{yy,ij} = -K_{ij}/2} and assigned the supplied -#' \code{interaction_prior}. A \code{normal_prior(scale = s)} therefore -#' constrains \eqn{K_{yy}} with standard deviation \eqn{s}, equivalent to a -#' \eqn{\textrm{Normal}(0, 2s)} prior on \eqn{K_{ij}} itself. The -#' diagonals \eqn{K_{ii}} are drawn from the supplied -#' \code{precision_scale_prior}. When \code{edge_indicators} is supplied, -#' off-diagonals at excluded positions are constrained to zero throughout -#' the chain. +#' The priors are specified on the partial-association scale +#' \eqn{K_{yy} = -K/2}: \code{interaction_prior} acts on +#' \eqn{K_{yy,ij} = -K_{ij}/2}, and \code{precision_scale_prior} acts on +#' \eqn{-K_{yy,ii} = K_{ii}/2}. The same convention is used by +#' \code{\link{bgm}} and by the continuous block of the mixed-MRF model, so +#' a prior argument passed here means the same distribution it would mean +#' there. Output samples are reported as entries of \eqn{K}; convert with +#' \eqn{K_{yy} = -K/2} if you want them on the partial-association scale. +#' +#' When \code{edge_indicators} is supplied, off-diagonals at excluded +#' positions are constrained to zero throughout the chain. #' #' @param p Integer. Dimension of the precision matrix (\eqn{p \ge 2}). #' @param n_samples Integer. Number of post-warmup draws to keep. #' @param n_warmup Integer. NUTS warmup iterations. Default \code{2000}. #' @param interaction_prior A \code{bgms_parameter_prior} for the -#' off-diagonal entries. Use \code{\link{cauchy_prior}()} or -#' \code{\link{normal_prior}()}; \code{\link{beta_prime_prior}()} is not -#' supported here. Default: \code{cauchy_prior(scale = 2.5)}. -#' @param precision_scale_prior A \code{bgms_scale_prior} for the diagonal -#' entries of \eqn{K}. Use \code{\link{gamma_prior}()} or -#' \code{\link{exponential_prior}()}. Default: \code{gamma_prior(1, 1)}. +#' partial-association off-diagonals \eqn{K_{yy,ij} = -K_{ij}/2}. Use +#' \code{\link{cauchy_prior}()} or \code{\link{normal_prior}()}; +#' \code{\link{beta_prime_prior}()} is not supported here. Default: +#' \code{cauchy_prior(scale = 2.5)} (i.e. \eqn{K_{ij}} has an implied +#' \eqn{\textrm{Cauchy}(0, 5)} prior). +#' @param precision_scale_prior A \code{bgms_scale_prior} for +#' \eqn{K_{ii}/2}. Use \code{\link{gamma_prior}()} or +#' \code{\link{exponential_prior}()}. Default: \code{gamma_prior(1, 1)}, +#' which implies \eqn{K_{ii}/2 \sim \textrm{Exp}(1)} and therefore +#' \eqn{K_{ii} \sim \textrm{Exp}(1/2)} (mean \eqn{2}). #' @param step_size Positive numeric. Initial NUTS step size used to seed #' dual-averaging adaptation. Default \code{0.1}. #' @param max_depth Integer. Maximum NUTS tree depth. Default \code{10}. @@ -33,6 +39,10 @@ #' @param edge_indicators Optional integer \eqn{p \times p} matrix with #' \code{1} = edge included, \code{0} = excluded. Must be symmetric with #' \code{1}s on the diagonal. Default: full graph (all edges included). +#' @param delta Non-negative numeric. Determinant-tilt exponent: multiplies +#' the prior by \eqn{|K|^{\delta}}, pushing the chain away from the +#' positive-definite cone boundary. \code{delta = 0} (default) recovers +#' the untilted prior. #' #' @return A list with components #' \describe{ @@ -61,7 +71,7 @@ #' @examples #' \donttest{ #' # Default Cauchy(0, 2.5) off-diagonal, Gamma(1, 1) diagonal, p = 4. -#' draws = sample_precision_prior( +#' draws = sample_ggm_prior( #' p = 4, n_samples = 200, n_warmup = 200, #' verbose = FALSE #' ) @@ -72,7 +82,7 @@ #' # Sparser graph: drop the (1, 4) edge. #' E = matrix(1L, 4, 4) #' E[1, 4] = E[4, 1] = 0L -#' draws = sample_precision_prior( +#' draws = sample_ggm_prior( #' p = 4, n_samples = 200, n_warmup = 200, #' edge_indicators = E, verbose = FALSE #' ) @@ -80,7 +90,7 @@ #' all(draws$K_offdiag[, "K_1_4"] == 0) # TRUE #' } #' @export -sample_precision_prior = function( +sample_ggm_prior = function( p, n_samples, n_warmup = 2e3, @@ -90,7 +100,8 @@ sample_precision_prior = function( max_depth = 10L, seed = 1L, verbose = TRUE, - edge_indicators = NULL + edge_indicators = NULL, + delta = 0 ) { validate_positive_integer(p, "p", min_value = 2L) validate_positive_integer(n_samples, "n_samples", min_value = 1L) @@ -98,6 +109,10 @@ sample_precision_prior = function( validate_positive_integer(max_depth, "max_depth", min_value = 1L) validate_finite_scalar(step_size, "step_size", positive = TRUE) validate_positive_integer(seed, "seed", min_value = 0L) + if(!is.numeric(delta) || length(delta) != 1L || is.na(delta) || + !is.finite(delta) || delta < 0) { + stop("'delta' must be a single finite non-negative numeric.") + } if(!is.logical(verbose) || length(verbose) != 1L || is.na(verbose)) { stop("'verbose' must be TRUE or FALSE.") } @@ -106,14 +121,14 @@ sample_precision_prior = function( if(identical(ip$interaction_prior_type, "beta-prime")) { stop( "beta_prime_prior() is not supported for 'interaction_prior' in ", - "sample_precision_prior(). Use cauchy_prior() or normal_prior()." + "sample_ggm_prior(). Use cauchy_prior() or normal_prior()." ) } sp = unpack_scale_prior(precision_scale_prior) edge_indicators = validate_ggm_prior_edge_indicators(edge_indicators, p) - sample_precision_prior_cpp( + sample_ggm_prior_cpp( p = as.integer(p), n_samples = as.integer(n_samples), n_warmup = as.integer(n_warmup), @@ -126,7 +141,8 @@ sample_precision_prior = function( max_depth = as.integer(max_depth), seed = as.integer(seed), verbose = verbose, - edge_indicators_nullable = edge_indicators + edge_indicators_nullable = edge_indicators, + delta = as.numeric(delta) ) } diff --git a/man/sample_precision_prior.Rd b/man/sample_ggm_prior.Rd similarity index 58% rename from man/sample_precision_prior.Rd rename to man/sample_ggm_prior.Rd index 4099c241..05f53f78 100644 --- a/man/sample_precision_prior.Rd +++ b/man/sample_ggm_prior.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/sample_precision_prior.R -\name{sample_precision_prior} -\alias{sample_precision_prior} -\title{Sample Precision Matrices from the GGM Prior} +% Please edit documentation in R/sample_ggm_prior.R +\name{sample_ggm_prior} +\alias{sample_ggm_prior} +\title{Sample from the GGM (Partial-Association) Prior} \usage{ -sample_precision_prior( +sample_ggm_prior( p, n_samples, n_warmup = 2000, @@ -14,7 +14,8 @@ sample_precision_prior( max_depth = 10L, seed = 1L, verbose = TRUE, - edge_indicators = NULL + edge_indicators = NULL, + delta = 0 ) } \arguments{ @@ -25,13 +26,17 @@ sample_precision_prior( \item{n_warmup}{Integer. NUTS warmup iterations. Default \code{2000}.} \item{interaction_prior}{A \code{bgms_parameter_prior} for the -off-diagonal entries. Use \code{\link{cauchy_prior}()} or -\code{\link{normal_prior}()}; \code{\link{beta_prime_prior}()} is not -supported here. Default: \code{cauchy_prior(scale = 2.5)}.} +partial-association off-diagonals \eqn{K_{yy,ij} = -K_{ij}/2}. Use +\code{\link{cauchy_prior}()} or \code{\link{normal_prior}()}; +\code{\link{beta_prime_prior}()} is not supported here. Default: +\code{cauchy_prior(scale = 2.5)} (i.e. \eqn{K_{ij}} has an implied +\eqn{\textrm{Cauchy}(0, 5)} prior).} -\item{precision_scale_prior}{A \code{bgms_scale_prior} for the diagonal -entries of \eqn{K}. Use \code{\link{gamma_prior}()} or -\code{\link{exponential_prior}()}. Default: \code{gamma_prior(1, 1)}.} +\item{precision_scale_prior}{A \code{bgms_scale_prior} for +\eqn{K_{ii}/2}. Use \code{\link{gamma_prior}()} or +\code{\link{exponential_prior}()}. Default: \code{gamma_prior(1, 1)}, +which implies \eqn{K_{ii}/2 \sim \textrm{Exp}(1)} and therefore +\eqn{K_{ii} \sim \textrm{Exp}(1/2)} (mean \eqn{2}).} \item{step_size}{Positive numeric. Initial NUTS step size used to seed dual-averaging adaptation. Default \code{0.1}.} @@ -45,6 +50,11 @@ dual-averaging adaptation. Default \code{0.1}.} \item{edge_indicators}{Optional integer \eqn{p \times p} matrix with \code{1} = edge included, \code{0} = excluded. Must be symmetric with \code{1}s on the diagonal. Default: full graph (all edges included).} + +\item{delta}{Non-negative numeric. Determinant-tilt exponent: multiplies +the prior by \eqn{|K|^{\delta}}, pushing the chain away from the +positive-definite cone boundary. \code{delta = 0} (default) recovers +the untilted prior.} } \value{ A list with components @@ -68,26 +78,28 @@ the columns of \code{K_diag}.} } } \description{ -Draws precision matrices \eqn{K} from the prior of a Gaussian graphical -model using the same constrained NUTS sampler that drives \code{\link{bgm}} -for continuous data. The likelihood is omitted (\eqn{n = 0}, -\eqn{S = 0}), so the chain targets the prior alone. +Draws from the prior of a Gaussian graphical model using the same +constrained NUTS sampler that drives \code{\link{bgm}} for continuous +data. The likelihood is omitted (\eqn{n = 0}, \eqn{S = 0}), so the chain +targets the prior alone. } \details{ -Off-diagonals are placed on the association scale -\eqn{K_{yy,ij} = -K_{ij}/2} and assigned the supplied -\code{interaction_prior}. A \code{normal_prior(scale = s)} therefore -constrains \eqn{K_{yy}} with standard deviation \eqn{s}, equivalent to a -\eqn{\textrm{Normal}(0, 2s)} prior on \eqn{K_{ij}} itself. The -diagonals \eqn{K_{ii}} are drawn from the supplied -\code{precision_scale_prior}. When \code{edge_indicators} is supplied, -off-diagonals at excluded positions are constrained to zero throughout -the chain. +The priors are specified on the partial-association scale +\eqn{K_{yy} = -K/2}: \code{interaction_prior} acts on +\eqn{K_{yy,ij} = -K_{ij}/2}, and \code{precision_scale_prior} acts on +\eqn{-K_{yy,ii} = K_{ii}/2}. The same convention is used by +\code{\link{bgm}} and by the continuous block of the mixed-MRF model, so +a prior argument passed here means the same distribution it would mean +there. Output samples are reported as entries of \eqn{K}; convert with +\eqn{K_{yy} = -K/2} if you want them on the partial-association scale. + +When \code{edge_indicators} is supplied, off-diagonals at excluded +positions are constrained to zero throughout the chain. } \examples{ \donttest{ # Default Cauchy(0, 2.5) off-diagonal, Gamma(1, 1) diagonal, p = 4. -draws = sample_precision_prior( +draws = sample_ggm_prior( p = 4, n_samples = 200, n_warmup = 200, verbose = FALSE ) @@ -98,7 +110,7 @@ head(draws$K_offdiag) # Sparser graph: drop the (1, 4) edge. E = matrix(1L, 4, 4) E[1, 4] = E[4, 1] = 0L -draws = sample_precision_prior( +draws = sample_ggm_prior( p = 4, n_samples = 200, n_warmup = 200, edge_indicators = E, verbose = FALSE ) diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index b07142dc..fe27010b 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -218,9 +218,9 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// sample_precision_prior -Rcpp::List sample_precision_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 edge_indicators_nullable); -RcppExport SEXP _bgms_sample_precision_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) { +// 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 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) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -237,7 +237,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< int >::type seed(seedSEXP); Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable >::type edge_indicators_nullable(edge_indicators_nullableSEXP); - rcpp_result_gen = Rcpp::wrap(sample_precision_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)); + 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)); return rcpp_result_gen; END_RCPP } @@ -798,7 +799,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bgms_ggm_test_project_momentum", (DL_FUNC) &_bgms_ggm_test_project_momentum, 4}, {"_bgms_ggm_test_leapfrog_constrained", (DL_FUNC) &_bgms_ggm_test_leapfrog_constrained, 9}, {"_bgms_ggm_test_leapfrog_constrained_checked", (DL_FUNC) &_bgms_ggm_test_leapfrog_constrained_checked, 10}, - {"_bgms_sample_precision_prior", (DL_FUNC) &_bgms_sample_precision_prior, 13}, + {"_bgms_sample_ggm_prior", (DL_FUNC) &_bgms_sample_ggm_prior, 14}, {"_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}, diff --git a/src/ggm_gradient_interface.cpp b/src/ggm_gradient_interface.cpp index bd58fc53..c92f877a 100644 --- a/src/ggm_gradient_interface.cpp +++ b/src/ggm_gradient_interface.cpp @@ -2,7 +2,7 @@ // // Exposes logp_and_gradient, forward_map, project_position, // project_momentum, and constrained leapfrog to R for validation. -// Also exposes sample_precision_prior() for sampling precision matrices +// Also exposes sample_ggm_prior() for sampling from the GGM prior // from the GGM prior using NUTS. #include @@ -344,7 +344,8 @@ Rcpp::List ggm_test_leapfrog_constrained_checked( // ----------------------------------------------------------------------------- -// sample_precision_prior: Sample precision matrices from the GGM prior using NUTS +// sample_ggm_prior: Sample from the GGM prior (K_yy = -K/2 partial-association +// scale; output reported as K) using NUTS. // ----------------------------------------------------------------------------- // Uses the Cholesky parameterization with NUTS. By setting n=0 and S=0, // the likelihood vanishes and the sampler targets the prior: @@ -371,8 +372,8 @@ Rcpp::List ggm_test_leapfrog_constrained_checked( // adaptation, Welford diagonal mass-matrix adaptation) as sample_ggm and bgm. // ----------------------------------------------------------------------------- -// [[Rcpp::export(name = "sample_precision_prior_cpp")]] -Rcpp::List sample_precision_prior( +// [[Rcpp::export(name = "sample_ggm_prior_cpp")]] +Rcpp::List sample_ggm_prior( int p, int n_samples, int n_warmup = 1000, @@ -442,7 +443,7 @@ Rcpp::List sample_precision_prior( if(results.empty() || results[0].error) { std::string msg = results.empty() ? "empty result" : results[0].error_msg; - Rcpp::stop("sample_precision_prior: chain failed (%s)", msg.c_str()); + Rcpp::stop("sample_ggm_prior: chain failed (%s)", msg.c_str()); } // Samples are stored as the upper triangle of K (p(p+1)/2 elements per diff --git a/tests/testthat/test-sample-precision-prior.R b/tests/testthat/test-sample-ggm-prior.R similarity index 95% rename from tests/testthat/test-sample-precision-prior.R rename to tests/testthat/test-sample-ggm-prior.R index d87f497e..47afa520 100644 --- a/tests/testthat/test-sample-precision-prior.R +++ b/tests/testthat/test-sample-ggm-prior.R @@ -1,9 +1,9 @@ -# Tests for sample_precision_prior(): structure, prior plumbing, edge constraints, +# Tests for sample_ggm_prior(): structure, prior plumbing, edge constraints, # and input validation. Sampling correctness (uniform SBC ranks) is covered # in test-sbc-ggm.R. short_run = function(p, n_samples = 50L, n_warmup = 50L, seed = 1L, ...) { - sample_precision_prior( + sample_ggm_prior( p = p, n_samples = n_samples, n_warmup = n_warmup, seed = seed, verbose = FALSE, ... ) @@ -12,7 +12,7 @@ short_run = function(p, n_samples = 50L, n_warmup = 50L, seed = 1L, ...) { # ---- Return value structure -------------------------------------------------- -test_that("sample_precision_prior returns the documented list shape", { +test_that("sample_ggm_prior returns the documented list shape", { p = 4L n_samples = 60L draws = short_run(p = p, n_samples = n_samples) @@ -139,7 +139,7 @@ test_that("invalid scalar arguments error early with informative messages", { "'max_depth' must be >= 1" ) expect_error( - sample_precision_prior(p = 3L, n_samples = 10L, n_warmup = 10L, verbose = NA), + sample_ggm_prior(p = 3L, n_samples = 10L, n_warmup = 10L, verbose = NA), "'verbose' must be TRUE or FALSE" ) }) @@ -208,3 +208,4 @@ test_that("malformed edge_indicators are rejected", { "0 or 1" ) }) + From 50e8625e8aa9821c99bdd089e1827e03b500f4e7 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Fri, 15 May 2026 08:49:03 +0200 Subject: [PATCH 03/10] feat(mixed-mrf): determinant-tilt prior for the Kyy NUTS gradient Parallel of D for GGM. Adds a non-negative tilt parameter delta to the mixed-MRF Kyy block: multiplies the prior by |Kyy|^delta, which in Cholesky coordinates is + 2*delta * sum(psi). Two added lines per gradient path (logp and psi_bar) in mixed_mrf_gradient.cpp; cost is O(q) per leapfrog step inside an O(q^3) gradient. MixedMRFModel gains a determinant_tilt_yy_ member (default 0.0) and a set_determinant_tilt_yy() setter. The value is carried in the copy constructor so per-chain clones see the intended posterior. Currently consumed only by the NUTS gradient paths; the continuous-block MH ratios in mixed_mrf_metropolis.cpp are not yet updated. mixed_test_logp_and_gradient_full grows a delta argument so the tilt is exercisable from R. Finite-difference tests added to test-mixed-gradient-pfaffian.R cover: - delta = 0 bit-identical to no-tilt call - delta > 0 on full Kyy - delta > 0 on sparse Kyy with random diagonal mass - analytic check: tilted logp - untilted logp == delta * log|Kyy| All pass at max relative error < 1e-5. RcppExports regenerated via Rcpp::compileAttributes(). No regressions in the full test suite (6708 pass, 0 fail). --- R/RcppExports.R | 4 +- src/RcppExports.cpp | 9 ++- src/mixed_gradient_interface.cpp | 4 +- src/models/mixed/mixed_mrf_gradient.cpp | 22 ++++- src/models/mixed/mixed_mrf_model.cpp | 1 + src/models/mixed/mixed_mrf_model.h | 16 ++++ tests/testthat/test-mixed-gradient-pfaffian.R | 81 +++++++++++++++++-- 7 files changed, 120 insertions(+), 17 deletions(-) diff --git a/R/RcppExports.R b/R/RcppExports.R index 4080c26d..3a946468 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -69,8 +69,8 @@ mixed_test_logp_and_gradient <- function(params, discrete_observations, continuo .Call(`_bgms_mixed_test_logp_and_gradient`, params, discrete_observations, continuous_observations, num_categories, is_ordinal_variable, baseline_category, edge_indicators, pairwise_scale, main_alpha, main_beta, interaction_prior_type, threshold_prior_type, threshold_scale, means_prior_type, means_scale, diagonal_prior_type, diagonal_shape, diagonal_rate) } -mixed_test_logp_and_gradient_full <- function(params, discrete_observations, continuous_observations, num_categories, is_ordinal_variable, baseline_category, edge_indicators, pairwise_scale, main_alpha = 1.0, main_beta = 1.0, interaction_prior_type = "cauchy", threshold_prior_type = "beta-prime", threshold_scale = 1.0, means_prior_type = "normal", means_scale = 1.0, diagonal_prior_type = "gamma", diagonal_shape = 1.0, diagonal_rate = 1.0, inv_mass_diag = NULL) { - .Call(`_bgms_mixed_test_logp_and_gradient_full`, params, discrete_observations, continuous_observations, num_categories, is_ordinal_variable, baseline_category, edge_indicators, pairwise_scale, main_alpha, main_beta, interaction_prior_type, threshold_prior_type, threshold_scale, means_prior_type, means_scale, diagonal_prior_type, diagonal_shape, diagonal_rate, inv_mass_diag) +mixed_test_logp_and_gradient_full <- function(params, discrete_observations, continuous_observations, num_categories, is_ordinal_variable, baseline_category, edge_indicators, pairwise_scale, main_alpha = 1.0, main_beta = 1.0, interaction_prior_type = "cauchy", threshold_prior_type = "beta-prime", threshold_scale = 1.0, means_prior_type = "normal", means_scale = 1.0, diagonal_prior_type = "gamma", diagonal_shape = 1.0, diagonal_rate = 1.0, inv_mass_diag = NULL, delta = 0.0) { + .Call(`_bgms_mixed_test_logp_and_gradient_full`, params, discrete_observations, continuous_observations, num_categories, is_ordinal_variable, baseline_category, edge_indicators, pairwise_scale, main_alpha, main_beta, interaction_prior_type, threshold_prior_type, threshold_scale, means_prior_type, means_scale, diagonal_prior_type, diagonal_shape, diagonal_rate, inv_mass_diag, delta) } mixed_test_project_position <- function(x, inv_mass, discrete_observations, continuous_observations, num_categories, is_ordinal_variable, baseline_category, edge_indicators, pairwise_scale, main_alpha = 1.0, main_beta = 1.0, interaction_prior_type = "cauchy", threshold_prior_type = "beta-prime", threshold_scale = 1.0) { diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index fe27010b..7bc2ba9f 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -304,8 +304,8 @@ BEGIN_RCPP END_RCPP } // mixed_test_logp_and_gradient_full -Rcpp::List mixed_test_logp_and_gradient_full(const arma::vec& params, const arma::imat& discrete_observations, const arma::mat& continuous_observations, const arma::ivec& num_categories, const arma::uvec& is_ordinal_variable, const arma::ivec& baseline_category, const arma::imat& edge_indicators, double pairwise_scale, double main_alpha, double main_beta, std::string interaction_prior_type, std::string threshold_prior_type, double threshold_scale, std::string means_prior_type, double means_scale, std::string diagonal_prior_type, double diagonal_shape, double diagonal_rate, Rcpp::Nullable inv_mass_diag); -RcppExport SEXP _bgms_mixed_test_logp_and_gradient_full(SEXP paramsSEXP, SEXP discrete_observationsSEXP, SEXP continuous_observationsSEXP, SEXP num_categoriesSEXP, SEXP is_ordinal_variableSEXP, SEXP baseline_categorySEXP, SEXP edge_indicatorsSEXP, SEXP pairwise_scaleSEXP, SEXP main_alphaSEXP, SEXP main_betaSEXP, SEXP interaction_prior_typeSEXP, SEXP threshold_prior_typeSEXP, SEXP threshold_scaleSEXP, SEXP means_prior_typeSEXP, SEXP means_scaleSEXP, SEXP diagonal_prior_typeSEXP, SEXP diagonal_shapeSEXP, SEXP diagonal_rateSEXP, SEXP inv_mass_diagSEXP) { +Rcpp::List mixed_test_logp_and_gradient_full(const arma::vec& params, const arma::imat& discrete_observations, const arma::mat& continuous_observations, const arma::ivec& num_categories, const arma::uvec& is_ordinal_variable, const arma::ivec& baseline_category, const arma::imat& edge_indicators, double pairwise_scale, double main_alpha, double main_beta, std::string interaction_prior_type, std::string threshold_prior_type, double threshold_scale, std::string means_prior_type, double means_scale, std::string diagonal_prior_type, double diagonal_shape, double diagonal_rate, Rcpp::Nullable inv_mass_diag, double delta); +RcppExport SEXP _bgms_mixed_test_logp_and_gradient_full(SEXP paramsSEXP, SEXP discrete_observationsSEXP, SEXP continuous_observationsSEXP, SEXP num_categoriesSEXP, SEXP is_ordinal_variableSEXP, SEXP baseline_categorySEXP, SEXP edge_indicatorsSEXP, SEXP pairwise_scaleSEXP, SEXP main_alphaSEXP, SEXP main_betaSEXP, SEXP interaction_prior_typeSEXP, SEXP threshold_prior_typeSEXP, SEXP threshold_scaleSEXP, SEXP means_prior_typeSEXP, SEXP means_scaleSEXP, SEXP diagonal_prior_typeSEXP, SEXP diagonal_shapeSEXP, SEXP diagonal_rateSEXP, SEXP inv_mass_diagSEXP, SEXP deltaSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -328,7 +328,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< double >::type diagonal_shape(diagonal_shapeSEXP); Rcpp::traits::input_parameter< double >::type diagonal_rate(diagonal_rateSEXP); Rcpp::traits::input_parameter< Rcpp::Nullable >::type inv_mass_diag(inv_mass_diagSEXP); - rcpp_result_gen = Rcpp::wrap(mixed_test_logp_and_gradient_full(params, discrete_observations, continuous_observations, num_categories, is_ordinal_variable, baseline_category, edge_indicators, pairwise_scale, main_alpha, main_beta, interaction_prior_type, threshold_prior_type, threshold_scale, means_prior_type, means_scale, diagonal_prior_type, diagonal_shape, diagonal_rate, inv_mass_diag)); + Rcpp::traits::input_parameter< double >::type delta(deltaSEXP); + rcpp_result_gen = Rcpp::wrap(mixed_test_logp_and_gradient_full(params, discrete_observations, continuous_observations, num_categories, is_ordinal_variable, baseline_category, edge_indicators, pairwise_scale, main_alpha, main_beta, interaction_prior_type, threshold_prior_type, threshold_scale, means_prior_type, means_scale, diagonal_prior_type, diagonal_shape, diagonal_rate, inv_mass_diag, delta)); return rcpp_result_gen; END_RCPP } @@ -804,7 +805,7 @@ static const R_CallMethodDef CallEntries[] = { {"_bgms_compute_rhat_cpp", (DL_FUNC) &_bgms_compute_rhat_cpp, 1}, {"_bgms_compute_indicator_ess_cpp", (DL_FUNC) &_bgms_compute_indicator_ess_cpp, 1}, {"_bgms_mixed_test_logp_and_gradient", (DL_FUNC) &_bgms_mixed_test_logp_and_gradient, 18}, - {"_bgms_mixed_test_logp_and_gradient_full", (DL_FUNC) &_bgms_mixed_test_logp_and_gradient_full, 19}, + {"_bgms_mixed_test_logp_and_gradient_full", (DL_FUNC) &_bgms_mixed_test_logp_and_gradient_full, 20}, {"_bgms_mixed_test_project_position", (DL_FUNC) &_bgms_mixed_test_project_position, 14}, {"_bgms_mixed_test_project_momentum", (DL_FUNC) &_bgms_mixed_test_project_momentum, 15}, {"_bgms_mixed_test_leapfrog_constrained", (DL_FUNC) &_bgms_mixed_test_leapfrog_constrained, 17}, diff --git a/src/mixed_gradient_interface.cpp b/src/mixed_gradient_interface.cpp index c113c3e0..00086c23 100644 --- a/src/mixed_gradient_interface.cpp +++ b/src/mixed_gradient_interface.cpp @@ -76,7 +76,8 @@ Rcpp::List mixed_test_logp_and_gradient_full( std::string diagonal_prior_type = "gamma", double diagonal_shape = 1.0, double diagonal_rate = 1.0, - Rcpp::Nullable inv_mass_diag = R_NilValue) + Rcpp::Nullable inv_mass_diag = R_NilValue, + double delta = 0.0) { size_t p = discrete_observations.n_cols; size_t q = continuous_observations.n_cols; @@ -102,6 +103,7 @@ Rcpp::List mixed_test_logp_and_gradient_full( if(inv_mass_diag.isNotNull()) { model.set_inv_mass(Rcpp::as(inv_mass_diag)); } + model.set_determinant_tilt_yy(delta); auto result = model.logp_and_gradient_full(params); diff --git a/src/models/mixed/mixed_mrf_gradient.cpp b/src/models/mixed/mixed_mrf_gradient.cpp index 2cac4a2c..1b6d657b 100644 --- a/src/models/mixed/mixed_mrf_gradient.cpp +++ b/src/models/mixed/mixed_mrf_gradient.cpp @@ -652,6 +652,11 @@ std::pair MixedMRFModel::logp_and_gradient( logp += static_cast(q_ + 1 - j) * std::log(temp_cholesky(j, j)); } + // Determinant tilt on the Kyy block: adds delta * log|Kyy| = 2*delta * sum(psi) + // to the log-prior. Pushes the continuous-block precision matrix away from + // the PD-cone boundary. delta = 0 recovers the untilted target. + logp += determinant_tilt_yy_ * temp_log_det; + const auto& cs = chol_constraint_structure_; arma::mat Aq_buf; std::vector G_chol(q_); @@ -708,14 +713,16 @@ std::pair MixedMRFModel::logp_and_gradient( } // Extract position gradient from R_bar with the unified weight (q+1-j) - // on the diagonal-psi entries. + // on the diagonal-psi entries. The determinant tilt adds +2*delta to the + // psi-bar entries (d/dpsi [delta * 2 * sum(psi)] = 2 * delta). size_t gidx = static_cast(chol_grad_offset_); for(size_t j = 0; j < q_; ++j) { double w_j = static_cast(q_ + 1 - j); for(size_t i = 0; i < j; ++i) { grad(gidx++) = R_bar(i, j); } - grad(gidx++) = R_bar(j, j) * temp_cholesky(j, j) + w_j; + grad(gidx++) = R_bar(j, j) * temp_cholesky(j, j) + w_j + + 2.0 * determinant_tilt_yy_; } return {logp, grad}; @@ -1177,6 +1184,11 @@ std::pair MixedMRFModel::logp_and_gradient_full( logp += static_cast(q_ + 1 - j) * std::log(temp_cholesky(j, j)); } + // Determinant tilt on the Kyy block: adds delta * log|Kyy| = 2*delta * sum(psi) + // to the log-prior. Pushes the continuous-block precision matrix away from + // the PD-cone boundary. delta = 0 recovers the untilted target. + logp += determinant_tilt_yy_ * temp_log_det; + const auto& cs = chol_constraint_structure_; const bool identity_mass = inv_mass_.is_empty(); arma::mat Aq_buf; @@ -1251,14 +1263,16 @@ std::pair MixedMRFModel::logp_and_gradient_full( } // Extract position gradient from R_bar with the unified weight (q+1-j) - // on the diagonal-psi entries. + // on the diagonal-psi entries. The determinant tilt adds +2*delta to the + // psi-bar entries (d/dpsi [delta * 2 * sum(psi)] = 2 * delta). size_t gidx = chol_offset; for(size_t j = 0; j < q_; ++j) { double w_j = static_cast(q_ + 1 - j); for(size_t i = 0; i < j; ++i) { grad(gidx++) = R_bar(i, j); } - grad(gidx++) = R_bar(j, j) * temp_cholesky(j, j) + w_j; + grad(gidx++) = R_bar(j, j) * temp_cholesky(j, j) + w_j + + 2.0 * determinant_tilt_yy_; } return {logp, grad}; diff --git a/src/models/mixed/mixed_mrf_model.cpp b/src/models/mixed/mixed_mrf_model.cpp index 55676699..71cd12cc 100644 --- a/src/models/mixed/mixed_mrf_model.cpp +++ b/src/models/mixed/mixed_mrf_model.cpp @@ -124,6 +124,7 @@ MixedMRFModel::MixedMRFModel( MixedMRFModel::MixedMRFModel(const MixedMRFModel& other) : BaseModel(other), target_accept_(other.target_accept_), + determinant_tilt_yy_(other.determinant_tilt_yy_), n_(other.n_), p_(other.p_), q_(other.q_), diff --git a/src/models/mixed/mixed_mrf_model.h b/src/models/mixed/mixed_mrf_model.h index c86d3c04..b39bd80f 100644 --- a/src/models/mixed/mixed_mrf_model.h +++ b/src/models/mixed/mixed_mrf_model.h @@ -134,6 +134,17 @@ class MixedMRFModel : public BaseModel { target_accept_ = target; } + /** + * Set the determinant-tilt exponent delta for the Kyy block. Adds + * delta * log|Kyy| to the log-prior, pushing the continuous-block + * precision matrix away from the PD-cone boundary. delta = 0 + * (default) recovers the untilted target. Currently consumed only + * by the NUTS gradient paths; the MH path is unchanged. + */ + void set_determinant_tilt_yy(double delta) { + determinant_tilt_yy_ = delta; + } + /** * Construct Robbins-Monro adaptation controllers for the per-iteration * MH proposal SDs. Called once by MetropolisSampler before warmup; under @@ -302,6 +313,11 @@ class MixedMRFModel : public BaseModel { // to 0.44 (componentwise random-walk Metropolis optimum). double target_accept_ = 0.44; + // Determinant-tilt exponent on the Kyy block (see set_determinant_tilt_yy). + // Adds determinant_tilt_yy_ * log|Kyy| to the NUTS log-prior; MH ratios + // are not yet adjusted. + double determinant_tilt_yy_ = 0.0; + /// Per-iteration adaptation controllers (MH mode only — under NUTS these /// stay null and the stage-3b path in tune_proposal_sd is used instead). /// One adapter per proposal-SD storage; off-diag and diag of the continuous diff --git a/tests/testthat/test-mixed-gradient-pfaffian.R b/tests/testthat/test-mixed-gradient-pfaffian.R index d39a4247..015fe76c 100644 --- a/tests/testthat/test-mixed-gradient-pfaffian.R +++ b/tests/testthat/test-mixed-gradient-pfaffian.R @@ -19,7 +19,7 @@ mixed_full_fd_gradient = function(params, x, y, num_cats, is_ord, base_cat, edge_ind, scale, inv_mass = NULL, - eps = 1e-5) { + delta = 0, eps = 1e-5) { n_total = length(params) fd = numeric(n_total) for(k in seq_len(n_total)) { @@ -30,12 +30,12 @@ mixed_full_fd_gradient = function(params, x, y, num_cats, is_ord, base_cat, fp = mixed_test_logp_and_gradient_full( p_plus, x, y, num_cats, as.integer(is_ord), base_cat, edge_ind, scale, - inv_mass_diag = inv_mass + inv_mass_diag = inv_mass, delta = delta )$value fm = mixed_test_logp_and_gradient_full( p_minus, x, y, num_cats, as.integer(is_ord), base_cat, edge_ind, scale, - inv_mass_diag = inv_mass + inv_mass_diag = inv_mass, delta = delta )$value fd[k] = (fp - fm) / (2 * eps) } @@ -44,15 +44,15 @@ mixed_full_fd_gradient = function(params, x, y, num_cats, is_ord, base_cat, mixed_full_check_gradient = function(params, x, y, num_cats, is_ord, base_cat, edge_ind, scale, inv_mass = NULL, - eps = 1e-5) { + delta = 0, eps = 1e-5) { res = mixed_test_logp_and_gradient_full( params, x, y, num_cats, as.integer(is_ord), base_cat, edge_ind, scale, - inv_mass_diag = inv_mass + inv_mass_diag = inv_mass, delta = delta ) fd = mixed_full_fd_gradient( params, x, y, num_cats, is_ord, base_cat, - edge_ind, scale, inv_mass, eps + edge_ind, scale, inv_mass, delta, eps ) ag = res$gradient denom = pmax(abs(ag), abs(fd), 1) @@ -148,3 +148,72 @@ test_that("gradient matches FD with two excluded Kyy edges (random mass)", { ) expect_lt(err, 1e-5) }) + + +# ---- Determinant tilt on the Kyy block -------------------------------------- # + +test_that("delta = 0 is a no-op (bit-identical to no-tilt call)", { + s = mixed_full_make_setup() + ref = mixed_test_logp_and_gradient_full( + s$params, s$x, s$y, s$num_cats, as.integer(s$is_ord), + s$base_cat, s$edge_ind, 2.5 + ) + with_zero = mixed_test_logp_and_gradient_full( + s$params, s$x, s$y, s$num_cats, as.integer(s$is_ord), + s$base_cat, s$edge_ind, 2.5, + delta = 0 + ) + expect_equal(with_zero$value, ref$value) + expect_equal(with_zero$gradient, ref$gradient) +}) + + +test_that("gradient matches FD with delta > 0 on full Kyy", { + s = mixed_full_make_setup() + err = mixed_full_check_gradient( + s$params, s$x, s$y, s$num_cats, s$is_ord, s$base_cat, + s$edge_ind, 2.5, + delta = 1.5 + ) + expect_lt(err, 1e-5) +}) + + +test_that("gradient matches FD with delta > 0 on sparse Kyy (random mass)", { + s = mixed_full_make_setup() + s$edge_ind[s$p + 1L, s$p + 3L] = 0L + s$edge_ind[s$p + 3L, s$p + 1L] = 0L + set.seed(21) + inv_mass = runif(s$total_dim, min = 0.5, max = 2.0) + err = mixed_full_check_gradient( + s$params, s$x, s$y, s$num_cats, s$is_ord, s$base_cat, + s$edge_ind, 2.5, + inv_mass = inv_mass, delta = 2.0 + ) + expect_lt(err, 1e-5) +}) + + +test_that("delta shifts logp by exactly delta * log|Kyy|", { + # Analytic check: tilting by |Kyy|^delta adds delta * log|Kyy| to logp. + s = mixed_full_make_setup() + base = mixed_test_logp_and_gradient_full( + s$params, s$x, s$y, s$num_cats, as.integer(s$is_ord), + s$base_cat, s$edge_ind, 2.5, delta = 0 + )$value + tilted = mixed_test_logp_and_gradient_full( + s$params, s$x, s$y, s$num_cats, as.integer(s$is_ord), + s$base_cat, s$edge_ind, 2.5, delta = 3.0 + )$value + + # log|Kyy| from psi entries (Cholesky-of-precision diagonals are exp(psi)). + n_main = sum(s$num_cats) + n_pw_xx = s$p * (s$p - 1L) / 2L + n_means = s$q + n_xy = s$p * s$q + chol_offset = n_main + n_pw_xx + n_means + n_xy + psi_positions = chol_offset + cumsum(seq_len(s$q)) + log_det_Kyy = 2 * sum(s$params[psi_positions]) + + expect_equal(tilted - base, 3.0 * log_det_Kyy, tolerance = 1e-12) +}) From 68ef39c2310856b089b4bf1cef0004889c52bea5 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Fri, 15 May 2026 09:15:48 +0200 Subject: [PATCH 04/10] feat(ggm): determinant-tilt prior for the MH path Companion to commit e4f24cb (NUTS-side GGM tilt). Adds the matching delta * (log|K_prop| - log|K_curr|) term to all six GGMModel MH-ratio sites: - update_edge_parameter - update_diagonal_parameter - update_edge_indicator_parameter_pair (both branches) - tune_proposal_sd (off-diagonal and diagonal passes) The log-determinant ratio is computed in O(p) for rank-2 edge updates and O(1) for rank-1 diagonal updates via the matrix-determinant lemma, reusing the cached covariance_matrix_ and the already-set-up precision_proposal_. Existing log_density_impl_edge/diag are refactored to share two new private helpers log_det_ratio_edge/diag; both helpers return exactly the logdet term that those functions previously computed inline, so the refactor is observationally inert at delta = 0. Each tilt addition is guarded by `if (determinant_tilt_ != 0.0)` so the existing default code path remains a single-pass MH ratio with no additional log() call. NUTS and MH paths now both consume determinant_tilt_ consistently; the chain targets the same posterior under either sampler. No regressions in the existing test suite (6708 pass, 0 fail). --- src/models/ggm/ggm_model.cpp | 86 ++++++++++++++++++++++++++---------- src/models/ggm/ggm_model.h | 17 +++++++ 2 files changed, 80 insertions(+), 23 deletions(-) diff --git a/src/models/ggm/ggm_model.cpp b/src/models/ggm/ggm_model.cpp index 8b9496f7..6e440e44 100644 --- a/src/models/ggm/ggm_model.cpp +++ b/src/models/ggm/ggm_model.cpp @@ -610,41 +610,49 @@ double GGMModel::log_density_impl(const arma::mat& omega, const arma::mat& phi) return log_likelihood; } -double GGMModel::log_density_impl_edge(size_t i, size_t j) const { - - // Log-likelihood ratio (not the full log-likelihood) - +double GGMModel::log_det_ratio_edge(size_t i, size_t j) const { + // Rank-2 matrix-determinant lemma: log|K_prop| - log|K_curr| where K_prop + // differs from K_curr at entries (i,j), (j,i), and (j,j). cc11, cc12, cc22 + // are the entries of the 2x2 update Gram matrix I + V^T Sigma U used by + // the lemma; |det(...)| is its determinant. double Ui2 = precision_matrix_(i, j) - precision_proposal_(i, j); double Uj2 = (precision_matrix_(j, j) - precision_proposal_(j, j)) / 2; - double cc11 = 0 + covariance_matrix_(j, j); - double cc12 = 1 - (covariance_matrix_(i, j) * Ui2 + covariance_matrix_(j, j) * Uj2); - double cc22 = 0 + Ui2 * Ui2 * covariance_matrix_(i, i) + 2 * Ui2 * Uj2 * covariance_matrix_(i, j) + Uj2 * Uj2 * covariance_matrix_(j, j); - - double logdet = MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); - // logdet - (logdet(aOmega_prop) - logdet(aOmega)) - - double trace_prod = -2 * (suf_stat_(j, j) * Uj2 + suf_stat_(i, j) * Ui2); - - double log_likelihood_ratio = (n_ * logdet - trace_prod) / 2; - return log_likelihood_ratio; + double cc11 = covariance_matrix_(j, j); + double cc12 = 1 - (covariance_matrix_(i, j) * Ui2 + + covariance_matrix_(j, j) * Uj2); + double cc22 = Ui2 * Ui2 * covariance_matrix_(i, i) + + 2 * Ui2 * Uj2 * covariance_matrix_(i, j) + + Uj2 * Uj2 * covariance_matrix_(j, j); + return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); } -double GGMModel::log_density_impl_diag(size_t j) const { - // same as above but for i == j, so Ui2 = 0 +double GGMModel::log_det_ratio_diag(size_t j) const { + // Rank-1 specialisation of log_det_ratio_edge (Ui2 = 0). double Uj2 = (precision_matrix_(j, j) - precision_proposal_(j, j)) / 2; - double cc11 = 0 + covariance_matrix_(j, j); + double cc11 = covariance_matrix_(j, j); double cc12 = 1 - covariance_matrix_(j, j) * Uj2; - double cc22 = 0 + Uj2 * Uj2 * covariance_matrix_(j, j); + double cc22 = Uj2 * Uj2 * covariance_matrix_(j, j); - double logdet = MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); - double trace_prod = -2 * suf_stat_(j, j) * Uj2; + return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); +} - double log_likelihood_ratio = (n_ * logdet - trace_prod) / 2; - return log_likelihood_ratio; +double GGMModel::log_density_impl_edge(size_t i, size_t j) const { + // Log-likelihood ratio (not the full log-likelihood). + double Ui2 = precision_matrix_(i, j) - precision_proposal_(i, j); + double Uj2 = (precision_matrix_(j, j) - precision_proposal_(j, j)) / 2; + double logdet = log_det_ratio_edge(i, j); + double trace_prod = -2 * (suf_stat_(j, j) * Uj2 + suf_stat_(i, j) * Ui2); + return (n_ * logdet - trace_prod) / 2; +} +double GGMModel::log_density_impl_diag(size_t j) const { + double Uj2 = (precision_matrix_(j, j) - precision_proposal_(j, j)) / 2; + double logdet = log_det_ratio_diag(j); + double trace_prod = -2 * suf_stat_(j, j) * Uj2; + return (n_ * logdet - trace_prod) / 2; } double GGMModel::update_edge_parameter(size_t i, size_t j) { @@ -672,6 +680,14 @@ double GGMModel::update_edge_parameter(size_t i, size_t j) { double ln_alpha = log_density_impl_edge(i, j); + // Determinant-tilt prior: |K|^delta contributes + // delta * (log|K_prop| - log|K_curr|) + // to the MH ratio. log_det_ratio_edge uses the rank-2 matrix-determinant + // lemma in O(p) via the cached covariance, so this is essentially free. + if (determinant_tilt_ != 0.0) { + ln_alpha += determinant_tilt_ * log_det_ratio_edge(i, j); + } + // Interaction prior on K_yy_{ij} = -0.5 * Omega_{ij} ln_alpha += interaction_prior_->logp(-0.5 * precision_proposal_(i, j)); ln_alpha -= interaction_prior_->logp(-0.5 * precision_matrix_(i, j)); @@ -750,6 +766,12 @@ double GGMModel::update_diagonal_parameter(size_t i) { double ln_alpha = log_density_impl_diag(i); + // Determinant-tilt prior: |K|^delta contributes delta * log_det_ratio + // to the MH ratio. Rank-1 update => O(1) via the cached covariance. + if (determinant_tilt_ != 0.0) { + ln_alpha += determinant_tilt_ * log_det_ratio_diag(i); + } + ln_alpha += diagonal_prior_->logp(0.5 * precision_proposal_(i, i)); ln_alpha -= diagonal_prior_->logp(0.5 * precision_matrix_(i, i)); ln_alpha += 2.0 * (theta_prop - theta_curr); // Jacobian: dK_ii/dtheta = 2*exp(2*theta) @@ -818,6 +840,11 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { // } // } + // Determinant-tilt prior: |K|^delta contributes delta * log_det_ratio + // to the MH ratio. The rank-2 update at (i,j),(j,j) makes this O(p). + if (determinant_tilt_ != 0.0) { + ln_alpha += determinant_tilt_ * log_det_ratio_edge(i, j); + } ln_alpha += MY_LOG(1.0 - inclusion_probability_(i, j)) - MY_LOG(inclusion_probability_(i, j)); @@ -877,6 +904,13 @@ void GGMModel::update_edge_indicator_parameter_pair(size_t i, size_t j) { // Rcpp::Rcout << "ln_alpha: " << ln_alpha << ", ln_alpha_ref: " << ln_alpha_ref << std::endl; // } // } + + // Determinant-tilt prior: |K|^delta contributes delta * log_det_ratio + // to the MH ratio. + if (determinant_tilt_ != 0.0) { + ln_alpha += determinant_tilt_ * log_det_ratio_edge(i, j); + } + ln_alpha += MY_LOG(inclusion_probability_(i, j)) - MY_LOG(1.0 - inclusion_probability_(i, j)); // Slab in K_yy coords; proposal in K_ij coords. Jacobian |dK_yy/dK_ij| = 1/2. @@ -1004,6 +1038,9 @@ void GGMModel::tune_proposal_sd(int iteration, const WarmupSchedule& schedule) { precision_proposal_(j, j) = omega_prop_qq; double ln_alpha = log_density_impl_edge(i, j); + if (determinant_tilt_ != 0.0) { + ln_alpha += determinant_tilt_ * log_det_ratio_edge(i, j); + } // Interaction prior on K_yy_{ij} = -0.5 * Omega_{ij} ln_alpha += interaction_prior_->logp(-0.5 * precision_proposal_(i, j)); ln_alpha -= interaction_prior_->logp(-0.5 * precision_matrix_(i, j)); @@ -1045,6 +1082,9 @@ void GGMModel::tune_proposal_sd(int iteration, const WarmupSchedule& schedule) { + MY_EXP(theta_prop) * MY_EXP(theta_prop); double ln_alpha = log_density_impl_diag(i); + if (determinant_tilt_ != 0.0) { + ln_alpha += determinant_tilt_ * log_det_ratio_diag(i); + } ln_alpha += diagonal_prior_->logp(0.5 * precision_proposal_(i, i)); ln_alpha -= diagonal_prior_->logp(0.5 * precision_matrix_(i, i)); ln_alpha += 2.0 * (theta_prop - theta_curr); // Jacobian: dK_ii/dtheta = 2*exp(2*theta) diff --git a/src/models/ggm/ggm_model.h b/src/models/ggm/ggm_model.h index f0cc5c5d..0e81bf99 100644 --- a/src/models/ggm/ggm_model.h +++ b/src/models/ggm/ggm_model.h @@ -674,6 +674,23 @@ class GGMModel : public BaseModel { */ double log_density_impl_diag(size_t j) const; + /** + * log|K_prop| - log|K_curr| for a rank-2 off-diagonal proposal at (i, j), + * computed via the matrix-determinant lemma in O(p). Reads + * precision_matrix_, precision_proposal_, and covariance_matrix_; assumes + * precision_proposal_ has already been set up at (i, j), (j, i), (j, j). + * Used to add the determinant-tilt term delta * (log|K_prop| - log|K_curr|) + * to MH ratios. + */ + double log_det_ratio_edge(size_t i, size_t j) const; + + /** + * log|K_prop| - log|K_curr| for a rank-1 diagonal proposal at j. + * Computed via the matrix-determinant lemma in O(1). Reads the same + * cached state as log_det_ratio_edge. + */ + double log_det_ratio_diag(size_t j) const; + /** From 083b5fae1c3bcd2ec55d22a0a666f2af145c73fc Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Fri, 15 May 2026 09:21:09 +0200 Subject: [PATCH 05/10] feat(mixed-mrf): determinant-tilt prior for the Kyy MH path Companion to commit 1bb74c7 (NUTS-side mixed-MRF tilt). Adds the matching delta * (log|Kyy_prop| - log|Kyy_curr|) term to the three continuous-block MH-ratio sites in MixedMRFModel: - update_pairwise_effects_continuous_offdiag - update_pairwise_effects_continuous_diag - update_edge_indicator_continuous Adds two private helpers log_det_ratio_yy_edge/diag that compute the log-determinant ratio via the matrix-determinant lemma in O(q) for rank-2 and O(1) for rank-1, reusing covariance_continuous_ and the already-set-up precision_proposal_. The helpers mirror GGMModel::log_det_ratio_edge/diag. Each tilt addition is guarded by `if (determinant_tilt_yy_ != 0.0)` so the default delta_yy = 0 path is bit-identical to the previous code. NUTS and MH now both consume determinant_tilt_yy_ consistently; the chain targets the same posterior under either sampler. No regressions in the existing test suite (6708 pass, 0 fail). --- src/models/mixed/mixed_mrf_metropolis.cpp | 52 +++++++++++++++++++++++ src/models/mixed/mixed_mrf_model.h | 13 ++++++ 2 files changed, 65 insertions(+) diff --git a/src/models/mixed/mixed_mrf_metropolis.cpp b/src/models/mixed/mixed_mrf_metropolis.cpp index d5c0e75d..50f49b7c 100644 --- a/src/models/mixed/mixed_mrf_metropolis.cpp +++ b/src/models/mixed/mixed_mrf_metropolis.cpp @@ -187,6 +187,41 @@ double MixedMRFModel::precision_constrained_diagonal(double x) const { // an O(nq) rank-2 shortcut. // ============================================================================= +double MixedMRFModel::log_det_ratio_yy_edge(int i, int j) const { + // Rank-2 matrix-determinant lemma applied to the Kyy block. Mirrors + // GGMModel::log_det_ratio_edge but uses pairwise_effects_continuous_ + // (Kyy = -2 * pairwise_effects_continuous_) and covariance_continuous_. + size_t ui = static_cast(i); + size_t uj = static_cast(j); + double precision_ij_curr = -2.0 * pairwise_effects_continuous_(ui, uj); + double precision_jj_curr = -2.0 * pairwise_effects_continuous_(uj, uj); + double Ui = precision_ij_curr - precision_proposal_(ui, uj); + double Uj = (precision_jj_curr - precision_proposal_(uj, uj)) / 2.0; + + double cc11 = covariance_continuous_(uj, uj); + double cc12 = 1.0 - (covariance_continuous_(ui, uj) * Ui + + covariance_continuous_(uj, uj) * Uj); + double cc22 = Ui * Ui * covariance_continuous_(ui, ui) + + 2.0 * Ui * Uj * covariance_continuous_(ui, uj) + + Uj * Uj * covariance_continuous_(uj, uj); + + return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); +} + +double MixedMRFModel::log_det_ratio_yy_diag(int i) const { + // Rank-1 specialisation of log_det_ratio_yy_edge (Ui = 0). + size_t ui = static_cast(i); + double precision_ii_curr = -2.0 * pairwise_effects_continuous_(ui, ui); + double Uj = (precision_ii_curr - precision_proposal_(ui, ui)) / 2.0; + + double cc11 = covariance_continuous_(ui, ui); + double cc12 = 1.0 - covariance_continuous_(ui, ui) * Uj; + double cc22 = Uj * Uj * covariance_continuous_(ui, ui); + + return MY_LOG(std::abs(cc11 * cc22 - cc12 * cc12)); +} + + double MixedMRFModel::log_ggm_ratio_edge(int i, int j, arma::mat& cov_prop_out) const { size_t ui = static_cast(i); size_t uj = static_cast(j); @@ -403,6 +438,13 @@ double MixedMRFModel::update_pairwise_effects_continuous_offdiag(int i, int j, s arma::mat cov_prop; double ln_alpha = log_ggm_ratio_edge(i, j, cov_prop); + // Determinant-tilt prior on the Kyy block: |Kyy|^delta contributes + // delta * (log|Kyy_prop| - log|Kyy_curr|) + // to the MH ratio. Rank-2 lemma, O(q) via the cached covariance. + if (determinant_tilt_yy_ != 0.0) { + ln_alpha += determinant_tilt_yy_ * log_det_ratio_yy_edge(i, j); + } + // OMRF ratio with proposed continuous interactions. The marginal // interactions M = A_xx + 2 A_xy Σ A_xy^T depend on Σ; when Kyy changes, // Σ changes too, so we must use Σ' (cov_prop) for the proposed-state @@ -490,6 +532,11 @@ double MixedMRFModel::update_pairwise_effects_continuous_diag(int i, std::option arma::mat cov_prop; double ln_alpha = log_ggm_ratio_diag(i, cov_prop); + // Determinant-tilt prior: rank-1 lemma, O(1) via the cached covariance. + if (determinant_tilt_yy_ != 0.0) { + ln_alpha += determinant_tilt_yy_ * log_det_ratio_yy_diag(i); + } + // OMRF ratio with proposed continuous interactions. Use Σ' (cov_prop) for // the proposed-state marginal_interactions, not the cached Σ. for(size_t s = 0; s < p_; ++s) @@ -683,6 +730,11 @@ void MixedMRFModel::update_edge_indicator_continuous(int i, int j) { arma::mat cov_prop; double ln_alpha = log_ggm_ratio_edge(i, j, cov_prop); + // Determinant-tilt prior: see update_pairwise_effects_continuous_offdiag. + if (determinant_tilt_yy_ != 0.0) { + ln_alpha += determinant_tilt_yy_ * log_det_ratio_yy_edge(i, j); + } + // OMRF ratio with proposed continuous interactions. Use Σ' (cov_prop) for // the proposed-state marginal_interactions, not the cached Σ. for(size_t s = 0; s < p_; ++s) diff --git a/src/models/mixed/mixed_mrf_model.h b/src/models/mixed/mixed_mrf_model.h index b39bd80f..9aa56522 100644 --- a/src/models/mixed/mixed_mrf_model.h +++ b/src/models/mixed/mixed_mrf_model.h @@ -547,6 +547,19 @@ class MixedMRFModel : public BaseModel { // proposed covariance Σ' (computed via Sherman-Morrison) to cov_prop_out. double log_ggm_ratio_diag(int i, arma::mat& cov_prop_out) const; + // log|Kyy_prop| - log|Kyy_curr| for a rank-2 off-diagonal proposal at + // (i, j), via the matrix-determinant lemma in O(q). Reads + // pairwise_effects_continuous_, precision_proposal_, and + // covariance_continuous_; assumes precision_proposal_ has the proposed + // Kyy at (i, j), (j, i), (j, j) already filled. Used to add the + // determinant-tilt term delta_yy * (log|Kyy_prop| - log|Kyy_curr|) to MH + // ratios. + double log_det_ratio_yy_edge(int i, int j) const; + + // log|Kyy_prop| - log|Kyy_curr| for a rank-1 diagonal proposal at i. + // Computed via the matrix-determinant lemma in O(1). + double log_det_ratio_yy_diag(int i) const; + // Rank-1 Cholesky update after accepting an off-diagonal precision change. void cholesky_update_after_precision_edge(double old_ij, double old_jj, int i, int j); From 8a9f7cbb6c8e3c8585f3ddba2fe2bb73322cae7e Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Fri, 15 May 2026 09:37:27 +0200 Subject: [PATCH 06/10] feat(bgm): expose determinant-tilt prior delta on bgm() Adds a delta argument to bgm() that propagates the determinant tilt to both the NUTS and adaptive-Metropolis paths for GGM and mixed-MRF models. delta = 0 is the default (untilted) and reproduces the current behaviour bit-for-bit. R-layer plumbing: bgm(delta = 0) -> bgm_spec(delta = 0) # validation + dispatch -> build_spec_{ggm,mixed_mrf} # stored on spec$prior$delta -> run_sampler_{ggm,mixed_mrf} # forwards to C++ entry point -> sample_{ggm,mixed_mrf}(delta = ...) C++-layer plumbing: sample_ggm calls model.set_determinant_tilt(delta) sample_mixed_mrf calls model.set_determinant_tilt_yy(delta) bgm_spec rejects delta > 0 for pure-ordinal models (model_type "omrf" or "compare") since there's no precision matrix to tilt. Negative, non-finite, or non-scalar delta values produce a clear error early in bgm_spec(). New test file test-bgm-delta.R covers: - NUTS path: delta > 0 shifts tr(K) upward vs delta = 0 (GGM) - MH path: delta > 0 shifts tr(K) upward vs delta = 0 (GGM) - Pure-ordinal model rejects delta > 0 with informative message - Invalid delta (negative, NA, non-scalar) rejected at the validator RcppExports regenerated via Rcpp::compileAttributes(). Full suite: 6714 pass, 0 fail. --- R/RcppExports.R | 8 +-- R/bgm.R | 10 +++ R/bgm_spec.R | 20 ++++++ R/run_sampler.R | 6 +- src/RcppExports.cpp | 18 ++--- src/sample_ggm.cpp | 8 ++- src/sample_mixed.cpp | 9 ++- tests/testthat/test-bgm-delta.R | 112 ++++++++++++++++++++++++++++++++ 8 files changed, 175 insertions(+), 16 deletions(-) create mode 100644 tests/testthat/test-bgm-delta.R diff --git a/R/RcppExports.R b/R/RcppExports.R index 3a946468..b6327205 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -141,12 +141,12 @@ ggm_test_logp_and_gradient_full_prior <- function(x, suf_stat, n, edge_indicator .Call(`_bgms_ggm_test_logp_and_gradient_full_prior`, x, suf_stat, n, edge_indicators, interaction_prior_type, interaction_scale, interaction_alpha, interaction_beta, diagonal_prior_type, diagonal_shape, diagonal_rate, inv_mass_diag) } -sample_ggm <- function(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback = NULL, edge_prior = "Bernoulli", beta_bernoulli_alpha = 1.0, beta_bernoulli_beta = 1.0, beta_bernoulli_alpha_between = 1.0, beta_bernoulli_beta_between = 1.0, dirichlet_alpha = 1.0, lambda = 1.0, target_acceptance = 0.8, max_tree_depth = 10L, na_impute = FALSE, missing_index_nullable = NULL) { - .Call(`_bgms_sample_ggm`, inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback, edge_prior, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, target_acceptance, max_tree_depth, na_impute, missing_index_nullable) +sample_ggm <- function(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback = NULL, edge_prior = "Bernoulli", beta_bernoulli_alpha = 1.0, beta_bernoulli_beta = 1.0, beta_bernoulli_alpha_between = 1.0, beta_bernoulli_beta_between = 1.0, dirichlet_alpha = 1.0, lambda = 1.0, target_acceptance = 0.8, max_tree_depth = 10L, na_impute = FALSE, missing_index_nullable = NULL, delta = 0.0) { + .Call(`_bgms_sample_ggm`, inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback, edge_prior, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, target_acceptance, max_tree_depth, na_impute, missing_index_nullable, delta) } -sample_mixed_mrf <- function(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type, progress_callback = NULL, edge_prior = "Bernoulli", beta_bernoulli_alpha = 1.0, beta_bernoulli_beta = 1.0, beta_bernoulli_alpha_between = 1.0, beta_bernoulli_beta_between = 1.0, dirichlet_alpha = 1.0, lambda = 1.0, sampler_type = "mh", target_acceptance = 0.80, max_tree_depth = 10L, na_impute = FALSE, missing_index_discrete_nullable = NULL, missing_index_continuous_nullable = NULL) { - .Call(`_bgms_sample_mixed_mrf`, inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type, progress_callback, edge_prior, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, sampler_type, target_acceptance, max_tree_depth, na_impute, missing_index_discrete_nullable, missing_index_continuous_nullable) +sample_mixed_mrf <- function(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type, progress_callback = NULL, edge_prior = "Bernoulli", beta_bernoulli_alpha = 1.0, beta_bernoulli_beta = 1.0, beta_bernoulli_alpha_between = 1.0, beta_bernoulli_beta_between = 1.0, dirichlet_alpha = 1.0, lambda = 1.0, sampler_type = "mh", target_acceptance = 0.80, max_tree_depth = 10L, na_impute = FALSE, missing_index_discrete_nullable = NULL, missing_index_continuous_nullable = NULL, delta = 0.0) { + .Call(`_bgms_sample_mixed_mrf`, inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type, progress_callback, edge_prior, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, sampler_type, target_acceptance, max_tree_depth, na_impute, missing_index_discrete_nullable, missing_index_continuous_nullable, delta) } sample_omrf <- function(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback = NULL, edge_prior = "Bernoulli", na_impute = FALSE, missing_index_nullable = NULL, beta_bernoulli_alpha = 1.0, beta_bernoulli_beta = 1.0, beta_bernoulli_alpha_between = 1.0, beta_bernoulli_beta_between = 1.0, dirichlet_alpha = 1.0, lambda = 1.0, target_acceptance = 0.8, max_tree_depth = 10L, pairwise_scaling_factors_nullable = NULL) { diff --git a/R/bgm.R b/R/bgm.R index 9d6ef990..e89b63b6 100644 --- a/R/bgm.R +++ b/R/bgm.R @@ -94,6 +94,14 @@ #' Ignored for pure ordinal models. #' Default: \code{gamma_prior(shape = 1, rate = 1)}. #' +#' @param delta Non-negative numeric. Determinant-tilt exponent on the +#' continuous-block precision matrix \eqn{K} (GGM) or \eqn{K_{yy}} +#' (mixed MRF): multiplies the prior by \eqn{|K|^{\delta}}, pushing the +#' chain away from the positive-definite cone boundary. Both NUTS and +#' adaptive-Metropolis update paths apply the tilt. \code{delta = 0} +#' (default) recovers the untilted prior. Not allowed for pure ordinal +#' models (no precision matrix to tilt). +#' #' @param pairwise_scale `r lifecycle::badge("deprecated")` Double. #' Scale of the Cauchy prior for pairwise #' interaction parameters. Use \code{interaction_prior} instead. @@ -325,6 +333,7 @@ bgm = function( threshold_prior = beta_prime_prior(alpha = 0.5, beta = 0.5), means_prior = normal_prior(scale = 1), precision_scale_prior = gamma_prior(shape = 1, rate = 1), + delta = 0, edge_selection = TRUE, edge_prior = bernoulli_prior(0.5), na_action = c("listwise", "impute"), @@ -492,6 +501,7 @@ bgm = function( scale_prior_type = sp$scale_prior_type, scale_shape = sp$scale_shape, scale_rate = sp$scale_rate, + delta = delta, standardize = standardize, edge_selection = edge_selection, edge_prior = edge_prior, diff --git a/R/bgm_spec.R b/R/bgm_spec.R index 2199024b..d99a28ed 100644 --- a/R/bgm_spec.R +++ b/R/bgm_spec.R @@ -265,6 +265,7 @@ bgm_spec = function(x, scale_prior_type = "gamma", scale_shape = 1, scale_rate = 1, + delta = 0, standardize = FALSE, edge_selection = TRUE, edge_prior = bernoulli_prior(0.5), @@ -344,6 +345,19 @@ bgm_spec = function(x, model_type = "mixed_mrf" } + # Validate determinant-tilt exponent and reject for pure-ordinal models + if(!is.numeric(delta) || length(delta) != 1L || is.na(delta) || + !is.finite(delta) || delta < 0) { + stop("'delta' must be a single finite non-negative numeric.") + } + if(delta > 0 && model_type %in% c("omrf", "compare")) { + stop( + "'delta' (determinant tilt) requires continuous variables; the ", + "current model_type is '", model_type, "', which has no precision ", + "matrix to tilt. Pass delta = 0 (default) or use continuous data." + ) + } + # --- Sampler (needs is_continuous and edge_selection early) ------------------ sampler = validate_sampler( update_method = update_method, @@ -415,6 +429,7 @@ bgm_spec = function(x, scale_prior_type = scale_prior_type, scale_shape = scale_shape, scale_rate = scale_rate, + delta = delta, edge_prior_flat = ep_flat ) } else if(model_type == "mixed_mrf") { @@ -438,6 +453,7 @@ bgm_spec = function(x, scale_prior_type = scale_prior_type, scale_shape = scale_shape, scale_rate = scale_rate, + delta = delta, standardize = standardize, edge_prior_flat = ep_flat ) @@ -505,6 +521,7 @@ build_spec_ggm = function(x, data_columnnames, num_variables, interaction_prior_type, pairwise_scale, interaction_alpha, interaction_beta, scale_prior_type, scale_shape, scale_rate, + delta = 0, edge_prior_flat) { # Missing data md = validate_missing_data( @@ -545,6 +562,7 @@ build_spec_ggm = function(x, data_columnnames, num_variables, scale_prior_type = scale_prior_type, scale_shape = scale_shape, scale_rate = scale_rate, + delta = delta, edge_selection = ep$edge_selection, edge_prior = ep$edge_prior, inclusion_probability = ep$inclusion_probability, @@ -683,6 +701,7 @@ build_spec_mixed_mrf = function(x, data_columnnames, num_variables, means_prior_type, means_scale, means_alpha, means_beta, scale_prior_type, scale_shape, scale_rate, + delta = 0, standardize, edge_prior_flat) { # Identify discrete vs continuous columns @@ -818,6 +837,7 @@ build_spec_mixed_mrf = function(x, data_columnnames, num_variables, scale_prior_type = scale_prior_type, scale_shape = scale_shape, scale_rate = scale_rate, + delta = delta, standardize = standardize, edge_selection = ep$edge_selection, edge_prior = ep$edge_prior, diff --git a/R/run_sampler.R b/R/run_sampler.R index d893818f..53af96b0 100644 --- a/R/run_sampler.R +++ b/R/run_sampler.R @@ -104,7 +104,8 @@ run_sampler_ggm = function(spec) { target_acceptance = s$target_accept, max_tree_depth = s$nuts_max_depth, na_impute = m$na_impute, - missing_index_nullable = m$missing_index + missing_index_nullable = m$missing_index, + delta = p$delta ) out_raw @@ -236,7 +237,8 @@ run_sampler_mixed_mrf = function(spec) { max_tree_depth = s$nuts_max_depth, na_impute = m$na_impute, missing_index_discrete_nullable = m$missing_index_discrete, - missing_index_continuous_nullable = m$missing_index_continuous + missing_index_continuous_nullable = m$missing_index_continuous, + delta = p$delta ) out_raw diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 7bc2ba9f..23dba4e2 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -672,8 +672,8 @@ BEGIN_RCPP END_RCPP } // sample_ggm -Rcpp::List sample_ggm(const Rcpp::List& inputFromR, const arma::mat& prior_inclusion_prob, const arma::imat& initial_edge_indicators, const int no_iter, const int no_warmup, const int no_chains, const bool edge_selection, const std::string& sampler_type, const int seed, const int no_threads, const int progress_type, SEXP progress_callback, const std::string& edge_prior, const double beta_bernoulli_alpha, const double beta_bernoulli_beta, const double beta_bernoulli_alpha_between, const double beta_bernoulli_beta_between, const double dirichlet_alpha, const double lambda, const double target_acceptance, const int max_tree_depth, const bool na_impute, const Rcpp::Nullable missing_index_nullable); -RcppExport SEXP _bgms_sample_ggm(SEXP inputFromRSEXP, SEXP prior_inclusion_probSEXP, SEXP initial_edge_indicatorsSEXP, SEXP no_iterSEXP, SEXP no_warmupSEXP, SEXP no_chainsSEXP, SEXP edge_selectionSEXP, SEXP sampler_typeSEXP, SEXP seedSEXP, SEXP no_threadsSEXP, SEXP progress_typeSEXP, SEXP progress_callbackSEXP, SEXP edge_priorSEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP beta_bernoulli_alpha_betweenSEXP, SEXP beta_bernoulli_beta_betweenSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP target_acceptanceSEXP, SEXP max_tree_depthSEXP, SEXP na_imputeSEXP, SEXP missing_index_nullableSEXP) { +Rcpp::List sample_ggm(const Rcpp::List& inputFromR, const arma::mat& prior_inclusion_prob, const arma::imat& initial_edge_indicators, const int no_iter, const int no_warmup, const int no_chains, const bool edge_selection, const std::string& sampler_type, const int seed, const int no_threads, const int progress_type, SEXP progress_callback, const std::string& edge_prior, const double beta_bernoulli_alpha, const double beta_bernoulli_beta, const double beta_bernoulli_alpha_between, const double beta_bernoulli_beta_between, const double dirichlet_alpha, const double lambda, const double target_acceptance, const int max_tree_depth, const bool na_impute, const Rcpp::Nullable missing_index_nullable, const double delta); +RcppExport SEXP _bgms_sample_ggm(SEXP inputFromRSEXP, SEXP prior_inclusion_probSEXP, SEXP initial_edge_indicatorsSEXP, SEXP no_iterSEXP, SEXP no_warmupSEXP, SEXP no_chainsSEXP, SEXP edge_selectionSEXP, SEXP sampler_typeSEXP, SEXP seedSEXP, SEXP no_threadsSEXP, SEXP progress_typeSEXP, SEXP progress_callbackSEXP, SEXP edge_priorSEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP beta_bernoulli_alpha_betweenSEXP, SEXP beta_bernoulli_beta_betweenSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP target_acceptanceSEXP, SEXP max_tree_depthSEXP, SEXP na_imputeSEXP, SEXP missing_index_nullableSEXP, SEXP deltaSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -700,13 +700,14 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const int >::type max_tree_depth(max_tree_depthSEXP); Rcpp::traits::input_parameter< const bool >::type na_impute(na_imputeSEXP); Rcpp::traits::input_parameter< const Rcpp::Nullable >::type missing_index_nullable(missing_index_nullableSEXP); - rcpp_result_gen = Rcpp::wrap(sample_ggm(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback, edge_prior, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, target_acceptance, max_tree_depth, na_impute, missing_index_nullable)); + Rcpp::traits::input_parameter< const double >::type delta(deltaSEXP); + rcpp_result_gen = Rcpp::wrap(sample_ggm(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, sampler_type, seed, no_threads, progress_type, progress_callback, edge_prior, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, target_acceptance, max_tree_depth, na_impute, missing_index_nullable, delta)); return rcpp_result_gen; END_RCPP } // sample_mixed_mrf -Rcpp::List sample_mixed_mrf(const Rcpp::List& inputFromR, const arma::mat& prior_inclusion_prob, const arma::imat& initial_edge_indicators, const int no_iter, const int no_warmup, const int no_chains, const bool edge_selection, const int seed, const int no_threads, const int progress_type, SEXP progress_callback, const std::string& edge_prior, const double beta_bernoulli_alpha, const double beta_bernoulli_beta, const double beta_bernoulli_alpha_between, const double beta_bernoulli_beta_between, const double dirichlet_alpha, const double lambda, const std::string& sampler_type, const double target_acceptance, const int max_tree_depth, const bool na_impute, const Rcpp::Nullable missing_index_discrete_nullable, const Rcpp::Nullable missing_index_continuous_nullable); -RcppExport SEXP _bgms_sample_mixed_mrf(SEXP inputFromRSEXP, SEXP prior_inclusion_probSEXP, SEXP initial_edge_indicatorsSEXP, SEXP no_iterSEXP, SEXP no_warmupSEXP, SEXP no_chainsSEXP, SEXP edge_selectionSEXP, SEXP seedSEXP, SEXP no_threadsSEXP, SEXP progress_typeSEXP, SEXP progress_callbackSEXP, SEXP edge_priorSEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP beta_bernoulli_alpha_betweenSEXP, SEXP beta_bernoulli_beta_betweenSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP sampler_typeSEXP, SEXP target_acceptanceSEXP, SEXP max_tree_depthSEXP, SEXP na_imputeSEXP, SEXP missing_index_discrete_nullableSEXP, SEXP missing_index_continuous_nullableSEXP) { +Rcpp::List sample_mixed_mrf(const Rcpp::List& inputFromR, const arma::mat& prior_inclusion_prob, const arma::imat& initial_edge_indicators, const int no_iter, const int no_warmup, const int no_chains, const bool edge_selection, const int seed, const int no_threads, const int progress_type, SEXP progress_callback, const std::string& edge_prior, const double beta_bernoulli_alpha, const double beta_bernoulli_beta, const double beta_bernoulli_alpha_between, const double beta_bernoulli_beta_between, const double dirichlet_alpha, const double lambda, const std::string& sampler_type, const double target_acceptance, const int max_tree_depth, const bool na_impute, const Rcpp::Nullable missing_index_discrete_nullable, const Rcpp::Nullable missing_index_continuous_nullable, const double delta); +RcppExport SEXP _bgms_sample_mixed_mrf(SEXP inputFromRSEXP, SEXP prior_inclusion_probSEXP, SEXP initial_edge_indicatorsSEXP, SEXP no_iterSEXP, SEXP no_warmupSEXP, SEXP no_chainsSEXP, SEXP edge_selectionSEXP, SEXP seedSEXP, SEXP no_threadsSEXP, SEXP progress_typeSEXP, SEXP progress_callbackSEXP, SEXP edge_priorSEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP beta_bernoulli_alpha_betweenSEXP, SEXP beta_bernoulli_beta_betweenSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP sampler_typeSEXP, SEXP target_acceptanceSEXP, SEXP max_tree_depthSEXP, SEXP na_imputeSEXP, SEXP missing_index_discrete_nullableSEXP, SEXP missing_index_continuous_nullableSEXP, SEXP deltaSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -734,7 +735,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< const bool >::type na_impute(na_imputeSEXP); Rcpp::traits::input_parameter< const Rcpp::Nullable >::type missing_index_discrete_nullable(missing_index_discrete_nullableSEXP); Rcpp::traits::input_parameter< const Rcpp::Nullable >::type missing_index_continuous_nullable(missing_index_continuous_nullableSEXP); - rcpp_result_gen = Rcpp::wrap(sample_mixed_mrf(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type, progress_callback, edge_prior, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, sampler_type, target_acceptance, max_tree_depth, na_impute, missing_index_discrete_nullable, missing_index_continuous_nullable)); + Rcpp::traits::input_parameter< const double >::type delta(deltaSEXP); + rcpp_result_gen = Rcpp::wrap(sample_mixed_mrf(inputFromR, prior_inclusion_prob, initial_edge_indicators, no_iter, no_warmup, no_chains, edge_selection, seed, no_threads, progress_type, progress_callback, edge_prior, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, sampler_type, target_acceptance, max_tree_depth, na_impute, missing_index_discrete_nullable, missing_index_continuous_nullable, delta)); return rcpp_result_gen; END_RCPP } @@ -823,8 +825,8 @@ static const R_CallMethodDef CallEntries[] = { {"_bgms_test_scale_prior", (DL_FUNC) &_bgms_test_scale_prior, 4}, {"_bgms_ggm_test_logp_and_gradient_prior", (DL_FUNC) &_bgms_ggm_test_logp_and_gradient_prior, 11}, {"_bgms_ggm_test_logp_and_gradient_full_prior", (DL_FUNC) &_bgms_ggm_test_logp_and_gradient_full_prior, 12}, - {"_bgms_sample_ggm", (DL_FUNC) &_bgms_sample_ggm, 23}, - {"_bgms_sample_mixed_mrf", (DL_FUNC) &_bgms_sample_mixed_mrf, 24}, + {"_bgms_sample_ggm", (DL_FUNC) &_bgms_sample_ggm, 24}, + {"_bgms_sample_mixed_mrf", (DL_FUNC) &_bgms_sample_mixed_mrf, 25}, {"_bgms_sample_omrf", (DL_FUNC) &_bgms_sample_omrf, 24}, {"_bgms_compute_Vn_mfm_sbm", (DL_FUNC) &_bgms_compute_Vn_mfm_sbm, 4}, {NULL, NULL, 0} diff --git a/src/sample_ggm.cpp b/src/sample_ggm.cpp index 9a7ab080..8a979d4d 100644 --- a/src/sample_ggm.cpp +++ b/src/sample_ggm.cpp @@ -37,7 +37,8 @@ Rcpp::List sample_ggm( const double target_acceptance = 0.8, const int max_tree_depth = 10, const bool na_impute = false, - const Rcpp::Nullable missing_index_nullable = R_NilValue + const Rcpp::Nullable missing_index_nullable = R_NilValue, + const double delta = 0.0 ) { // Create parameter priors from R input @@ -79,6 +80,11 @@ Rcpp::List sample_ggm( const double mh_target = (sampler_type == "nuts") ? 0.44 : target_acceptance; model.set_metropolis_target_accept(mh_target); + // Determinant-tilt prior on |K|: shifts both NUTS and MH targets by + // delta * log|K|. delta = 0 is the default (untilted). Consumed by + // both gradient paths and all four MH ratios in GGMModel. + model.set_determinant_tilt(delta); + // Set up missing data imputation (same pattern as OMRF) if (na_impute && missing_index_nullable.isNotNull()) { arma::imat missing_index = Rcpp::as( diff --git a/src/sample_mixed.cpp b/src/sample_mixed.cpp index 557caf17..2c7a11bf 100644 --- a/src/sample_mixed.cpp +++ b/src/sample_mixed.cpp @@ -75,7 +75,8 @@ Rcpp::List sample_mixed_mrf( const int max_tree_depth = 10, const bool na_impute = false, const Rcpp::Nullable missing_index_discrete_nullable = R_NilValue, - const Rcpp::Nullable missing_index_continuous_nullable = R_NilValue + const Rcpp::Nullable missing_index_continuous_nullable = R_NilValue, + const double delta = 0.0 ) { // Extract model inputs from R list arma::imat discrete_obs = Rcpp::as(inputFromR["discrete_observations"]); @@ -147,6 +148,12 @@ Rcpp::List sample_mixed_mrf( const double mh_target = (sampler_type == "nuts") ? 0.44 : target_acceptance; model.set_metropolis_target_accept(mh_target); + // Determinant-tilt prior on |Kyy|: shifts both NUTS and MH targets by + // delta * log|Kyy|. delta = 0 is the default (untilted). Consumed by + // both gradient paths and the continuous-block MH ratios in + // MixedMRFModel. + model.set_determinant_tilt_yy(delta); + // Set up missing data imputation if(na_impute) { arma::imat missing_disc, missing_cont; diff --git a/tests/testthat/test-bgm-delta.R b/tests/testthat/test-bgm-delta.R new file mode 100644 index 00000000..e1bbb1f3 --- /dev/null +++ b/tests/testthat/test-bgm-delta.R @@ -0,0 +1,112 @@ +# --------------------------------------------------------------------------- # +# End-to-end tests for the bgm(delta = ...) plumbing. +# +# Verifies that: +# - delta = 0 (default) reproduces the untilted posterior +# - delta > 0 shifts the K_ii posterior upward (away from the PD-cone +# boundary), per the manuscript prediction Pr[lambda_min < x] <= C x^(delta+1) +# - delta > 0 is rejected for pure-ordinal models (no precision matrix) +# - delta < 0 / NA / non-numeric is rejected +# +# Tests run with NUTS to keep cost low; the MH path is exercised in a separate +# block with adaptive-metropolis. +# --------------------------------------------------------------------------- # + +# Small continuous dataset shared across tests. +gen_continuous = function(n = 60, p = 3, seed = 1L) { + set.seed(seed) + matrix(rnorm(n * p), nrow = n, ncol = p) +} + + +# Sum of K_ii from a fitted bgms object. Larger sum = K diagonal mass +# concentrated further from zero (further from the PD-cone boundary). +trace_K = function(out) { + sum(diag(extract_precision(out))) +} + + +# ---- delta > 0 pushes K diagonal upward ------------------------------------ # + +test_that("delta > 0 shifts NUTS K_ii posterior mean upward (GGM)", { + x = gen_continuous(n = 80, p = 3, seed = 2L) + base = bgm( + x, variable_type = "continuous", + iter = 300L, warmup = 300L, chains = 1L, + edge_selection = FALSE, verbose = FALSE, seed = 2L, + delta = 0 + ) + tilted = bgm( + x, variable_type = "continuous", + iter = 300L, warmup = 300L, chains = 1L, + edge_selection = FALSE, verbose = FALSE, seed = 2L, + delta = 5 + ) + expect_gt(trace_K(tilted), trace_K(base)) +}) + + +test_that("delta > 0 shifts MH K_ii posterior mean upward (GGM)", { + x = gen_continuous(n = 80, p = 3, seed = 3L) + base = bgm( + x, variable_type = "continuous", + update_method = "adaptive-metropolis", + iter = 400L, warmup = 400L, chains = 1L, + edge_selection = FALSE, verbose = FALSE, seed = 3L, + delta = 0 + ) + tilted = bgm( + x, variable_type = "continuous", + update_method = "adaptive-metropolis", + iter = 400L, warmup = 400L, chains = 1L, + edge_selection = FALSE, verbose = FALSE, seed = 3L, + delta = 5 + ) + expect_gt(trace_K(tilted), trace_K(base)) +}) + + +# ---- Validation ------------------------------------------------------------ # + +test_that("delta > 0 is rejected for pure-ordinal models", { + set.seed(4) + x = matrix(sample(0:1, 40 * 3, replace = TRUE), nrow = 40, ncol = 3) + expect_error( + bgm( + x, variable_type = "ordinal", + iter = 20L, warmup = 20L, chains = 1L, + edge_selection = FALSE, verbose = FALSE, + delta = 1 + ), + "no precision matrix to tilt" + ) +}) + + +test_that("invalid delta values are rejected with a clear message", { + x = gen_continuous(n = 40, p = 2) + expect_error( + bgm( + x, variable_type = "continuous", + iter = 20L, warmup = 20L, chains = 1L, + verbose = FALSE, delta = -1 + ), + "non-negative" + ) + expect_error( + bgm( + x, variable_type = "continuous", + iter = 20L, warmup = 20L, chains = 1L, + verbose = FALSE, delta = NA_real_ + ), + "finite" + ) + expect_error( + bgm( + x, variable_type = "continuous", + iter = 20L, warmup = 20L, chains = 1L, + verbose = FALSE, delta = c(0, 1) + ), + "single" + ) +}) From 9980c269c6fc362c14360d38314ae29b05c123b8 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Sat, 16 May 2026 11:32:15 +0200 Subject: [PATCH 07/10] fix(ggm): silence -Wreorder-ctor by matching declaration order The copy constructor initialised pcg_lambda_cache_ before determinant_tilt_, but determinant_tilt_ is declared first in the class. C++ always initialises in declaration order, so the listing was misleading and triggered -Wreorder-ctor across every translation unit that includes ggm_model.h. --- src/models/ggm/ggm_model.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/models/ggm/ggm_model.h b/src/models/ggm/ggm_model.h index 0e81bf99..658e6aa0 100644 --- a/src/models/ggm/ggm_model.h +++ b/src/models/ggm/ggm_model.h @@ -152,8 +152,8 @@ class GGMModel : public BaseModel { constraint_dirty_(other.constraint_dirty_), theta_valid_(other.theta_valid_), theta_(other.theta_), - pcg_lambda_cache_(other.pcg_lambda_cache_), - determinant_tilt_(other.determinant_tilt_) + determinant_tilt_(other.determinant_tilt_), + pcg_lambda_cache_(other.pcg_lambda_cache_) {} /** @return true (GGM supports NUTS via free-element Cholesky gradient). */ From c8b08e987a3c072ca57890eeb22396d8e0184d11 Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Sat, 16 May 2026 11:57:28 +0200 Subject: [PATCH 08/10] fix(ggm): move determinant_tilt_ init to match its declaration position The previous reorder fix swapped the trailing pair but left determinant_tilt_ near the end of the init list. determinant_tilt_ is declared just after target_accept_ near the top of the private section, so the init list still triggered -Wreorder-ctor against theta_. Move determinant_tilt_ to immediately after target_accept_ to match. --- src/models/ggm/ggm_model.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/models/ggm/ggm_model.h b/src/models/ggm/ggm_model.h index 658e6aa0..98e85c50 100644 --- a/src/models/ggm/ggm_model.h +++ b/src/models/ggm/ggm_model.h @@ -123,6 +123,7 @@ class GGMModel : public BaseModel { GGMModel(const GGMModel& other) : BaseModel(other), target_accept_(other.target_accept_), + determinant_tilt_(other.determinant_tilt_), n_(other.n_), p_(other.p_), dim_(other.dim_), @@ -152,7 +153,6 @@ class GGMModel : public BaseModel { constraint_dirty_(other.constraint_dirty_), theta_valid_(other.theta_valid_), theta_(other.theta_), - determinant_tilt_(other.determinant_tilt_), pcg_lambda_cache_(other.pcg_lambda_cache_) {} From c78ef0b6e52301d2af6cd251fb5611817394cbfd Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Tue, 19 May 2026 10:24:42 +0200 Subject: [PATCH 09/10] test: cover determinant-tilt path with SBC and NUTS-vs-MH at delta=1 MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three new slow tests gated behind BGMS_RUN_SLOW_TESTS: - GGM NUTS SBC at delta=1 (p=3, R=200), with a tilt-aware prior sampler using the shifted-shape Gamma proposal + (|K|/prod(K_ii))^delta acceptance from the spikeslab manuscript §sec:tilt-empirics. - GGM NUTS-vs-MH posterior moment concordance at delta=1 (p=4). - Mixed-MRF NUTS-vs-MH posterior concordance at delta=1 (Kyy block, no edge selection). Regenerate man/bgm.Rd via devtools::document() to pick up the delta parameter exposed by 8a9f7cb. --- man/bgm.Rd | 9 +++ tests/testthat/test-ggm-nuts.R | 62 +++++++++++++++++ tests/testthat/test-mixed-nuts.R | 39 +++++++++++ tests/testthat/test-sbc-ggm.R | 114 +++++++++++++++++++++++++++++++ 4 files changed, 224 insertions(+) diff --git a/man/bgm.Rd b/man/bgm.Rd index b476cbd7..10c3ed35 100644 --- a/man/bgm.Rd +++ b/man/bgm.Rd @@ -14,6 +14,7 @@ bgm( threshold_prior = beta_prime_prior(alpha = 0.5, beta = 0.5), means_prior = normal_prior(scale = 1), precision_scale_prior = gamma_prior(shape = 1, rate = 1), + delta = 0, edge_selection = TRUE, edge_prior = bernoulli_prior(0.5), na_action = c("listwise", "impute"), @@ -115,6 +116,14 @@ Only used for models with continuous variables (GGM and mixed MRF). Ignored for pure ordinal models. Default: \code{gamma_prior(shape = 1, rate = 1)}.} +\item{delta}{Non-negative numeric. Determinant-tilt exponent on the +continuous-block precision matrix \eqn{K} (GGM) or \eqn{K_{yy}} +(mixed MRF): multiplies the prior by \eqn{|K|^{\delta}}, pushing the +chain away from the positive-definite cone boundary. Both NUTS and +adaptive-Metropolis update paths apply the tilt. \code{delta = 0} +(default) recovers the untilted prior. Not allowed for pure ordinal +models (no precision matrix to tilt).} + \item{edge_selection}{Logical. Whether to perform Bayesian edge selection. If \code{FALSE}, the model estimates all edges. Default: \code{TRUE}.} diff --git a/tests/testthat/test-ggm-nuts.R b/tests/testthat/test-ggm-nuts.R index f0830bde..f6cc0463 100644 --- a/tests/testthat/test-ggm-nuts.R +++ b/tests/testthat/test-ggm-nuts.R @@ -374,3 +374,65 @@ test_that("NUTS diagnostics are clean for well-specified model", { expect_equal(fit$nuts_diag$summary$max_tree_depth_hits, 0) expect_gt(fit$nuts_diag$summary$min_ebfmi, 0.3) }) + + +# ---- 3.5 Tilt: NUTS vs MH posterior moments at delta = 1 ------------------- + +test_that("NUTS and MH posteriors agree under tilt (p=4, delta=1)", { + skip_unless_slow() + + dat = generate_ggm_data(p = 4, n = 200, seed = 36) + n_iter = 5000 + n_warmup = 2000 + + fit_mh = bgm( + x = dat$x, variable_type = "continuous", + update_method = "adaptive-metropolis", + edge_selection = FALSE, + delta = 1, + iter = n_iter, warmup = n_warmup, chains = 2, + seed = 102, display_progress = "none" + ) + + fit_nuts = bgm( + x = dat$x, variable_type = "continuous", + update_method = "nuts", + edge_selection = FALSE, + delta = 1, + iter = n_iter, warmup = n_warmup, chains = 2, + seed = 202, display_progress = "none" + ) + + mh_pairwise = do.call(rbind, fit_mh$raw_samples$pairwise) + nuts_pairwise = do.call(rbind, fit_nuts$raw_samples$pairwise) + mh_main = do.call(rbind, fit_mh$raw_samples$main) + nuts_main = do.call(rbind, fit_nuts$raw_samples$main) + + for(j in seq_len(ncol(mh_pairwise))) { + mh_mean = mean(mh_pairwise[, j]) + nuts_mean = mean(nuts_pairwise[, j]) + pooled_sd = sqrt((var(mh_pairwise[, j]) + var(nuts_pairwise[, j])) / 2) + expect_lt( + abs(mh_mean - nuts_mean) / pooled_sd, 2, + label = paste0("pairwise[", j, "] mean within 2 SDs (delta=1)") + ) + } + + for(j in seq_len(ncol(mh_main))) { + mh_mean = mean(mh_main[, j]) + nuts_mean = mean(nuts_main[, j]) + pooled_sd = sqrt((var(mh_main[, j]) + var(nuts_main[, j])) / 2) + expect_lt( + abs(mh_mean - nuts_mean) / pooled_sd, 2, + label = paste0("main[", j, "] mean within 2 SDs (delta=1)") + ) + } + + for(j in seq_len(ncol(mh_pairwise))) { + ratio = var(nuts_pairwise[, j]) / var(mh_pairwise[, j]) + expect_gt(ratio, 0.7, + label = paste0("pairwise[", j, "] var ratio > 0.7 (delta=1)")) + expect_lt(ratio, 1.4, + label = paste0("pairwise[", j, "] var ratio < 1.4 (delta=1)")) + } +}) diff --git a/tests/testthat/test-mixed-nuts.R b/tests/testthat/test-mixed-nuts.R index da616bec..ef00fa3a 100644 --- a/tests/testthat/test-mixed-nuts.R +++ b/tests/testthat/test-mixed-nuts.R @@ -584,3 +584,42 @@ test_that("M.2J: NUTS diagnostics are clean for mixed MRF", { expect_equal(fit$nuts_diag$summary$max_tree_depth_hits, 0) expect_gt(fit$nuts_diag$summary$min_ebfmi, 0.3) }) + + +# ---- Condition T: determinant tilt, conditional PL, no ES, grouped ---------- + +test_that("M.2T: NUTS vs MH agree under tilt (conditional PL, no ES, delta=1)", { + skip_unless_slow_mixed() + + dat = generate_mixed_test_data(seed = 2036) + vtype = c(rep("ordinal", 3), rep("continuous", 2)) + + fit_nuts = bgm( + dat, + variable_type = vtype, + iter = n_iter, warmup = n_warmup, chains = n_chains, + edge_selection = FALSE, update_method = "nuts", + pairwise_scale = pw_scale, main_alpha = main_a, main_beta = main_b, + delta = 1, + display_progress = "none", seed = 104 + ) + + fit_mh = bgm( + dat, + variable_type = vtype, + iter = n_iter, warmup = n_warmup, chains = n_chains, + edge_selection = FALSE, update_method = "adaptive-metropolis", + pairwise_scale = pw_scale, main_alpha = main_a, main_beta = main_b, + delta = 1, + display_progress = "none", seed = 204 + ) + + compare_pairwise_samples(fit_nuts, fit_mh, "M.2T") + compare_main_samples(fit_nuts, fit_mh, "M.2T") + + diff = abs(fit_nuts$posterior_mean_pairwise - fit_mh$posterior_mean_pairwise) + expect_lt(max(diff), 0.1, + label = "M.2T max abs diff in associations < 0.1 (delta=1)") + + expect_equal(fit_nuts$nuts_diag$summary$total_divergences, 0) +}) diff --git a/tests/testthat/test-sbc-ggm.R b/tests/testthat/test-sbc-ggm.R index a8790375..7ba720ee 100644 --- a/tests/testthat/test-sbc-ggm.R +++ b/tests/testthat/test-sbc-ggm.R @@ -386,3 +386,117 @@ test_that("SBC: GGM MH produces uniform diagonal ranks (p=3, edge selection)", { info = sprintf("SBC global chi-squared p=%.4f (edge selection, diag)", chisq_p) ) }) + + +# ---- Tilted prior sampler (delta > 0) ---------------------------------------- + +# Draw K from the determinant-tilted GGM prior: +# p(K) propto |K|^delta * slab(K_ij) * diag(K_ii) * 1{K in M+} +# using the rejection scheme from the spikeslab manuscript sec:tilt-empirics. +# Proposal shifts the diagonal Gamma shape up by delta on the K_ii scale; +# acceptance probability r = (|K| / prod(K_ii))^delta lies in [0, 1] on PD +# support by Hadamard's inequality. Off-diagonals are unchanged from the +# untilted Cauchy-slab proposal. +draw_prior_K_tilted = function(p, scale = 2.5, delta = 1, + max_tries = 100000) { + for(attempt in seq_len(max_tries)) { + K = matrix(0, p, p) + + # Off-diagonal: Cauchy slab, K_ij = -2 * Omega_ij as in untilted sampler. + for(i in seq_len(p - 1)) { + for(j in (i + 1):p) { + omega_ij = rcauchy(1, 0, scale) + K[i, j] = -2 * omega_ij + K[j, i] = K[i, j] + } + } + + # Diagonal: shifted-shape Gamma proposal. Untilted draws K_ii ~ + # 2 * Gamma(1, 1) = Gamma(1, 1/2) on K_ii; tilted proposal draws K_ii ~ + # 2 * Gamma(1 + delta, 1) = Gamma(1 + delta, 1/2). The factor K_ii^delta + # from the bound is absorbed into the proposal density. + for(i in seq_len(p)) { + K[i, i] = 2 * rgamma(1, shape = 1 + delta, rate = 1) + } + + ev = eigen(K, symmetric = TRUE, only.values = TRUE)$values + if(!all(ev > 1e-8)) next + + # Tilt acceptance: log r = delta * (sum log eigenvalues - sum log diag). + log_r = delta * (sum(log(ev)) - sum(log(diag(K)))) + if(log(runif(1)) < log_r) { + return(K) + } + } + stop("draw_prior_K_tilted: failed in ", max_tries, " attempts") +} + + +# ---- SBC test: NUTS under determinant tilt ---------------------------------- + +test_that("SBC: GGM NUTS produces uniform ranks under tilt (p=3, delta=1)", { + skip_unless_slow_sbc() + + p = 3 + n = 100 + R = 200 + L = 999 + scale = 2.5 + delta = 1 + + set.seed(2029) + + prior_draws = vector("list", R) + for(r in seq_len(R)) { + prior_draws[[r]] = draw_prior_K_tilted(p, scale, delta) + } + + n_off = p * (p - 1) / 2 + n_params = n_off + p + ranks = matrix(NA_real_, nrow = R, ncol = n_params) + + for(r in seq_len(R)) { + K_true = prior_draws[[r]] + Sigma = solve(K_true) + X = MASS::mvrnorm(n, mu = rep(0, p), Sigma = Sigma) + dat = as.data.frame(X) + colnames(dat) = paste0("V", seq_len(p)) + + fit = bgm(dat, + variable_type = "continuous", + iter = L, warmup = 1000, chains = 1, + edge_selection = FALSE, update_method = "nuts", + pairwise_scale = scale, delta = delta, + display_progress = "none", seed = 2029L + r + ) + + ranks[r, ] = compute_sbc_ranks(K_true, p, fit) + if(r == 1) colnames(ranks) = names(ranks[1, ]) + } + + n_fail_ks = 0 + for(j in seq_len(ncol(ranks))) { + u = ranks[, j] / (L + 1) + p_val = suppressWarnings(ks.test(u, "punif")$p.value) + if(p_val <= 0.01) n_fail_ks = n_fail_ks + 1 + } + + max_fail = max(1, ceiling(n_params * 0.01 * 2)) + expect_true(n_fail_ks <= max_fail, + info = sprintf( + "SBC KS (delta=1): %d/%d parameters failed (limit %d)", + n_fail_ks, n_params, max_fail + ) + ) + + all_ranks = as.vector(ranks) + bins = cut(all_ranks / (L + 1), + breaks = seq(0, 1, length.out = 21), + include.lowest = TRUE + ) + counts = tabulate(bins, nbins = 20) + chisq_p = chisq.test(counts)$p.value + expect_true(chisq_p > 0.001, + info = sprintf("SBC global chi-squared (delta=1) p=%.4f", chisq_p) + ) +}) From 243602d09e1bcf0f40cfd4971f727535579a732c Mon Sep 17 00:00:00 2001 From: Maarten Marsman Date: Tue, 19 May 2026 11:20:55 +0200 Subject: [PATCH 10/10] fix(pkgdown): update reference to renamed sample_ggm_prior Commit 987cc8b renamed sample_precision_prior to sample_ggm_prior but left the topic reference in _pkgdown.yml pointing at the old name, causing build_reference_index() to fail on the missing topic. --- _pkgdown.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/_pkgdown.yml b/_pkgdown.yml index 6ddff3aa..bae7e810 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -32,7 +32,7 @@ reference: - simulate.bgmCompare - predict.bgms - predict.bgmCompare - - sample_precision_prior + - sample_ggm_prior - title: Prior constructors desc: Specify priors for model parameters, precision scale, and edge inclusion.