Skip to content

two sets of shape hyperparameters for the Beta hyperprior in the SBM #68

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
MaartenMarsman merged 6 commits into from
Nov 12, 2025
Merged

two sets of shape hyperparameters for the Beta hyperprior in the SBM #68

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
c8af72d
two sets of shape hyperparameters for the Beta hyperprior in the SBM
sekulovskin Oct 14, 2025
727f8cd
Merge branch 'Bayesian-Graphical-Modelling-Lab:main' into SBM
sekulovskin Oct 27, 2025
bf3ec25
Update documentation for the extra shape parameters
sekulovskin Oct 27, 2025
8182afd
Merge branch 'SBM' of https://github.com/sekulovskin/bgms into SBM
sekulovskin Oct 27, 2025
a462165
update documentation
sekulovskin Oct 27, 2025
dc4de41
adress review comments
Nov 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions R/RcppExports.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,8 @@ run_bgmCompare_parallel <- function(observations, num_groups, counts_per_categor
.Call(`_bgms_run_bgmCompare_parallel`, observations, num_groups, counts_per_category, blume_capel_stats, pairwise_stats, num_categories, main_alpha, main_beta, pairwise_scale, difference_scale, difference_selection_alpha, difference_selection_beta, difference_prior, iter, warmup, na_impute, missing_data_indices, is_ordinal_variable, baseline_category, difference_selection, main_effect_indices, pairwise_effect_indices, target_accept, nuts_max_depth, learn_mass_matrix, projection, group_membership, group_indices, interaction_index_matrix, inclusion_probability, num_chains, nThreads, seed, update_method, hmc_num_leapfrogs, progress_type)
}

