minor code simplification

QingLicsaggie · Oct 3, 2017 · c010b58 · c010b58
1 parent c0ab66f
commit c010b58
Show file tree

Hide file tree

Showing 4 changed files with 23 additions and 59 deletions.
diff --git a/src/de.cpp b/src/de.cpp
@@ -133,11 +133,12 @@ mcmc::de_int(const arma::vec& initial_vals, arma::cube& draws_out, std::function
                     prop_kernel_val = BIG_NEG_VAL;
                 }
 
-                double comp_val = prop_kernel_val - target_vals(i);
-
                 //
+
+                double comp_val = prop_kernel_val - target_vals(i);
+                double z = arma::as_scalar(arma::randu(1,1));
 
-                if (comp_val > temperature_j * std::log(arma::as_scalar(arma::randu(1)))) {
+                if (comp_val > temperature_j * std::log(z)) {
                     X.row(i) = X_prop;
 
                     target_vals(i) = prop_kernel_val;

diff --git a/src/hmc.cpp b/src/hmc.cpp
@@ -122,8 +122,7 @@ mcmc::hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function
 
     //
 
-    int n_accept = 0;    
-    double comp_val, rand_val;
+    int n_accept = 0;
     arma::vec krand(n_vals);
 
     for (int jj = 0; jj < n_draws_keep + n_draws_burnin; jj++) {
@@ -151,9 +150,10 @@ mcmc::hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function
 
         //
 
-        comp_val = - prop_U - prop_K + prev_U + prev_K;
-
-        if (comp_val > 0.0) { // the '> exp(0)' case; works around taking exp of big values and receiving an error
+        double comp_val = std::min(0.0,- prop_U - prop_K + prev_U + prev_K);
+        double z = arma::as_scalar(arma::randu(1));
+
+        if (z < std::exp(comp_val)) {
             prev_draw = new_draw;
             prev_U = prop_U;
             prev_K = prop_K;
@@ -163,21 +163,8 @@ mcmc::hmc_int(const arma::vec& initial_vals, arma::mat& draws_out, std::function
                 n_accept++;
             }
         } else {
-            rand_val = arma::as_scalar(arma::randu(1));
-
-            if (rand_val < std::exp(comp_val)) {
-                prev_draw = new_draw;
-                prev_U = prop_U;
-                prev_K = prop_K;
-
-                if (jj >= n_draws_burnin) {
-                    draws_out.row(jj - n_draws_burnin) = new_draw.t();
-                    n_accept++;
-                }
-            } else {
-                if (jj >= n_draws_burnin) {
-                    draws_out.row(jj - n_draws_burnin) = prev_draw.t();
-                }
+            if (jj >= n_draws_burnin) {
+                draws_out.row(jj - n_draws_burnin) = prev_draw.t();
             }
         }
     }

diff --git a/src/mala.cpp b/src/mala.cpp
@@ -117,7 +117,6 @@ mcmc::mala_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
     //
 
     int n_accept = 0;    
-    double comp_val, rand_val;
     arma::vec krand(n_vals);
 
     for (int jj = 0; jj < n_draws_keep + n_draws_burnin; jj++) {
@@ -132,9 +131,10 @@ mcmc::mala_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
 
         //
 
-        comp_val = prop_LP - prev_LP + mala_prop_adjustment(new_draw, prev_draw, step_size, vals_bound, mala_mean_fn, target_data);
-
-        if (comp_val > 0.0) { // the '> exp(0)' case; works around taking exp of big values and receiving an error
+        double comp_val = std::min(0.0, prop_LP - prev_LP + mala_prop_adjustment(new_draw, prev_draw, step_size, vals_bound, mala_mean_fn, target_data));
+        double z = arma::as_scalar(arma::randu(1));
+
+        if (z < std::exp(comp_val)) {
             prev_draw = new_draw;
             prev_LP = prop_LP;
 
@@ -143,20 +143,8 @@ mcmc::mala_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
                 n_accept++;
             }
         } else {
-            rand_val = arma::as_scalar(arma::randu(1));
-
-            if (rand_val < std::exp(comp_val)) {
-                prev_draw = new_draw;
-                prev_LP = prop_LP;
-
-                if (jj >= n_draws_burnin) {
-                    draws_out.row(jj - n_draws_burnin) = new_draw.t();
-                    n_accept++;
-                }
-            } else {
-                if (jj >= n_draws_burnin) {
-                    draws_out.row(jj - n_draws_burnin) = prev_draw.t();
-                }
+            if (jj >= n_draws_burnin) {
+                draws_out.row(jj - n_draws_burnin) = prev_draw.t();
             }
         }
     }

diff --git a/src/rwmh.cpp b/src/rwmh.cpp
@@ -89,7 +89,6 @@ mcmc::rwmh_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
     //
 
     int n_accept = 0;    
-    double comp_val, rand_val;
     arma::vec krand(n_vals);
 
     for (int jj = 0; jj < n_draws_keep + n_draws_burnin; jj++) {
@@ -104,9 +103,10 @@ mcmc::rwmh_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
 
         //
 
-        comp_val = prop_LP - prev_LP;
-
-        if (comp_val > 0.0) { // the '> exp(0)' case; works around taking exp of big values and receiving an error
+        double comp_val = std::min(0.0,prop_LP - prev_LP);
+        double z = arma::as_scalar(arma::randu(1));
+
+        if (z < std::exp(comp_val)) {
             prev_draw = new_draw;
             prev_LP = prop_LP;
 
@@ -115,20 +115,8 @@ mcmc::rwmh_int(const arma::vec& initial_vals, arma::mat& draws_out, std::functio
                 n_accept++;
             }
         } else {
-            rand_val = arma::as_scalar(arma::randu(1));
-
-            if (rand_val < std::exp(comp_val)) {
-                prev_draw = new_draw;
-                prev_LP = prop_LP;
-
-                if (jj >= n_draws_burnin) {
-                    draws_out.row(jj - n_draws_burnin) = new_draw.t();
-                    n_accept++;
-                }
-            } else {
-                if (jj >= n_draws_burnin) {
-                    draws_out.row(jj - n_draws_burnin) = prev_draw.t();
-                }
+            if (jj >= n_draws_burnin) {
+                draws_out.row(jj - n_draws_burnin) = prev_draw.t();
             }
         }
     }