Merge pull request #1 from wellington36/fix_mcmc

Fix Bounding pairs approach
wellington36 · Oct 15, 2024 · 0b4e64c · 0b4e64c
2 parents aa26b9e + ed393e1
commit 0b4e64c
Show file tree

Hide file tree

Showing 3 changed files with 10 additions and 34 deletions.
diff --git a/stan/ErrorBounding.stan b/stan/ErrorBounding.stan
@@ -15,6 +15,6 @@ array[] real log_Z_COMP_EB(real log_mu, real nu, real eps, int M) {
     reject("nu must be finite");
   }
   return(
-    infiniteErrorBoundingPairs({log_mu, nu}, eps, M, log(0), 0)
+    infiniteErrorBoundingPairs({log_mu, nu}, eps, M, 0)
   );
 }
diff --git a/stan/MCMC_COMP_2_ErrorBounding.stan b/stan/MCMC_COMP_2_ErrorBounding.stan
@@ -14,16 +14,13 @@ data{
   real<lower=0> eps;
   int<lower=0> M;
 }
-transformed data{
-  real logL = log(0);
-}
 parameters{
   real<lower=0> mu;
   real<lower=0> nu;
 }
 transformed parameters{
   real log_mu = log(mu);
-  array[2] real log_norm_const = infiniteErrorBoundingPairs({log_mu, nu}, eps, M, logL, 0);
+  array[2] real log_norm_const = infiniteErrorBoundingPairs({log_mu, nu}, eps, M, 0);
 }
 model{
   mu ~ gamma(s_mu, r_mu);

diff --git a/stan/infiniteErrorBoundingPairs.stan b/stan/infiniteErrorBoundingPairs.stan
@@ -1,29 +1,9 @@
-// Internal use
-int adaptiveConvergenceCheck(real oldT, real newT, real lepsilon, real log1mL) {
-  // if L = 0, there is a simpler check
-  if (log1mL == 0) return oldT - newT < log1p_exp(newT - lepsilon);
-
-  real logZ = newT + oldT - log_diff_exp(oldT, newT);
-  real ls = newT - log1mL;
-
-  if (logZ > ls){
-    return log_diff_exp(logZ, ls) >= lepsilon;
-  }else{
-    return log_diff_exp(ls, logZ) >= lepsilon;
-  }
-}
-
-// Adaptive inifinite sum algorithm
-// Requires definition of logFunction with two arguments:
-// int k and real[] parameters
-array[] real infiniteErrorBoundingPairs(array[] real p, real epsilon, int maxIter, real logL, int n0) {
+array[] real infiniteErrorBoundingPairs(array[] real p, real epsilon, int maxIter, int n0) {
   vector[maxIter + 1] storeVal;
-  real leps = log(epsilon) + log2();
+  real leps = log(epsilon);
   int n = 2;
   int n0_ = n0;
-  real temp;
-  real log1mL = log_diff_exp(0, logL);
-
+
   // Setting up first iterations
   storeVal[1] = logFunction(n0_, p);
   n0_ += 1;
@@ -36,15 +16,14 @@ array[] real infiniteErrorBoundingPairs(array[] real p, real epsilon, int maxIte
     storeVal[n] = logFunction(n0_, p);
     if (n >= maxIter) return({log_sum_exp(storeVal[:n]), 1. * n}); // Return if maxIter is reached
   }
-  // Start testing convergence after the maximum
-  while ( adaptiveConvergenceCheck(storeVal[n - 1], storeVal[n], leps, log1mL) )  {
+
+  while (storeVal[n] - log(-expm1(storeVal[n] - storeVal[n - 1])) >= leps + log2()) {
     n0_ += 1;
     n += 1;
     storeVal[n] = logFunction(n0_, p);
     if (n >= maxIter) break;
   }
-  temp = storeVal[n];
-  storeVal[n] = temp - log_diff_exp(0, temp - storeVal[n - 1]) - log2();
-  storeVal[n + 1] = temp - log1mL - log2();
+
+  storeVal[n + 1] = storeVal[n] - log(-expm1(storeVal[n] - storeVal[n - 1])) - log2();
   return {log_sum_exp(sort_asc(storeVal[:(n + 1)])), 1. * n};
-}
+}