run_bgm_parallel <- function(observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type) {
.Call(`_bgms_run_bgm_parallel`, observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type)
run_bgm_parallel <- function(observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type) {
.Call(`_bgms_run_bgm_parallel`, observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type)
}

get_explog_switch <- function() {
Expand Down
30 changes: 28 additions & 2 deletions R/bgm.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -210,6 +210,13 @@
#' priors. Must be positive. Defaults: \code{beta_bernoulli_alpha = 1} and
#' \code{beta_bernoulli_beta = 1}.
#'
#' @param beta_bernoulli_alpha_between,beta_bernoulli_beta_between Double.
#' Second set of shape parameters for the Stochastic-Block prior. Default: \code{NULL}.
#' If the users supply values for these parameters, the model will use different
#' between-cluster and within-cluster edge inclusion probabilities. In this case
#' the beta_bernoulli_alpha and beta_bernoulli_beta will be used as the within-cluster
#' parameters.
#'
#' @param dirichlet_alpha Double. Concentration parameter of the Dirichlet
#' prior on block assignments (used with the Stochastic Block model).
#' Default: \code{1}.
Expand Down Expand Up @@ -239,7 +246,7 @@
#' @param target_accept Numeric between 0 and 1. Target acceptance rate for
#' the sampler. Defaults are set automatically if not supplied:
#' \code{0.44} for adaptive Metropolis, \code{0.65} for HMC,
#' and \code{0.80} for NUTS.
#' and \code{0.60} for NUTS.
#'
#' @param hmc_num_leapfrogs Integer. Number of leapfrog steps for Hamiltonian
#' Monte Carlo. Must be positive. Default: \code{100}.
Expand Down Expand Up @@ -359,6 +366,8 @@ bgm = function(
inclusion_probability = 0.5,
beta_bernoulli_alpha = 1,
beta_bernoulli_beta = 1,
beta_bernoulli_alpha_between = NULL,
beta_bernoulli_beta_between = NULL,
dirichlet_alpha = 1,
lambda = 1,
na_action = c("listwise", "impute"),
Expand Down Expand Up @@ -418,7 +427,7 @@ bgm = function(
} else if(update_method == "hamiltonian-mc") {
target_accept = 0.65
} else if(update_method == "nuts") {
target_accept = 0.80
target_accept = 0.60
}
}

Expand All @@ -444,9 +453,21 @@ bgm = function(
inclusion_probability = inclusion_probability,
beta_bernoulli_alpha = beta_bernoulli_alpha,
beta_bernoulli_beta = beta_bernoulli_beta,
beta_bernoulli_alpha_between = beta_bernoulli_alpha_between,
beta_bernoulli_beta_between = beta_bernoulli_beta_between,
dirichlet_alpha = dirichlet_alpha,
lambda = lambda)

# check hyperparameters input
# If user left them NULL, pass -1 to C++ (means: ignore between prior)
if (is.null(beta_bernoulli_alpha_between) && is.null(beta_bernoulli_beta_between)) {
beta_bernoulli_alpha_between <- -1.0
beta_bernoulli_beta_between <- -1.0
} else if (is.null(beta_bernoulli_alpha_between) || is.null(beta_bernoulli_beta_between)) {
stop("If you wish to specify different between and within cluster probabilites,
provide both beta_bernoulli_alpha_between and beta_bernoulli_beta_between,
otherwise leave both NULL.")
}
Copy link
Collaborator

@MaartenMarsman MaartenMarsman Nov 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If I read this correctly, this will stop the analysis if the between parameters are left at their defaults. This means users cannot use bgms, except when they set other values for these defaults, irrespective if they use an SBM prior or not.

Copy link
Contributor Author

@sekulovskin sekulovskin Nov 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No, this is in case the user specifies only one of the between parameters, but not both. Otherwise running the code below did work. But it doesn't matter now, since as you requested, both sets of hyperparameters are always provided.

fit <- bgm(Wenchuan[c(1:100), c(1:5)], edge_prior = "Stochastic-Block",
beta_bernoulli_alpha = 1,
beta_bernoulli_beta = 1,
beta_bernoulli_alpha_between = NULL,
beta_bernoulli_beta_between = NULL,
update_method = "adaptive-metropolis",
iter = 1000)

# ----------------------------------------------------------------------------
# The vector variable_type is now coded as boolean.
# Ordinal (variable_bool == TRUE) or Blume-Capel (variable_bool == FALSE)
Expand Down Expand Up @@ -572,6 +593,8 @@ bgm = function(
inclusion_probability = inclusion_probability,
beta_bernoulli_alpha = beta_bernoulli_alpha,
beta_bernoulli_beta = beta_bernoulli_beta,
beta_bernoulli_alpha_between = beta_bernoulli_alpha_between,
beta_bernoulli_beta_between = beta_bernoulli_beta_between,
dirichlet_alpha = dirichlet_alpha, lambda = lambda,
interaction_index_matrix = interaction_index_matrix, iter = iter,
warmup = warmup, counts_per_category = counts_per_category,
Expand Down Expand Up @@ -603,6 +626,7 @@ bgm = function(
na_action = na_action, na_impute = na_impute,
edge_selection = edge_selection, edge_prior = edge_prior, inclusion_probability = inclusion_probability,
beta_bernoulli_alpha = beta_bernoulli_alpha, beta_bernoulli_beta = beta_bernoulli_beta,
beta_bernoulli_alpha_between = beta_bernoulli_alpha_between, beta_bernoulli_beta_between = beta_bernoulli_beta_between,
dirichlet_alpha = dirichlet_alpha, lambda = lambda,
variable_type = variable_type,
update_method = update_method,
Expand Down Expand Up @@ -634,6 +658,8 @@ bgm = function(
edge_selection = edge_selection, edge_prior = edge_prior, inclusion_probability = inclusion_probability,
beta_bernoulli_alpha = beta_bernoulli_alpha,
beta_bernoulli_beta = beta_bernoulli_beta,
beta_bernoulli_alpha_between = beta_bernoulli_alpha_between,
beta_bernoulli_beta_between = beta_bernoulli_beta_between,
dirichlet_alpha = dirichlet_alpha, lambda = lambda,
variable_type = variable_type,
update_method = update_method,
Expand Down
43 changes: 33 additions & 10 deletions R/function_input_utils.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,8 @@ check_model = function(x,
inclusion_probability = 0.5,
beta_bernoulli_alpha = 1,
beta_bernoulli_beta = 1,
beta_bernoulli_alpha_between = NULL,
beta_bernoulli_beta_between = NULL,
dirichlet_alpha = dirichlet_alpha,
lambda = lambda) {

Expand Down Expand Up @@ -204,18 +206,39 @@ check_model = function(x,
is.null(beta_bernoulli_alpha) || is.null(beta_bernoulli_beta))
stop("Values for both scale parameters of the beta distribution need to be specified.")
}

if(edge_prior == "Stochastic-Block") {
theta = matrix(0.5, nrow = ncol(x), ncol = ncol(x))
if(beta_bernoulli_alpha <= 0 || beta_bernoulli_beta <= 0 || dirichlet_alpha <= 0 || lambda <= 0)
stop("The scale parameters of the beta and Dirichlet distribution need to be positive.")
if(!is.finite(beta_bernoulli_alpha) || !is.finite(beta_bernoulli_beta) || !is.finite(dirichlet_alpha) || !is.finite(lambda))
stop("The scale parameters of the beta distribution, the concentration parameter of the Dirichlet distribution, and the rate parameter of the Poisson distribution need to be finite.")
if(is.na(beta_bernoulli_alpha) || is.na(beta_bernoulli_beta) ||
is.null(beta_bernoulli_alpha) || is.null(beta_bernoulli_beta) ||
is.null(dirichlet_alpha) || is.null(dirichlet_alpha) || is.null(lambda) || is.null(lambda))
stop("Values for both scale parameters of the beta distribution, the concentration parameter of the Dirichlet distribution, and the rate parameter of the Poisson distribution need to be specified.")
theta = matrix(0.5, nrow = ncol(x), ncol = ncol(x))

## resolve optional "between" pair: both-or-neither; if neither, send -1 to the cpp functions
if (is.null(beta_bernoulli_alpha_between) && is.null(beta_bernoulli_beta_between)) {
.have_between <- FALSE
beta_bernoulli_alpha_between <- -1.0
beta_bernoulli_beta_between <- -1.0
} else if (is.null(beta_bernoulli_alpha_between) || is.null(beta_bernoulli_beta_between)) {
stop("If you wish to specify different between and within cluster probabilites,
provide both beta_bernoulli_alpha_between and beta_bernoulli_beta_between,
otherwise leave both NULL.")
} else {
.have_between <- TRUE
}
} else {

if (beta_bernoulli_alpha <= 0 || beta_bernoulli_beta <= 0 ||
(.have_between && (beta_bernoulli_alpha_between <= 0 || beta_bernoulli_beta_between <= 0)) ||
dirichlet_alpha <= 0 || lambda <= 0)
stop("The scale parameters of the beta and Dirichlet distribution need to be positive.")

if (!is.finite(beta_bernoulli_alpha) || !is.finite(beta_bernoulli_beta) ||
(.have_between && (!is.finite(beta_bernoulli_alpha_between) || !is.finite(beta_bernoulli_beta_between))) ||
!is.finite(dirichlet_alpha) || !is.finite(lambda))
stop("The scale parameters of the beta distribution, the concentration parameter of the Dirichlet distribution, and the rate parameter of the Poisson distribution need to be finite.")

if (is.na(beta_bernoulli_alpha) || is.na(beta_bernoulli_beta) ||
is.null(beta_bernoulli_alpha) || is.null(beta_bernoulli_beta) ||
is.null(dirichlet_alpha) || is.null(lambda))
stop("Values for both scale parameters of the beta distribution, the concentration parameter of the Dirichlet distribution, and the rate parameter of the Poisson distribution need to be specified.")
}
}else {
theta = matrix(0.5, nrow = 1, ncol = 1)
edge_prior = "Not Applicable"
}
Expand Down
5 changes: 4 additions & 1 deletion R/output_utils.R
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,8 @@ prepare_output_bgm = function(
out, x, num_categories, iter, data_columnnames, is_ordinal_variable,
warmup, pairwise_scale, main_alpha, main_beta,
na_action, na_impute, edge_selection, edge_prior, inclusion_probability,
beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda,
beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between,
beta_bernoulli_beta_between,dirichlet_alpha, lambda,
variable_type, update_method, target_accept, hmc_num_leapfrogs,
nuts_max_depth, learn_mass_matrix, num_chains
) {
Expand All @@ -22,6 +23,8 @@ prepare_output_bgm = function(
inclusion_probability = inclusion_probability,
beta_bernoulli_alpha = beta_bernoulli_alpha,
beta_bernoulli_beta = beta_bernoulli_beta,
beta_bernoulli_alpha_between = beta_bernoulli_alpha_between,
beta_bernoulli_beta_between = beta_bernoulli_beta_between,
dirichlet_alpha = dirichlet_alpha,
lambda = lambda,
na_action = na_action,
Expand Down
11 changes: 10 additions & 1 deletion man/bgm.Rd
View file
Open in desktop

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

10 changes: 6 additions & 4 deletions src/RcppExports.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,8 +58,8 @@ BEGIN_RCPP
END_RCPP
}
// run_bgm_parallel
Rcpp::List run_bgm_parallel(const arma::imat& observations, const arma::ivec& num_categories, double pairwise_scale, const std::string& edge_prior, const arma::mat& inclusion_probability, double beta_bernoulli_alpha, double beta_bernoulli_beta, double dirichlet_alpha, double lambda, const arma::imat& interaction_index_matrix, int iter, int warmup, const arma::imat& counts_per_category, const arma::imat& blume_capel_stats, double main_alpha, double main_beta, bool na_impute, const arma::imat& missing_index, const arma::uvec& is_ordinal_variable, const arma::ivec& baseline_category, bool edge_selection, const std::string& update_method, const arma::imat& pairwise_effect_indices, double target_accept, const arma::imat& pairwise_stats, int hmc_num_leapfrogs, int nuts_max_depth, bool learn_mass_matrix, int num_chains, int nThreads, int seed, int progress_type);
RcppExport SEXP _bgms_run_bgm_parallel(SEXP observationsSEXP, SEXP num_categoriesSEXP, SEXP pairwise_scaleSEXP, SEXP edge_priorSEXP, SEXP inclusion_probabilitySEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP interaction_index_matrixSEXP, SEXP iterSEXP, SEXP warmupSEXP, SEXP counts_per_categorySEXP, SEXP blume_capel_statsSEXP, SEXP main_alphaSEXP, SEXP main_betaSEXP, SEXP na_imputeSEXP, SEXP missing_indexSEXP, SEXP is_ordinal_variableSEXP, SEXP baseline_categorySEXP, SEXP edge_selectionSEXP, SEXP update_methodSEXP, SEXP pairwise_effect_indicesSEXP, SEXP target_acceptSEXP, SEXP pairwise_statsSEXP, SEXP hmc_num_leapfrogsSEXP, SEXP nuts_max_depthSEXP, SEXP learn_mass_matrixSEXP, SEXP num_chainsSEXP, SEXP nThreadsSEXP, SEXP seedSEXP, SEXP progress_typeSEXP) {
Rcpp::List run_bgm_parallel(const arma::imat& observations, const arma::ivec& num_categories, double pairwise_scale, const std::string& edge_prior, const arma::mat& inclusion_probability, double beta_bernoulli_alpha, double beta_bernoulli_beta, double beta_bernoulli_alpha_between, double beta_bernoulli_beta_between, double dirichlet_alpha, double lambda, const arma::imat& interaction_index_matrix, int iter, int warmup, const arma::imat& counts_per_category, const arma::imat& blume_capel_stats, double main_alpha, double main_beta, bool na_impute, const arma::imat& missing_index, const arma::uvec& is_ordinal_variable, const arma::ivec& baseline_category, bool edge_selection, const std::string& update_method, const arma::imat& pairwise_effect_indices, double target_accept, const arma::imat& pairwise_stats, int hmc_num_leapfrogs, int nuts_max_depth, bool learn_mass_matrix, int num_chains, int nThreads, int seed, int progress_type);
RcppExport SEXP _bgms_run_bgm_parallel(SEXP observationsSEXP, SEXP num_categoriesSEXP, SEXP pairwise_scaleSEXP, SEXP edge_priorSEXP, SEXP inclusion_probabilitySEXP, SEXP beta_bernoulli_alphaSEXP, SEXP beta_bernoulli_betaSEXP, SEXP beta_bernoulli_alpha_betweenSEXP, SEXP beta_bernoulli_beta_betweenSEXP, SEXP dirichlet_alphaSEXP, SEXP lambdaSEXP, SEXP interaction_index_matrixSEXP, SEXP iterSEXP, SEXP warmupSEXP, SEXP counts_per_categorySEXP, SEXP blume_capel_statsSEXP, SEXP main_alphaSEXP, SEXP main_betaSEXP, SEXP na_imputeSEXP, SEXP missing_indexSEXP, SEXP is_ordinal_variableSEXP, SEXP baseline_categorySEXP, SEXP edge_selectionSEXP, SEXP update_methodSEXP, SEXP pairwise_effect_indicesSEXP, SEXP target_acceptSEXP, SEXP pairwise_statsSEXP, SEXP hmc_num_leapfrogsSEXP, SEXP nuts_max_depthSEXP, SEXP learn_mass_matrixSEXP, SEXP num_chainsSEXP, SEXP nThreadsSEXP, SEXP seedSEXP, SEXP progress_typeSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Expand All @@ -70,6 +70,8 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< const arma::mat& >::type inclusion_probability(inclusion_probabilitySEXP);
Rcpp::traits::input_parameter< double >::type beta_bernoulli_alpha(beta_bernoulli_alphaSEXP);
Rcpp::traits::input_parameter< double >::type beta_bernoulli_beta(beta_bernoulli_betaSEXP);
Rcpp::traits::input_parameter< double >::type beta_bernoulli_alpha_between(beta_bernoulli_alpha_betweenSEXP);
Rcpp::traits::input_parameter< double >::type beta_bernoulli_beta_between(beta_bernoulli_beta_betweenSEXP);
Rcpp::traits::input_parameter< double >::type dirichlet_alpha(dirichlet_alphaSEXP);
Rcpp::traits::input_parameter< double >::type lambda(lambdaSEXP);
Rcpp::traits::input_parameter< const arma::imat& >::type interaction_index_matrix(interaction_index_matrixSEXP);
Expand All @@ -95,7 +97,7 @@ BEGIN_RCPP
Rcpp::traits::input_parameter< int >::type nThreads(nThreadsSEXP);
Rcpp::traits::input_parameter< int >::type seed(seedSEXP);
Rcpp::traits::input_parameter< int >::type progress_type(progress_typeSEXP);
rcpp_result_gen = Rcpp::wrap(run_bgm_parallel(observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type));
rcpp_result_gen = Rcpp::wrap(run_bgm_parallel(observations, num_categories, pairwise_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, beta_bernoulli_alpha_between, beta_bernoulli_beta_between, dirichlet_alpha, lambda, interaction_index_matrix, iter, warmup, counts_per_category, blume_capel_stats, main_alpha, main_beta, na_impute, missing_index, is_ordinal_variable, baseline_category, edge_selection, update_method, pairwise_effect_indices, target_accept, pairwise_stats, hmc_num_leapfrogs, nuts_max_depth, learn_mass_matrix, num_chains, nThreads, seed, progress_type));
return rcpp_result_gen;
END_RCPP
}
Expand Down Expand Up @@ -182,7 +184,7 @@ END_RCPP

static const R_CallMethodDef CallEntries[] = {
{"_bgms_run_bgmCompare_parallel", (DL_FUNC) &_bgms_run_bgmCompare_parallel, 36},
{"_bgms_run_bgm_parallel", (DL_FUNC) &_bgms_run_bgm_parallel, 32},
{"_bgms_run_bgm_parallel", (DL_FUNC) &_bgms_run_bgm_parallel, 34},
{"_bgms_get_explog_switch", (DL_FUNC) &_bgms_get_explog_switch, 0},
{"_bgms_rcpp_ieee754_exp", (DL_FUNC) &_bgms_rcpp_ieee754_exp, 1},
{"_bgms_rcpp_ieee754_log", (DL_FUNC) &_bgms_rcpp_ieee754_log, 1},
Expand Down
Loading