Standardize code for updating proposal sd.

Simplify recoding for "Collapse" condition.

Author: MaartenMarsman

R/utility_functions.R

@@ -828,22 +828,13 @@ compare_reformat_data = function(x,
         zx = x[, node]
         zy = y[, node]
 
-        cntr = 0
+        cntr = -1
         for(value in unq_vls) {
           #Collapse categories for one group when not observed in the other.
-          if(!any(zx == value) | !any(zy == value)) {
-            check = FALSE
-            if(cntr > 0) {
-              cntr = cntr - 1
-            }
-          } else {
-            check = TRUE
-          }
-          x[zx == value, node] = cntr
-          y[zy == value, node] = cntr
-
-          if(!(check == FALSE & cntr == 0))
+          if(any(zx == value) & any(zy == value))
             cntr = cntr + 1
+          x[zx == value, node] = max(0, cntr)
+          y[zy == value, node] = max(0, cntr)
         }
       } else {
         z = x[, node]

src/gibbs_functions.cpp

@@ -284,14 +284,16 @@ void metropolis_thresholds_blumecapel(NumericMatrix thresholds,
     log_prob = std::exp(log_prob);
   }
 
-  proposal_sd_blumecapel(variable, 0) = proposal_sd_blumecapel(variable, 0) +
+  double update_proposal_sd = proposal_sd_blumecapel(variable, 0) +
     (log_prob - target_ar) * std::exp(-log(t) * phi);
-  if(proposal_sd_blumecapel(variable, 0) < epsilon_lo) {
-    proposal_sd_blumecapel(variable, 0) = epsilon_lo;
-  } else if (proposal_sd_blumecapel(variable, 0) > epsilon_hi) {
-    proposal_sd_blumecapel(variable, 0) = epsilon_hi;
+
+  if(std::isnan(update_proposal_sd) == true) {
+    update_proposal_sd = 1.0;
   }
 
+  update_proposal_sd = std::clamp(update_proposal_sd, epsilon_lo, epsilon_hi);
+  proposal_sd_blumecapel(variable, 0) = update_proposal_sd;
+
   //----------------------------------------------------------------------------
   //Adaptive Metropolis for the quadratic Blume-Capel parameter
   //----------------------------------------------------------------------------
@@ -363,13 +365,16 @@ void metropolis_thresholds_blumecapel(NumericMatrix thresholds,
     log_prob = std::exp(log_prob);
   }
 
-  proposal_sd_blumecapel(variable, 1) = proposal_sd_blumecapel(variable, 1) +
+  double update_proposal_sd = proposal_sd_blumecapel(variable, 1) +
     (log_prob - target_ar) * std::exp(-log(t) * phi);
-  if(proposal_sd_blumecapel(variable, 1) < epsilon_lo) {
-    proposal_sd_blumecapel(variable, 1) = epsilon_lo;
-  } else if (proposal_sd_blumecapel(variable, 1) > epsilon_hi) {
-    proposal_sd_blumecapel(variable, 1) = epsilon_hi;
+
+  if(std::isnan(update_proposal_sd) == true) {
+    update_proposal_sd = 1.0;
   }
+
+  update_proposal_sd = std::clamp(update_proposal_sd, epsilon_lo, epsilon_hi);
+  proposal_sd_blumecapel(variable, 1) = update_proposal_sd;
+
 }
 
 // ----------------------------------------------------------------------------|
@@ -556,13 +561,17 @@ void metropolis_interactions(NumericMatrix interactions,
         } else {
           log_prob = std::exp(log_prob);
         }
-        proposal_sd(variable1, variable2) = proposal_sd(variable1, variable2) +
+
+        double update_proposal_sd = proposal_sd(variable1, variable2) +
           (log_prob - target_ar) * std::exp(-log(t) * phi);
-        if(proposal_sd(variable1, variable2) < epsilon_lo) {
-          proposal_sd(variable1, variable2) = epsilon_lo;
-        } else if (proposal_sd(variable1, variable2) > epsilon_hi) {
-          proposal_sd(variable1, variable2) = epsilon_hi;
+
+        if(std::isnan(update_proposal_sd) == true) {
+          update_proposal_sd = 1.0;
         }
+
+        update_proposal_sd = std::clamp(update_proposal_sd, epsilon_lo, epsilon_hi);
+        proposal_sd(variable1, variable2) = update_proposal_sd;
+
       }
     }
   }
@@ -809,10 +818,10 @@ List gibbs_sampler(IntegerMatrix observations,
   IntegerMatrix index(no_interactions, 3);
 
   //Parameters of adaptive proposals -------------------------------------------
-  double phi =  0.750;
+  double phi =  0.75;
   double target_ar =  0.234;
-  double epsilon_lo = 1.0 / no_persons;
-  double epsilon_hi = 2.000;
+  double epsilon_lo = 1.0 / static_cast<double>(no_persons);
+  double epsilon_hi = 2.0;
 
   //The resizing based on ``save'' could probably be prettier ------------------
   int nrow = no_variables;

src/gibbs_functions_compare.cpp

@@ -2595,13 +2595,13 @@ List compare_gibbs_sampler(IntegerMatrix observations_gr1,
   std::fill(proposal_sd_blumecapel_gr2.begin(), proposal_sd_blumecapel_gr2.end(), 1.0);
 
   //Parameters for the Robbins-Monro approach for adaptive Metropolis ----------
-  double phi =        0.750;
-  double target_ar =  0.234;
-  double epsilon_lo = 1.0 / no_persons_gr1;
+  double phi = 0.75;
+  double target_ar = 0.234;
+  double epsilon_lo = 1.0 / static_cast<double>(no_persons_gr1);
   if(no_persons_gr1 > no_persons_gr2) {
-    epsilon_lo = 1.0 / no_persons_gr2;
+    epsilon_lo = 1.0 / static_cast<double>(no_persons_gr2);
   }
-  double epsilon_hi = 2.000;
+  double epsilon_hi = 2.0;
 
   //Randomized index for the pairwise updates ----------------------------------
   IntegerVector v = seq(0, no_interactions - 1);