From b048dd0d3b436142b286113fb4f288203a202a00 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Sun, 22 Jun 2025 12:01:05 -0700 Subject: [PATCH 1/5] Add info on target node when issuing logProb warning in sampling. --- packages/nimble/R/MCMC_samplers.R | 62 ++++++++++----------- packages/nimble/R/MCMC_utils.R | 6 +- packages/nimble/inst/CppCode/Utils.cpp | 4 +- packages/nimble/inst/include/nimble/Utils.h | 4 +- 4 files changed, 38 insertions(+), 38 deletions(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 3914f25ce..81fa12023 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -124,13 +124,13 @@ sampler_binary <- nimbleFunction( if(!model$isBinary(target)) stop('can only use binary sampler on discrete 0/1 (binary) nodes') }, run = function() { - currentLogProb <- checkLogProb(model$getLogProb(calcNodes)) + currentLogProb <- checkLogProb(model$getLogProb(calcNodes), target) model[[target]] <<- 1 - model[[target]] - otherLogProbPrior <- checkLogProb(model$calculate(target)) + otherLogProbPrior <- checkLogProb(model$calculate(target), target) if(otherLogProbPrior == -Inf) { otherLogProb <- otherLogProbPrior } else { - otherLogProb <- otherLogProbPrior + checkLogProb(model$calculate(calcNodesNoSelf)) + otherLogProb <- otherLogProbPrior + checkLogProb(model$calculate(calcNodesNoSelf), target) } if(currentLogProb == -Inf & otherLogProb == -Inf) stop("in binary sampler, all log probability density values are negative infinity and sampling cannot proceed") @@ -187,15 +187,15 @@ sampler_categorical <- nimbleFunction( }, run = function() { currentValue <- model[[target]] - logProbs[currentValue] <<- checkLogProb(model$getLogProb(calcNodes)) + logProbs[currentValue] <<- checkLogProb(model$getLogProb(calcNodes), target) for(i in 1:k) { if(i != currentValue) { model[[target]] <<- i - logProbPrior <- checkLogProb(model$calculate(target)) + logProbPrior <- checkLogProb(model$calculate(target), target) if(logProbPrior == -Inf) { logProbs[i] <<- -Inf } else { - logProbs[i] <<- logProbPrior + checkLogProb(model$calculate(calcNodesNoSelf)) + logProbs[i] <<- logProbPrior + checkLogProb(model$calculate(calcNodesNoSelf), target) } } } @@ -363,12 +363,12 @@ sampler_RW <- nimbleFunction( } } model[[target]] <<- propValue - logMHR <- checkLogProb(model$calculateDiff(target)) + logMHR <- checkLogProb(model$calculateDiff(target), target) if(logMHR == -Inf) { jump <- FALSE nimCopy(from = mvSaved, to = model, row = 1, nodes = target, logProb = TRUE) } else { - logMHR <- logMHR + checkLogProb(model$calculateDiff(calcNodesNoSelf)) + propLogScale + logMHR <- logMHR + checkLogProb(model$calculateDiff(calcNodesNoSelf), target) + propLogScale jump <- decide(logMHR) if(jump) { ##model$calculate(calcNodesPPomitted) @@ -520,14 +520,14 @@ sampler_RW_noncentered <- nimbleFunction( } model[[target]] <<- propValue - logMHR <- checkLogProb(model$calculateDiff(target)) + logMHR <- checkLogProb(model$calculateDiff(target), target) if(logMHR == -Inf) { jump <- FALSE nimCopy(from = mvSaved, to = model, row = 1, nodes = target, logProb = TRUE) } else { ## Shift effects and add log-determinant of Jacobian of transformation. logMHR <- logMHR + updateNoncentered(propValue, currentValue) - logMHR <- logMHR + checkLogProb(model$calculateDiff(calcNodesNoSelf)) + propLogScale + logMHR <- logMHR + checkLogProb(model$calculateDiff(calcNodesNoSelf), target) + propLogScale jump <- decide(logMHR) if(jump) { ##model$calculate(calcNodesPPomitted) @@ -689,12 +689,12 @@ sampler_RW_block <- nimbleFunction( for(i in 1:tries) { propValueVector <- generateProposalVector() values(model, targetAsScalar) <<- propValueVector - lpD <- checkLogProb(model$calculateDiff(calcNodesProposalStage)) + lpD <- checkLogProb(model$calculateDiff(calcNodesProposalStage), target) if(lpD == -Inf) { jump <- FALSE nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodesProposalStage, logProb = TRUE) } else { - lpD <- lpD + checkLogProb(model$calculateDiff(calcNodesDepStage)) + lpD <- lpD + checkLogProb(model$calculateDiff(calcNodesDepStage), target) jump <- decide(lpD) if(jump) { ##model$calculate(calcNodesPPomitted) @@ -941,9 +941,9 @@ sampler_slice <- nimbleFunction( setAndCalculateTarget = function(value = double()) { if(discrete) value <- floor(value) model[[target]] <<- value - lp <- checkLogProb(model$calculate(target)) + lp <- checkLogProb(model$calculate(target), target) if(lp == -Inf) return(-Inf) - lp <- lp + checkLogProb(model$calculate(calcNodesNoSelf)) + lp <- lp + checkLogProb(model$calculate(calcNodesNoSelf), target) returnType(double()) return(lp) }, @@ -1082,10 +1082,10 @@ sampler_slice_noncentered <- nimbleFunction( setAndCalculateTarget = function(value = double()) { if(discrete) value <- floor(value) model[[target]] <<- value - lp <- checkLogProb(model$calculate(target)) + lp <- checkLogProb(model$calculate(target), target) if(lp == -Inf) return(-Inf) lp <- lp + updateNoncentered(value) - lp <- lp + checkLogProb(model$calculate(calcNodesNoSelf)) + lp <- lp + checkLogProb(model$calculate(calcNodesNoSelf), target) returnType(double()) return(lp) }, @@ -1215,7 +1215,7 @@ sampler_ess <- nimbleFunction( theta_min <- theta - 2*Pi theta_max <- theta values(model, target) <<- f[1:d]*cos(theta) + nu[1:d]*sin(theta) + target_mean[1:d] - lp <- checkLogProb(model$calculate(calcNodesNoSelf)) + lp <- checkLogProb(model$calculate(calcNodesNoSelf), target) numContractions <- 0 while((is.nan(lp) | lp < u) & theta_max - theta_min > eps & numContractions < maxContractions) { # must be is.nan() ## The checks for theta_max - theta_min small and max number of contractions are @@ -1224,7 +1224,7 @@ sampler_ess <- nimbleFunction( if(theta < 0) theta_min <- theta else theta_max <- theta theta <- runif(1, theta_min, theta_max) values(model, target) <<- f[1:d]*cos(theta) + nu[1:d]*sin(theta) + target_mean[1:d] - lp <- checkLogProb(model$calculate(calcNodesNoSelf)) + lp <- checkLogProb(model$calculate(calcNodesNoSelf), target) numContractions <- numContractions + 1 } if(theta_max - theta_min <= eps | numContractions == maxContractions) { @@ -1370,9 +1370,9 @@ sampler_AF_slice <- nimbleFunction( for(i in 1:d) if(discrete[i] == 1) targetValues[i] <- floor(targetValues[i]) values(model, target) <<- targetValues - lp <- checkLogProb(model$calculate(calcNodesProposalStage)) + lp <- checkLogProb(model$calculate(calcNodesProposalStage), target) if(lp == -Inf) return(lp) - lp <- lp + checkLogProb(model$calculate(calcNodesDepStage)) + lp <- lp + checkLogProb(model$calculate(calcNodesDepStage), target) returnType(double()) return(lp) }, @@ -1499,7 +1499,7 @@ sampler_crossLevel <- nimbleFunction( propLP0 <- 0 for(iSF in seq_along(lowConjugateGetLogDensityFunctions)) { propLP0 <- propLP0 + lowConjugateGetLogDensityFunctions[[iSF]]$run() } propValueVector <- topRWblockSamplerFunction$generateProposalVector() - topLP <- checkLogProb(my_setAndCalculateTop$run(propValueVector)) + topLP <- checkLogProb(my_setAndCalculateTop$run(propValueVector), target) if(topLP == -Inf) { nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE) } @@ -1510,7 +1510,7 @@ sampler_crossLevel <- nimbleFunction( propLP1 <- 0 for(iSF in seq_along(lowConjugateGetLogDensityFunctions)) propLP1 <- propLP1 + lowConjugateGetLogDensityFunctions[[iSF]]$run() - logMHR <- checkLogProb(modelLP1) - checkLogProb(modelLP0) - checkLogProb(propLP1) + checkLogProb(propLP0) + logMHR <- checkLogProb(modelLP1, target) - checkLogProb(modelLP0, target) - checkLogProb(propLP1, target) + checkLogProb(propLP0, target) jump <- decide(logMHR) if(jump) { nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) @@ -1722,7 +1722,7 @@ sampler_RW_dirichlet <- nimbleFunction( thetaVecProp <- thetaVec thetaVecProp[i] <- propValue values(model, target) <<- thetaVecProp / sum(thetaVecProp) - logMHR <- alphaVec[i]*propLogScale + currentValue - propValue + checkLogProb(model$calculateDiff(calcNodesNoSelf)) + logMHR <- alphaVec[i]*propLogScale + currentValue - propValue + checkLogProb(model$calculateDiff(calcNodesNoSelf), target) jump <- decide(logMHR) } else jump <- FALSE if(adaptive & jump) timesAcceptedVec[i] <<- timesAcceptedVec[i] + 1 @@ -1842,7 +1842,7 @@ sampler_RW_wishart <- nimbleFunction( ## matrix multiply to get proposal value (matrix) model[[target]] <<- t(propValue_chol) %*% propValue_chol ## decide and jump - logMHR <- checkLogProb(model$calculateDiff(calcNodes)) + logMHR <- checkLogProb(model$calculateDiff(calcNodes), target) deltaDiag <- thetaVec_prop[1:d]-thetaVec[1:d] for(i in 1:d) logMHR <- logMHR + (d+2-i)*deltaDiag[i] ## took me quite a while to derive this jump <- decide(logMHR) @@ -1965,8 +1965,8 @@ sampler_RW_lkj_corr_cholesky <- nimbleFunction( propValue[jprime] <<- z[jprime, i] * sqrt(partialSumsProp[jprime]) } model[[target]][j:i, i] <<- propValue[j:i] - logMHR <- checkLogProb(calculateDiff(model, calcNodesNoSelf)) + - checkLogProb(calculateDiff(model, target)) + logMHR <- checkLogProb(calculateDiff(model, calcNodesNoSelf), target) + + checkLogProb(calculateDiff(model, target), target) ## Adjust MHR to account for non-symmetric proposal by adjusting prior on U to transformed scale (i.e., y). ## cosh component is for dz/dy and other component is for du/dz where 'u' is the corr matrix. ## This follows Stan reference manual Section 10.12 (for version 2.27). @@ -2138,12 +2138,12 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction( ## Adjust for log determinant term from initial values logMHR <- logMHR - logDetJac - lpD <- checkLogProb(calculateDiff(model, calcNodesProposalStage)) + lpD <- checkLogProb(calculateDiff(model, calcNodesProposalStage), target) if(lpD == -Inf) { nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodesProposalStage, logProb = TRUE) jump <- FALSE } else { - logMHR <- logMHR + lpD + checkLogProb(calculateDiff(model, calcNodesDepStage)) + logMHR <- logMHR + lpD + checkLogProb(calculateDiff(model, calcNodesDepStage), target) jump <- decide(logMHR) if(jump) { nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) @@ -2418,7 +2418,7 @@ CAR_scalar_RW <- nimbleFunction( propValue <- rnorm(1, mean = model[[targetScalar]], sd = scale) model[[targetScalar]] <<- propValue lp1 <- dcarList[[1]]$run() + model$calculate(depNodes) - logMHR <- checkLogProb(lp1) - checkLogProb(lp0) + logMHR <- checkLogProb(lp1, target) - checkLogProb(lp0, target) jump <- decide(logMHR) if(jump) { model$calculate(targetScalar) @@ -3359,8 +3359,8 @@ sampler_barker <- nimbleFunction( gradProposed <<- gradient(proposal) newLogDetJacobian <- my_parameterTransform$logDetJacobian(proposal) - lpD <- checkLogProb(model$calculateDiff(calcNodes)) + checkLogProb(newLogDetJacobian) - - checkLogProb(oldLogDetJacobian) + checkLogProb(calculateLogHastingsRatio(diff)) + lpD <- checkLogProb(model$calculateDiff(calcNodes), target) + checkLogProb(newLogDetJacobian, target) - + checkLogProb(oldLogDetJacobian, target) + checkLogProb(calculateLogHastingsRatio(diff), target) jump <- decide(lpD) if(jump) nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) diff --git a/packages/nimble/R/MCMC_utils.R b/packages/nimble/R/MCMC_utils.R index 07cf58c96..eddfcf28c 100644 --- a/packages/nimble/R/MCMC_utils.R +++ b/packages/nimble/R/MCMC_utils.R @@ -56,7 +56,7 @@ decideAndJump <- nimbleFunction( }, run = function(modelLP1 = double(), modelLP0 = double(), propLP1 = double(), propLP0 = double()) { ## Check each one individually to catch case like `3 - Inf`. - logMHR <- checkLogProb(modelLP1) - checkLogProb(modelLP0) - checkLogProb(propLP1) + checkLogProb(propLP0) + logMHR <- checkLogProb(modelLP1, target) - checkLogProb(modelLP0, target) - checkLogProb(propLP1, target) + checkLogProb(propLP0, target) jump <- decide(logMHR) if(jump) { nimCopy(from = model, to = mvSaved, row = 1, nodes = target, logProb = TRUE) @@ -72,11 +72,11 @@ decideAndJump <- nimbleFunction( } ) -checkLogProb <- function(logProb) { +checkLogProb <- function(logProb, target) { if(is.na(logProb)) return(-Inf) if(logProb == Inf) - print("MCMC sampling encountered a log probability density value of infinity. Results of sampling may not be valid.") + print("MCMC sampling of ", target, " encountered a log probability density value of infinity. Results of sampling may not be valid.") return(logProb) } diff --git a/packages/nimble/inst/CppCode/Utils.cpp b/packages/nimble/inst/CppCode/Utils.cpp index 0962815be..3ead1ace4 100644 --- a/packages/nimble/inst/CppCode/Utils.cpp +++ b/packages/nimble/inst/CppCode/Utils.cpp @@ -63,8 +63,8 @@ bool decide(double lMHr) { // simple function accept or reject based on log Metr return(false); } -void checkLogProbWarn() { - _nimble_global_output<<"MCMC sampling encountered a log probability density value of infinity. Results of sampling may not be valid.\n"; +void checkLogProbWarn(string target) { + _nimble_global_output<<"MCMC sampling of " << target << " encountered a log probability density value of infinity. Results of sampling may not be valid.\n"; nimble_print_to_R(_nimble_global_output); } diff --git a/packages/nimble/inst/include/nimble/Utils.h b/packages/nimble/inst/include/nimble/Utils.h index f46cfc19e..eb382b45a 100644 --- a/packages/nimble/inst/include/nimble/Utils.h +++ b/packages/nimble/inst/include/nimble/Utils.h @@ -162,9 +162,9 @@ class nimbleTimerClass_ { bool decide(double lMHr); //void allocate(vector< vector > *vv, int sampleSize, int variableSize); -void checkLogProbWarn(); +void checkLogProbWarn(string target); -static inline double checkLogProb(double logProb) { if(ISNAN(logProb)) return(-std::numeric_limits::infinity()); if(logProb == std::numeric_limits::infinity()) checkLogProbWarn(); return(logProb);} +static inline double checkLogProb(double logProb, string target) { if(ISNAN(logProb)) return(-std::numeric_limits::infinity()); if(logProb == std::numeric_limits::infinity()) checkLogProbWarn(target); return(logProb);} void nimStop(string msg); From 4019e2458b95f096b30e65eed23e925b3516bb41 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Wed, 25 Jun 2025 12:04:08 -0700 Subject: [PATCH 2/5] Fix use of target in checkLogProb case. --- packages/nimble/R/MCMC_samplers.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 81fa12023..81428697b 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -2418,7 +2418,7 @@ CAR_scalar_RW <- nimbleFunction( propValue <- rnorm(1, mean = model[[targetScalar]], sd = scale) model[[targetScalar]] <<- propValue lp1 <- dcarList[[1]]$run() + model$calculate(depNodes) - logMHR <- checkLogProb(lp1, target) - checkLogProb(lp0, target) + logMHR <- checkLogProb(lp1, targetScalar) - checkLogProb(lp0, targetScalar) jump <- decide(logMHR) if(jump) { model$calculate(targetScalar) From b2b8a29ddced7be6258134720bd36b74a4300749 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Wed, 25 Jun 2025 13:06:17 -0700 Subject: [PATCH 3/5] Fix checkLogProb messaging. --- packages/nimble/R/MCMC_utils.R | 2 +- packages/nimble/inst/CppCode/Utils.cpp | 2 +- packages/nimble/inst/include/nimble/Utils.h | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/packages/nimble/R/MCMC_utils.R b/packages/nimble/R/MCMC_utils.R index eddfcf28c..c043331d8 100644 --- a/packages/nimble/R/MCMC_utils.R +++ b/packages/nimble/R/MCMC_utils.R @@ -76,7 +76,7 @@ checkLogProb <- function(logProb, target) { if(is.na(logProb)) return(-Inf) if(logProb == Inf) - print("MCMC sampling of ", target, " encountered a log probability density value of infinity. Results of sampling may not be valid.") + cat("MCMC sampling of ", target, " encountered a log probability density value of infinity. Results of sampling may not be valid.\n") return(logProb) } diff --git a/packages/nimble/inst/CppCode/Utils.cpp b/packages/nimble/inst/CppCode/Utils.cpp index 3ead1ace4..b02968fc1 100644 --- a/packages/nimble/inst/CppCode/Utils.cpp +++ b/packages/nimble/inst/CppCode/Utils.cpp @@ -63,7 +63,7 @@ bool decide(double lMHr) { // simple function accept or reject based on log Metr return(false); } -void checkLogProbWarn(string target) { +void checkLogProbWarn(std::string target) { _nimble_global_output<<"MCMC sampling of " << target << " encountered a log probability density value of infinity. Results of sampling may not be valid.\n"; nimble_print_to_R(_nimble_global_output); } diff --git a/packages/nimble/inst/include/nimble/Utils.h b/packages/nimble/inst/include/nimble/Utils.h index eb382b45a..2152be672 100644 --- a/packages/nimble/inst/include/nimble/Utils.h +++ b/packages/nimble/inst/include/nimble/Utils.h @@ -162,9 +162,9 @@ class nimbleTimerClass_ { bool decide(double lMHr); //void allocate(vector< vector > *vv, int sampleSize, int variableSize); -void checkLogProbWarn(string target); +void checkLogProbWarn(std::string target); -static inline double checkLogProb(double logProb, string target) { if(ISNAN(logProb)) return(-std::numeric_limits::infinity()); if(logProb == std::numeric_limits::infinity()) checkLogProbWarn(target); return(logProb);} +static inline double checkLogProb(double logProb, std::string target) { if(ISNAN(logProb)) return(-std::numeric_limits::infinity()); if(logProb == std::numeric_limits::infinity()) checkLogProbWarn(target); return(logProb);} void nimStop(string msg); From dfbd56dd244c896e1047b94224679d212c7ca2a7 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Thu, 26 Jun 2025 07:47:29 -0700 Subject: [PATCH 4/5] Fix a multi-node case in checkLogProb. --- packages/nimble/R/MCMC_utils.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/packages/nimble/R/MCMC_utils.R b/packages/nimble/R/MCMC_utils.R index c043331d8..d1d3a64ae 100644 --- a/packages/nimble/R/MCMC_utils.R +++ b/packages/nimble/R/MCMC_utils.R @@ -53,10 +53,11 @@ decideAndJump <- nimbleFunction( setup = function(model, mvSaved, target, UNUSED) { ## should remove UNUSED argument, after next release of nimbleSMC -DT July 2024 ccList <- mcmc_determineCalcAndCopyNodes(model, target) copyNodesDeterm <- ccList$copyNodesDeterm; copyNodesStoch <- ccList$copyNodesStoch # not used: calcNodes, calcNodesNoSelf + targetOne <- target[1] }, run = function(modelLP1 = double(), modelLP0 = double(), propLP1 = double(), propLP0 = double()) { ## Check each one individually to catch case like `3 - Inf`. - logMHR <- checkLogProb(modelLP1, target) - checkLogProb(modelLP0, target) - checkLogProb(propLP1, target) + checkLogProb(propLP0, target) + logMHR <- checkLogProb(modelLP1, targetOne) - checkLogProb(modelLP0, targetOne) - checkLogProb(propLP1, targetOne) + checkLogProb(propLP0, targetOne) jump <- decide(logMHR) if(jump) { nimCopy(from = model, to = mvSaved, row = 1, nodes = target, logProb = TRUE) From 3c72628dfbebe8190fde34e3562b3978ee2e7911 Mon Sep 17 00:00:00 2001 From: Christopher Paciorek Date: Thu, 26 Jun 2025 09:34:57 -0700 Subject: [PATCH 5/5] Fix checkLogProb msg about target name for multi-node cases. --- packages/nimble/R/MCMC_samplers.R | 41 +++++++++++++++++++------------ packages/nimble/R/MCMC_utils.R | 10 +++++--- 2 files changed, 32 insertions(+), 19 deletions(-) diff --git a/packages/nimble/R/MCMC_samplers.R b/packages/nimble/R/MCMC_samplers.R index 81428697b..a399b76f1 100644 --- a/packages/nimble/R/MCMC_samplers.R +++ b/packages/nimble/R/MCMC_samplers.R @@ -684,17 +684,18 @@ sampler_RW_block <- nimbleFunction( if(!all(dim(propCov) == d)) stop('propCov matrix must have dimension ', d, 'x', d, '\n') if(!isSymmetric(propCov)) stop('propCov matrix must be symmetric') if(adaptInterval < 2) stop('sampler_RW_block: `adaptInterval` must be at least 2') + targetNames <- createNamesString(target) }, run = function() { for(i in 1:tries) { propValueVector <- generateProposalVector() values(model, targetAsScalar) <<- propValueVector - lpD <- checkLogProb(model$calculateDiff(calcNodesProposalStage), target) + lpD <- checkLogProb(model$calculateDiff(calcNodesProposalStage), targetNames) if(lpD == -Inf) { jump <- FALSE nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodesProposalStage, logProb = TRUE) } else { - lpD <- lpD + checkLogProb(model$calculateDiff(calcNodesDepStage), target) + lpD <- lpD + checkLogProb(model$calculateDiff(calcNodesDepStage), targetNames) jump <- decide(lpD) if(jump) { ##model$calculate(calcNodesPPomitted) @@ -1204,6 +1205,7 @@ sampler_ess <- nimbleFunction( ## checks if(length(target) > 1) stop('elliptical slice sampler only applies to one target node') if(!(model$getDistribution(target) %in% c('dnorm', 'dmnorm'))) stop('elliptical slice sampler only applies to normal distributions') + targetNames <- createNamesString(target) }, run = function() { u <- model$getLogProb(calcNodesNoSelf) - rexp(1, 1) @@ -1215,7 +1217,7 @@ sampler_ess <- nimbleFunction( theta_min <- theta - 2*Pi theta_max <- theta values(model, target) <<- f[1:d]*cos(theta) + nu[1:d]*sin(theta) + target_mean[1:d] - lp <- checkLogProb(model$calculate(calcNodesNoSelf), target) + lp <- checkLogProb(model$calculate(calcNodesNoSelf), targetNames) numContractions <- 0 while((is.nan(lp) | lp < u) & theta_max - theta_min > eps & numContractions < maxContractions) { # must be is.nan() ## The checks for theta_max - theta_min small and max number of contractions are @@ -1224,7 +1226,7 @@ sampler_ess <- nimbleFunction( if(theta < 0) theta_min <- theta else theta_max <- theta theta <- runif(1, theta_min, theta_max) values(model, target) <<- f[1:d]*cos(theta) + nu[1:d]*sin(theta) + target_mean[1:d] - lp <- checkLogProb(model$calculate(calcNodesNoSelf), target) + lp <- checkLogProb(model$calculate(calcNodesNoSelf), targetNames) numContractions <- numContractions + 1 } if(theta_max - theta_min <= eps | numContractions == maxContractions) { @@ -1300,6 +1302,7 @@ sampler_AF_slice <- nimbleFunction( stop('sliceWidths must be a numeric vector') if(length(widthVec) != d) stop('sliceWidths must have length = ', d) if(adaptFactorInterval < 2) stop('sampler_AF_slice: `adaptFactorInterval` must be at least 2') + targetNames <- createNamesString(target) }, run = function() { maxContractionsReached <- FALSE @@ -1370,9 +1373,9 @@ sampler_AF_slice <- nimbleFunction( for(i in 1:d) if(discrete[i] == 1) targetValues[i] <- floor(targetValues[i]) values(model, target) <<- targetValues - lp <- checkLogProb(model$calculate(calcNodesProposalStage), target) + lp <- checkLogProb(model$calculate(calcNodesProposalStage), targetNames) if(lp == -Inf) return(lp) - lp <- lp + checkLogProb(model$calculate(calcNodesDepStage), target) + lp <- lp + checkLogProb(model$calculate(calcNodesDepStage), targetNames) returnType(double()) return(lp) }, @@ -1493,13 +1496,14 @@ sampler_crossLevel <- nimbleFunction( lowConjugateGetLogDensityFunctions[[iLN]] <- getPosteriorDensityFromConjSampler(lowConjugateSamplerFunctions[[iLN]]) } my_setAndCalculateTop <- setAndCalculate(model, target) + targetNames <- createNamesString(target) }, run = function() { modelLP0 <- model$getLogProb(calcNodes) propLP0 <- 0 for(iSF in seq_along(lowConjugateGetLogDensityFunctions)) { propLP0 <- propLP0 + lowConjugateGetLogDensityFunctions[[iSF]]$run() } propValueVector <- topRWblockSamplerFunction$generateProposalVector() - topLP <- checkLogProb(my_setAndCalculateTop$run(propValueVector), target) + topLP <- checkLogProb(my_setAndCalculateTop$run(propValueVector), targetNames) if(topLP == -Inf) { nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE) } @@ -1510,7 +1514,7 @@ sampler_crossLevel <- nimbleFunction( propLP1 <- 0 for(iSF in seq_along(lowConjugateGetLogDensityFunctions)) propLP1 <- propLP1 + lowConjugateGetLogDensityFunctions[[iSF]]$run() - logMHR <- checkLogProb(modelLP1, target) - checkLogProb(modelLP0, target) - checkLogProb(propLP1, target) + checkLogProb(propLP0, target) + logMHR <- checkLogProb(modelLP1, targetNames) - checkLogProb(modelLP0, targetNames) - checkLogProb(propLP1, targetNames) + checkLogProb(propLP0, targetNames) jump <- decide(logMHR) if(jump) { nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) @@ -1710,6 +1714,7 @@ sampler_RW_dirichlet <- nimbleFunction( ## checks if(length(model$expandNodeNames(target)) > 1) stop('RW_dirichlet sampler only applies to one target node') if(model$getDistribution(target) != 'ddirch') stop('can only use RW_dirichlet sampler for dirichlet distributions') + targetNames <- createNamesString(target) }, run = function() { if(thetaVec[1] == 0) thetaVec <<- values(model, target) ## initialization @@ -1722,7 +1727,7 @@ sampler_RW_dirichlet <- nimbleFunction( thetaVecProp <- thetaVec thetaVecProp[i] <- propValue values(model, target) <<- thetaVecProp / sum(thetaVecProp) - logMHR <- alphaVec[i]*propLogScale + currentValue - propValue + checkLogProb(model$calculateDiff(calcNodesNoSelf), target) + logMHR <- alphaVec[i]*propLogScale + currentValue - propValue + checkLogProb(model$calculateDiff(calcNodesNoSelf), targetNames) jump <- decide(logMHR) } else jump <- FALSE if(adaptive & jump) timesAcceptedVec[i] <<- timesAcceptedVec[i] + 1 @@ -1815,6 +1820,7 @@ sampler_RW_wishart <- nimbleFunction( if(!all(dim(propCov) == nTheta)) stop('propCov matrix must have dimension ', d, 'x', d) if(!isSymmetric(propCov)) stop('propCov matrix must be symmetric') if(adaptInterval < 2) stop('sampler_RW_wishart: `adaptInterval` must be at least 2') + targetNames <- createNamesString(target) }, run = function() { currentValue <<- model[[target]] @@ -1842,7 +1848,7 @@ sampler_RW_wishart <- nimbleFunction( ## matrix multiply to get proposal value (matrix) model[[target]] <<- t(propValue_chol) %*% propValue_chol ## decide and jump - logMHR <- checkLogProb(model$calculateDiff(calcNodes), target) + logMHR <- checkLogProb(model$calculateDiff(calcNodes), targetNames) deltaDiag <- thetaVec_prop[1:d]-thetaVec[1:d] for(i in 1:d) logMHR <- logMHR + (d+2-i)*deltaDiag[i] ## took me quite a while to derive this jump <- decide(logMHR) @@ -1940,6 +1946,7 @@ sampler_RW_lkj_corr_cholesky <- nimbleFunction( if(dist != 'dlkj_corr_cholesky') stop('RW_lkj_corr_cholesky sampler can only be used with the dlkj_corr_cholesky distribution.') if(d < 2) stop('RW_lkj_corr_cholesky sampler requires target node dimension to be at least 2x2.') if(adaptFactorExponent < 0) stop('Cannot use RW_lkj_corr_cholesky sampler with adaptFactorExponent control parameter less than 0.') + targetNames <- createNamesString(target) }, run = function() { ## calculate transformed values (in unconstrained space) and partial sums in each column @@ -1965,8 +1972,8 @@ sampler_RW_lkj_corr_cholesky <- nimbleFunction( propValue[jprime] <<- z[jprime, i] * sqrt(partialSumsProp[jprime]) } model[[target]][j:i, i] <<- propValue[j:i] - logMHR <- checkLogProb(calculateDiff(model, calcNodesNoSelf), target) + - checkLogProb(calculateDiff(model, target), target) + logMHR <- checkLogProb(calculateDiff(model, calcNodesNoSelf), targetNames) + + checkLogProb(calculateDiff(model, target), targetNames) ## Adjust MHR to account for non-symmetric proposal by adjusting prior on U to transformed scale (i.e., y). ## cosh component is for dz/dy and other component is for du/dz where 'u' is the corr matrix. ## This follows Stan reference manual Section 10.12 (for version 2.27). @@ -2102,6 +2109,7 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction( if(adaptFactorExponent < 0) stop('Cannot use RW_block_lkj_corr_cholesky sampler with adaptFactorExponent control parameter less than 0.') if(adaptInterval < 2) stop('sampler_RW_block_lkj_corr_cholesky: `adaptInterval` must be at least 2') + targetNames <- createNamesString(target) }, run = function() { transform(model[[target]]) # compute z and partialSums @@ -2138,12 +2146,12 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction( ## Adjust for log determinant term from initial values logMHR <- logMHR - logDetJac - lpD <- checkLogProb(calculateDiff(model, calcNodesProposalStage), target) + lpD <- checkLogProb(calculateDiff(model, calcNodesProposalStage), targetNames) if(lpD == -Inf) { nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodesProposalStage, logProb = TRUE) jump <- FALSE } else { - logMHR <- logMHR + lpD + checkLogProb(calculateDiff(model, calcNodesDepStage), target) + logMHR <- logMHR + lpD + checkLogProb(calculateDiff(model, calcNodesDepStage), targetNames) jump <- decide(logMHR) if(jump) { nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) @@ -3325,6 +3333,7 @@ sampler_barker <- nimbleFunction( ## checks if(!isTRUE(nimbleOptions('enableDerivs'))) stop('must enable NIMBLE derivatives, set nimbleOptions(enableDerivs = TRUE)', call. = FALSE) if(!isTRUE(model$modelDef[['buildDerivs']])) stop('must set buildDerivs = TRUE when building model', call. = FALSE) + targetNames <- createNamesString(target) }, run = function() { current <- my_parameterTransform$transform(values(model, targetNodes)) @@ -3359,8 +3368,8 @@ sampler_barker <- nimbleFunction( gradProposed <<- gradient(proposal) newLogDetJacobian <- my_parameterTransform$logDetJacobian(proposal) - lpD <- checkLogProb(model$calculateDiff(calcNodes), target) + checkLogProb(newLogDetJacobian, target) - - checkLogProb(oldLogDetJacobian, target) + checkLogProb(calculateLogHastingsRatio(diff), target) + lpD <- checkLogProb(model$calculateDiff(calcNodes), targetNames) + checkLogProb(newLogDetJacobian, targetNames) - + checkLogProb(oldLogDetJacobian, targetNames) + checkLogProb(calculateLogHastingsRatio(diff), targetNames) jump <- decide(lpD) if(jump) nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE) diff --git a/packages/nimble/R/MCMC_utils.R b/packages/nimble/R/MCMC_utils.R index d1d3a64ae..65da14490 100644 --- a/packages/nimble/R/MCMC_utils.R +++ b/packages/nimble/R/MCMC_utils.R @@ -53,11 +53,11 @@ decideAndJump <- nimbleFunction( setup = function(model, mvSaved, target, UNUSED) { ## should remove UNUSED argument, after next release of nimbleSMC -DT July 2024 ccList <- mcmc_determineCalcAndCopyNodes(model, target) copyNodesDeterm <- ccList$copyNodesDeterm; copyNodesStoch <- ccList$copyNodesStoch # not used: calcNodes, calcNodesNoSelf - targetOne <- target[1] + targetNames <- createNamesString(target) }, run = function(modelLP1 = double(), modelLP0 = double(), propLP1 = double(), propLP0 = double()) { ## Check each one individually to catch case like `3 - Inf`. - logMHR <- checkLogProb(modelLP1, targetOne) - checkLogProb(modelLP0, targetOne) - checkLogProb(propLP1, targetOne) + checkLogProb(propLP0, targetOne) + logMHR <- checkLogProb(modelLP1, targetNames) - checkLogProb(modelLP0, targetNames) - checkLogProb(propLP1, targetNames) + checkLogProb(propLP0, targetNames) jump <- decide(logMHR) if(jump) { nimCopy(from = model, to = mvSaved, row = 1, nodes = target, logProb = TRUE) @@ -551,7 +551,11 @@ mcmc_checkTargetAD <- function(model, targetNodes, samplerType) { paste0(dists[!ADok], collapse = ', ')) } - +createNamesString <- function(target) { + if(length(target) == 1) return(target) + if(length(target) > 4) target <- c(target[1:4], "...") + return(paste0("{", paste(target, collapse = ', '), "}")) +}