From a8038d6fb14371c7a32808f93ec6f5fcf8430d62 Mon Sep 17 00:00:00 2001
From: jonlachmann <jon@lachmann.nu>
Date: Mon, 6 Mar 2023 18:51:17 +0100
Subject: [PATCH 1/3] Add debug code to find a strange error...

---
 R/mjmcmc.R | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/R/mjmcmc.R b/R/mjmcmc.R
index 80bb0b2..1542487 100644
--- a/R/mjmcmc.R
+++ b/R/mjmcmc.R
@@ -201,6 +201,13 @@ mjmcmc.prop <- function (data, loglik.pi, model.cur, complex, pip_estimate, prob
   # Calculate acceptance probability for proposed model
   proposal$alpha <- min(0, (proposal$crit + model.cur$prob) - (model.cur$crit + proposal$prob))
 
+  ### DEBUG CODE
+  if (is.null(proposal$alpha) || is.na(proposal$alpha) || is.nan(proposal$alpha)) {
+    state <- list(mget(ls()))
+    number <- sample.int(1, 10000)
+    save(state, file = paste0("state", number, ".RData"))
+  }
+
   ### Format results and return them
   proposal$swap <- NULL; proposal$S <- NULL
   proposal$coefs <- proposal.res$coefs
-- 
GitLab


From a8317a9fe601039fa7cda79473694a048c9b23d8 Mon Sep 17 00:00:00 2001
From: jonlachmann <jon@lachmann.nu>
Date: Mon, 6 Mar 2023 20:33:32 +0100
Subject: [PATCH 2/3] ...

---
 R/mjmcmc.R | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/R/mjmcmc.R b/R/mjmcmc.R
index 1542487..a14c591 100644
--- a/R/mjmcmc.R
+++ b/R/mjmcmc.R
@@ -204,7 +204,7 @@ mjmcmc.prop <- function (data, loglik.pi, model.cur, complex, pip_estimate, prob
   ### DEBUG CODE
   if (is.null(proposal$alpha) || is.na(proposal$alpha) || is.nan(proposal$alpha)) {
     state <- list(mget(ls()))
-    number <- sample.int(1, 10000)
+    number <- sample.int(10000, 1)
     save(state, file = paste0("state", number, ".RData"))
   }
 
-- 
GitLab


From 04f5aed8a3e53dd5ee10bd46be397c38918686da Mon Sep 17 00:00:00 2001
From: jonlachmann <jon@lachmann.nu>
Date: Mon, 6 Mar 2023 20:55:49 +0100
Subject: [PATCH 3/3] Alpha cannot be calculated if the current and proposed
 models have crit which are -Inf or Inf Bug fixed.

---
 R/gmjmcmc_support.R |  5 +++++
 R/mjmcmc.R          | 11 ++---------
 2 files changed, 7 insertions(+), 9 deletions(-)

diff --git a/R/gmjmcmc_support.R b/R/gmjmcmc_support.R
index ce54548..1d63a51 100644
--- a/R/gmjmcmc_support.R
+++ b/R/gmjmcmc_support.R
@@ -84,6 +84,11 @@ loglik.pre <- function (loglik.pi, model, complex, data, params) {
   model.res <- loglik.pi(data[,1], data[,-1], c(T,model), complex, params)
   # Check that the critical value is acceptable
   if (!is.numeric(model.res$crit) || is.nan(model.res$crit)) model.res$crit <- -.Machine$double.xmax
+  # Alpha cannot be calculated if the current and proposed models have crit which are -Inf or Inf
+  if (is.infinite(model.res$crit)) {
+    if (model.res$crit > 0)  model.res$crit <- .Machine$double.xmax
+    else model.res$crit <- -.Machine$double.xmax
+  }
   return(model.res)
 }
 
diff --git a/R/mjmcmc.R b/R/mjmcmc.R
index a14c591..003a792 100644
--- a/R/mjmcmc.R
+++ b/R/mjmcmc.R
@@ -179,7 +179,7 @@ mjmcmc.prop <- function (data, loglik.pi, model.cur, complex, pip_estimate, prob
     # Select MH kernel
     q.g <- sample.int(n = 6, size = 1, prob = probs$mh)
     # Generate the proposal
-    proposal <- gen.proposal(model.cur$model, params$mh, q.g, NULL, pip_estimate, prob=T)
+    proposal <- gen.proposal(model.cur$model, params$mh, q.g, NULL, pip_estimate, prob = T)
     proposal$model <- xor(proposal$swap, model.cur$model)
 
     # Calculate current model probability given proposal
@@ -201,15 +201,8 @@ mjmcmc.prop <- function (data, loglik.pi, model.cur, complex, pip_estimate, prob
   # Calculate acceptance probability for proposed model
   proposal$alpha <- min(0, (proposal$crit + model.cur$prob) - (model.cur$crit + proposal$prob))
 
-  ### DEBUG CODE
-  if (is.null(proposal$alpha) || is.na(proposal$alpha) || is.nan(proposal$alpha)) {
-    state <- list(mget(ls()))
-    number <- sample.int(10000, 1)
-    save(state, file = paste0("state", number, ".RData"))
-  }
-
   ### Format results and return them
   proposal$swap <- NULL; proposal$S <- NULL
   proposal$coefs <- proposal.res$coefs
   return(proposal)
-}
\ No newline at end of file
+}
-- 
GitLab