Created functions for computing marginal inclusion probabilities, stashed in archived code/gibbs_functions_moms.cpp

Author: MaartenMarsman

R/RcppExports.R

@@ -9,10 +9,6 @@ sample_bcomrf_gibbs <- function(no_states, no_variables, no_categories, interact
     .Call(`_bgms_sample_bcomrf_gibbs`, no_states, no_variables, no_categories, interactions, thresholds, variable_type, reference_category, iter)
 }
 
-optimize_log_pseudoposterior_interaction <- function(initial_value, pairwise_effects, main_effects, observations, num_categories, num_persons, variable1, variable2, proposed_state, current_state, residual_matrix, is_ordinal_variable, reference_category, interaction_scale) {
-    .Call(`_bgms_optimize_log_pseudoposterior_interaction`, initial_value, pairwise_effects, main_effects, observations, num_categories, num_persons, variable1, variable2, proposed_state, current_state, residual_matrix, is_ordinal_variable, reference_category, interaction_scale)
-}
-
 run_gibbs_sampler_for_bgm <- function(observations, num_categories, interaction_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, burnin, num_obs_categories, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, is_ordinal_variable, reference_category, save_main = FALSE, save_pairwise = FALSE, save_indicator = FALSE, display_progress = FALSE, edge_selection = TRUE, update_method = "adaptive-metropolis") {
     .Call(`_bgms_run_gibbs_sampler_for_bgm`, observations, num_categories, interaction_scale, edge_prior, inclusion_probability, beta_bernoulli_alpha, beta_bernoulli_beta, dirichlet_alpha, lambda, interaction_index_matrix, iter, burnin, num_obs_categories, sufficient_blume_capel, threshold_alpha, threshold_beta, na_impute, missing_index, is_ordinal_variable, reference_category, save_main, save_pairwise, save_indicator, display_progress, edge_selection, update_method)
 }

src/RcppExports.cpp

@@ -45,30 +45,6 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
-// optimize_log_pseudoposterior_interaction
-double optimize_log_pseudoposterior_interaction(const double initial_value, arma::mat& pairwise_effects, const arma::mat& main_effects, const arma::imat& observations, const arma::ivec& num_categories, const int num_persons, const int variable1, const int variable2, const double proposed_state, const double current_state, const arma::mat& residual_matrix, const arma::uvec& is_ordinal_variable, const arma::ivec& reference_category, const double interaction_scale);
-RcppExport SEXP _bgms_optimize_log_pseudoposterior_interaction(SEXP initial_valueSEXP, SEXP pairwise_effectsSEXP, SEXP main_effectsSEXP, SEXP observationsSEXP, SEXP num_categoriesSEXP, SEXP num_personsSEXP, SEXP variable1SEXP, SEXP variable2SEXP, SEXP proposed_stateSEXP, SEXP current_stateSEXP, SEXP residual_matrixSEXP, SEXP is_ordinal_variableSEXP, SEXP reference_categorySEXP, SEXP interaction_scaleSEXP) {
-BEGIN_RCPP
-    Rcpp::RObject rcpp_result_gen;
-    Rcpp::RNGScope rcpp_rngScope_gen;
-    Rcpp::traits::input_parameter< const double >::type initial_value(initial_valueSEXP);
-    Rcpp::traits::input_parameter< arma::mat& >::type pairwise_effects(pairwise_effectsSEXP);
-    Rcpp::traits::input_parameter< const arma::mat& >::type main_effects(main_effectsSEXP);
-    Rcpp::traits::input_parameter< const arma::imat& >::type observations(observationsSEXP);
-    Rcpp::traits::input_parameter< const arma::ivec& >::type num_categories(num_categoriesSEXP);
-    Rcpp::traits::input_parameter< const int >::type num_persons(num_personsSEXP);
-    Rcpp::traits::input_parameter< const int >::type variable1(variable1SEXP);
-    Rcpp::traits::input_parameter< const int >::type variable2(variable2SEXP);
-    Rcpp::traits::input_parameter< const double >::type proposed_state(proposed_stateSEXP);
-    Rcpp::traits::input_parameter< const double >::type current_state(current_stateSEXP);
-    Rcpp::traits::input_parameter< const arma::mat& >::type residual_matrix(residual_matrixSEXP);
-    Rcpp::traits::input_parameter< const arma::uvec& >::type is_ordinal_variable(is_ordinal_variableSEXP);
-    Rcpp::traits::input_parameter< const arma::ivec& >::type reference_category(reference_categorySEXP);
-    Rcpp::traits::input_parameter< const double >::type interaction_scale(interaction_scaleSEXP);
-    rcpp_result_gen = Rcpp::wrap(optimize_log_pseudoposterior_interaction(initial_value, pairwise_effects, main_effects, observations, num_categories, num_persons, variable1, variable2, proposed_state, current_state, residual_matrix, is_ordinal_variable, reference_category, interaction_scale));
-    return rcpp_result_gen;
-END_RCPP
-}
 // run_gibbs_sampler_for_bgm
 List run_gibbs_sampler_for_bgm(arma::imat& observations, const arma::ivec& num_categories, const double interaction_scale, const String& edge_prior, arma::mat& inclusion_probability, const double beta_bernoulli_alpha, const double beta_bernoulli_beta, const double dirichlet_alpha, const double lambda, const arma::imat& interaction_index_matrix, const int iter, const int burnin, arma::imat& num_obs_categories, arma::imat& sufficient_blume_capel, const double threshold_alpha, const double threshold_beta, const bool na_impute, const arma::imat& missing_index, const arma::uvec& is_ordinal_variable, const arma::ivec& reference_category, const bool save_main, const bool save_pairwise, const bool save_indicator, const bool display_progress, bool edge_selection, const std::string& update_method);
 RcppExport SEXP _bgms_run_gibbs_sampler_for_bgm(SEXP observationsSEXP, SEXP num_categoriesSEXP, SEXP interaction_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 burninSEXP, SEXP num_obs_categoriesSEXP, SEXP sufficient_blume_capelSEXP, SEXP threshold_alphaSEXP, SEXP threshold_betaSEXP, SEXP na_imputeSEXP, SEXP missing_indexSEXP, SEXP is_ordinal_variableSEXP, SEXP reference_categorySEXP, SEXP save_mainSEXP, SEXP save_pairwiseSEXP, SEXP save_indicatorSEXP, SEXP display_progressSEXP, SEXP edge_selectionSEXP, SEXP update_methodSEXP) {
@@ -167,7 +143,6 @@ END_RCPP
 static const R_CallMethodDef CallEntries[] = {
     {"_bgms_sample_omrf_gibbs", (DL_FUNC) &_bgms_sample_omrf_gibbs, 6},
     {"_bgms_sample_bcomrf_gibbs", (DL_FUNC) &_bgms_sample_bcomrf_gibbs, 8},
-    {"_bgms_optimize_log_pseudoposterior_interaction", (DL_FUNC) &_bgms_optimize_log_pseudoposterior_interaction, 14},
     {"_bgms_run_gibbs_sampler_for_bgm", (DL_FUNC) &_bgms_run_gibbs_sampler_for_bgm, 26},
     {"_bgms_compare_anova_gibbs_sampler", (DL_FUNC) &_bgms_compare_anova_gibbs_sampler, 34},
     {"_bgms_compute_Vn_mfm_sbm", (DL_FUNC) &_bgms_compute_Vn_mfm_sbm, 4},

src/gibbs_functions.cpp

@@ -9,7 +9,6 @@
 using namespace Rcpp;
 
 
-
 /**
  * Adapts the log step size using dual averaging during MCMC burn-in.
  *
@@ -1398,110 +1397,6 @@ double gradient_log_pseudoposterior_interaction_single (
 }
 
 
-double hessian_log_pseudoposterior_interaction_single (
-    int var1,
-    int var2,
-    const arma::mat& pairwise_effects,
-    const arma::mat& main_effects,
-    const arma::imat& observations,
-    const arma::ivec& num_categories,
-    const arma::uvec& is_ordinal_variable,
-    const arma::ivec& reference_category,
-    const double interaction_scale
-) {
-  const int num_persons = observations.n_rows;
-
-  // Extract observed score vectors for each variable
-  arma::vec x_var1 = arma::conv_to<arma::vec>::from (observations.col (var1));
-  arma::vec x_var2 = arma::conv_to<arma::vec>::from (observations.col (var2));
-
-  // First-order gradient from data
-  double hessian = 0.0;
-
-  // --- Contribution from var1
-  int num_categories_var1 = num_categories (var1);
-  arma::vec rest_scores_var1 = observations * pairwise_effects.col (var1);  // β_{var1,var1} = 0
-  arma::vec numerator_var1_E (num_persons, arma::fill::zeros);
-  arma::vec denominator_var1 (num_persons, arma::fill::zeros);
-  arma::vec numerator_var1_E2 (num_persons, arma::fill::zeros);
-  arma::vec bounds_var1 = arma::max (rest_scores_var1, arma::zeros<arma::vec> (num_persons)) * num_categories_var1;
-
-  if (is_ordinal_variable (var1)) {
-    denominator_var1 += arma::exp ( -bounds_var1 );
-    for (int category = 0; category < num_categories_var1; category++) {
-      arma::vec exponent = main_effects (var1, category) + (category + 1) * rest_scores_var1 - bounds_var1;
-      arma::vec weight = arma::exp (exponent);
-      denominator_var1 += weight;
-      numerator_var1_E += (category + 1) * x_var2 % weight;
-      numerator_var1_E2 += (category + 1) * (category + 1) * x_var2 % x_var2 % weight;
-    }
-  } else {
-    const int ref_cat = reference_category (var1);
-    for (int category = 0; category <= num_categories_var1; category++) {
-      int centered = category - ref_cat;
-      double lin_term = main_effects (var1, 0) * category;
-      double quad_term = main_effects (var1, 1) * centered * centered;
-      arma::vec exponent = lin_term + quad_term + category * rest_scores_var1 - bounds_var1;
-      arma::vec weight = arma::exp (exponent);
-      denominator_var1 += weight;
-      numerator_var1_E += category * x_var2 % weight;
-      numerator_var1_E2 += category * category * x_var2 % x_var2 % weight;
-    }
-  }
-  //- E((XiXj)^2)
-  hessian -= arma::accu (numerator_var1_E2 / denominator_var1);
-
-  //+E(XiXj)^2
-  arma::vec expectation = numerator_var1_E / denominator_var1;
-  hessian += arma::accu(arma::square(expectation));
-
-  // --- Contribution from var2
-  int num_categories_var2 = num_categories (var2);
-  arma::vec rest_scores_var2 = observations * pairwise_effects.col (var2);
-  arma::vec numerator_var2_E (num_persons, arma::fill::zeros);
-  arma::vec numerator_var2_E2 (num_persons, arma::fill::zeros);
-  arma::vec denominator_var2 (num_persons, arma::fill::zeros);
-  arma::vec bounds_var2 = arma::max (rest_scores_var2, arma::zeros<arma::vec> (num_persons)) * num_categories_var2;
-
-  if (is_ordinal_variable (var2)) {
-    denominator_var2 += arma::exp ( -bounds_var2 );
-    for (int category = 0; category < num_categories_var2; category++) {
-      arma::vec exponent = main_effects (var2, category) + (category + 1) * rest_scores_var2 - bounds_var2;
-      arma::vec weight = arma::exp (exponent);
-      denominator_var2 += weight;
-      numerator_var2_E += (category + 1) * x_var1 % weight;
-      numerator_var2_E2 += (category + 1) * (category + 1) * x_var1 % x_var1 % weight;
-    }
-  } else {
-    const int ref_cat = reference_category (var2);
-    for (int category = 0; category <= num_categories_var2; category++) {
-      int centered = category - ref_cat;
-      double lin_term = main_effects (var2, 0) * category;
-      double quad_term = main_effects (var2, 1) * centered * centered;
-      arma::vec exponent = lin_term + quad_term + category * rest_scores_var2 - bounds_var2;
-      arma::vec weight = arma::exp (exponent);
-      denominator_var2 += weight;
-      numerator_var2_E += category * x_var1 % weight;
-      numerator_var2_E2 += category * category * x_var1 % x_var1 % weight;
-    }
-  }
-
-  //- E((XiXj)^2)
-  hessian -= arma::accu (numerator_var2_E2 / denominator_var2);
-
-  //+E(XiXj)^2
-  expectation = numerator_var2_E / denominator_var2;
-  hessian += arma::accu(arma::square(expectation));
-
-
-  // --- Cauchy prior derivative
-  double beta = pairwise_effects (var1, var2) * pairwise_effects (var1, var2);
-  double s = interaction_scale * interaction_scale;
-  hessian += 2.0 * (beta - s) / ((beta + s) * (beta + s));
-
-  return hessian;
-}
-
 
 /**
  * Function: log_pseudoposterior_interactions
@@ -1542,7 +1437,7 @@ double log_pseudoposterior_interactions (
   arma::mat real_observations = arma::conv_to<arma::mat>::from (observations);
 
   // Leading term: trace(X * B * X^T)
-  double log_pseudo_likelihood = arma::trace (real_observations * pairwise_effects * real_observations.t ());
+  double log_pseudo_posterior = arma::trace (real_observations * pairwise_effects * real_observations.t ());
 
   for (int var = 0; var < num_variables; var++) {
     int num_categories_var = num_categories (var);
@@ -1572,105 +1467,25 @@ double log_pseudoposterior_interactions (
     }
 
     // Subtract log partition function and bounds adjustment
-    log_pseudo_likelihood -= arma::accu (arma::log (denominator));
-    log_pseudo_likelihood -= arma::accu (bounds);
+    log_pseudo_posterior -= arma::accu (arma::log (denominator));
+    log_pseudo_posterior -= arma::accu (bounds);
   }
 
   // Add Cauchy prior terms for included pairwise effects
   for (int var1 = 0; var1 < num_variables - 1; var1++) {
     for (int var2 = var1 + 1; var2 < num_variables; var2++) {
       if (inclusion_indicator (var1, var2) == 1) {
-        log_pseudo_likelihood += R::dcauchy (pairwise_effects (var1, var2), 0.0, interaction_scale, true);
+        log_pseudo_posterior += R::dcauchy (pairwise_effects (var1, var2), 0.0, interaction_scale, true);
       }
     }
   }
 
-  return log_pseudo_likelihood;
+  return log_pseudo_posterior;
 }
 
 
-/**
- *
- *
- *
- */
-//[[Rcpp::export]]
-double optimize_log_pseudoposterior_interaction (
-    const double initial_value,
-    arma::mat& pairwise_effects,
-    const arma::mat& main_effects,
-    const arma::imat& inclusion_indicator,
-    const arma::imat& observations,
-    const arma::ivec& num_categories,
-    const int num_persons,
-    const int variable1,
-    const int variable2,
-    const double proposed_state,
-    const double current_state,
-    const arma::mat& residual_matrix,
-    const arma::uvec& is_ordinal_variable,
-    const arma::ivec& reference_category,
-    const double interaction_scale
-) {
-
-  double x = initial_value;
 
-  const int    max_steps = 10;
-  const double tolerance = 1e-6;//sqrt (std::numeric_limits<double>::epsilon ());
 
-  const double x0 = pairwise_effects(variable1, variable2);
-  double hessian_at_x;
-  // find mode
-  for (int t = 0; t < max_steps; t++) {
-
-    // TODO: need to assign x to pairwise_effects[variable1, variable2]
-    pairwise_effects(variable1, variable2) = x;
-    pairwise_effects(variable2, variable1) = x;
-    Rcpp::Rcout << "t: " << t << " x: " << x << std::endl;
-    double gradient_at_x = gradient_log_pseudoposterior_interaction_single (
-      variable1, variable2, pairwise_effects, main_effects, observations,
-      num_categories, is_ordinal_variable, reference_category, interaction_scale
-    );
-
-    Rcpp::Rcout << "hessian_at_x" << std::endl;
-    hessian_at_x = hessian_log_pseudoposterior_interaction_single (
-      variable1, variable2, pairwise_effects, main_effects, observations,
-      num_categories, is_ordinal_variable, reference_category, interaction_scale
-    );
-
-    // double x_new = x - gradient_at_x / hessian_at_x;
-    double x_new = x - gradient_at_x / hessian_at_x;
-
-    if (std::abs(x_new - x) < tolerance) {
-      x = x_new;
-      break;
-    }
-    x = x_new;
-
-  }
-
-  pairwise_effects(variable1, variable2) = x;
-  pairwise_effects(variable2, variable1) = x;
-
-  const double fx = log_pseudoposterior_interactions(
-    pairwise_effects,
-    main_effects,
-    observations,
-    num_categories,
-    inclusion_indicator,
-    is_ordinal_variable,
-    reference_category,
-    interaction_scale
-  );
-
-  pairwise_effects(variable1, variable2) = x0;
-  pairwise_effects(variable2, variable1) = x0;
-
-  // @maarten not sure if you need or want both?
-  const double log_integral = fx + (log(2 * M_PI) - log(-hessian_at_x)) / 2;
-  return x;
-
-}
 
 
 
@@ -2797,6 +2612,7 @@ void update_indicator_interaction_pair_with_fisher_mala (
 }
 
 
+
 /**
  * Performs a single iteration of the Gibbs sampler for graphical model parameters.
  *
@@ -2894,7 +2710,8 @@ void gibbs_update_step_for_graphical_model_parameters (
     arma::mat& sqrt_inv_fisher_pairwise,
     const std::string& update_method,
     arma::vec& cached_interaction_gradient,
-    bool& gradient_valid
+    bool& gradient_valid,
+    arma::vec& posterior_prob
 ) {
   // --- Robbins-Monro weight for adaptive Metropolis updates
   const double exp_neg_log_t_rm_adaptation_rate =
@@ -3088,11 +2905,11 @@ List run_gibbs_sampler_for_bgm (
   const int num_main = count_num_main_effects(num_categories, is_ordinal_variable);
   arma::mat* main_effect_samples = nullptr;
   arma::mat* pairwise_effect_samples = nullptr;
-  arma::imat* indicator_samples = nullptr;
+  arma::mat* indicator_samples = nullptr;
 
   if (save_main) main_effect_samples = new arma::mat(iter, num_main);
   if (save_pairwise) pairwise_effect_samples = new arma::mat(iter, num_pairwise);
-  if (save_indicator) indicator_samples = new arma::imat(iter, num_pairwise);
+  if (save_indicator) indicator_samples = new arma::mat(iter, num_pairwise);
 
   // Initialize proposal SDs and MALA tracking
   arma::mat proposal_sd_main(num_main, 2, arma::fill::ones);
@@ -3175,7 +2992,7 @@ List run_gibbs_sampler_for_bgm (
   }
   arma::vec cached_interaction_gradient;  // will hold a cached gradient vector
   bool gradient_valid = false;            // indicates whether the cache is valid
-
+  arma::vec posterior_prob(num_pairwise);
 
   // --- Set up total number of iterations (burn-in + sampling)
   bool enable_edge_selection = edge_selection;
@@ -3225,7 +3042,7 @@ List run_gibbs_sampler_for_bgm (
         dual_averaging_main, total_burnin, initial_step_size_main,
         sqrt_inv_fisher_main, step_size_pairwise, dual_averaging_pairwise,
         initial_step_size_pairwise, sqrt_inv_fisher_pairwise, update_method,
-        cached_interaction_gradient, gradient_valid
+        cached_interaction_gradient, gradient_valid, posterior_prob
     );
 
     // --- Update edge probabilities under the prior (if edge selection is active)
@@ -3295,11 +3112,11 @@ List run_gibbs_sampler_for_bgm (
       }
 
       if (save_indicator) {
-        arma::ivec vectorized_indicator(num_pairwise);
-        for (int i = 0; i < num_pairwise; i++) {
-          vectorized_indicator(i) = inclusion_indicator(interaction_index_matrix(i, 1), interaction_index_matrix(i, 2));
-        }
-        indicator_samples->row(sample_index) = vectorized_indicator.t();
+        //arma::ivec vectorized_indicator(num_pairwise);
+        //for (int i = 0; i < num_pairwise; i++) {
+        //  vectorized_indicator(i) = inclusion_indicator(interaction_index_matrix(i, 1), interaction_index_matrix(i, 2));
+        //}
+        indicator_samples->row(sample_index) = posterior_prob.t();//vectorized_indicator.t();
       }
 
       if (edge_prior == "Stochastic-Block") {