Skip to content

Commit 7d2f839

Browse files
paciorekpaciorek
authored
tell user which node is being sampled when checkLogProb warning invoked (#1559)
* Add info on target node when issuing logProb warning in sampling. * Fix use of target in checkLogProb case. * Fix checkLogProb messaging. * Fix a multi-node case in checkLogProb. * Fix checkLogProb msg about target name for multi-node cases.
1 parent 5417a45 commit 7d2f839

4 files changed

Lines changed: 54 additions & 38 deletions

File tree

‎packages/nimble/R/MCMC_samplers.R‎

Lines changed: 40 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -124,13 +124,13 @@ sampler_binary <- nimbleFunction(
124124
if(!model$isBinary(target)) stop('can only use binary sampler on discrete 0/1 (binary) nodes')
125125
},
126126
run = function() {
127-
currentLogProb <- checkLogProb(model$getLogProb(calcNodes))
127+
currentLogProb <- checkLogProb(model$getLogProb(calcNodes), target)
128128
model[[target]] <<- 1 - model[[target]]
129-
otherLogProbPrior <- checkLogProb(model$calculate(target))
129+
otherLogProbPrior <- checkLogProb(model$calculate(target), target)
130130
if(otherLogProbPrior == -Inf) {
131131
otherLogProb <- otherLogProbPrior
132132
} else {
133-
otherLogProb <- otherLogProbPrior + checkLogProb(model$calculate(calcNodesNoSelf))
133+
otherLogProb <- otherLogProbPrior + checkLogProb(model$calculate(calcNodesNoSelf), target)
134134
}
135135
if(currentLogProb == -Inf & otherLogProb == -Inf)
136136
stop("in binary sampler, all log probability density values are negative infinity and sampling cannot proceed")
@@ -187,15 +187,15 @@ sampler_categorical <- nimbleFunction(
187187
},
188188
run = function() {
189189
currentValue <- model[[target]]
190-
logProbs[currentValue] <<- checkLogProb(model$getLogProb(calcNodes))
190+
logProbs[currentValue] <<- checkLogProb(model$getLogProb(calcNodes), target)
191191
for(i in 1:k) {
192192
if(i != currentValue) {
193193
model[[target]] <<- i
194-
logProbPrior <- checkLogProb(model$calculate(target))
194+
logProbPrior <- checkLogProb(model$calculate(target), target)
195195
if(logProbPrior == -Inf) {
196196
logProbs[i] <<- -Inf
197197
} else {
198-
logProbs[i] <<- logProbPrior + checkLogProb(model$calculate(calcNodesNoSelf))
198+
logProbs[i] <<- logProbPrior + checkLogProb(model$calculate(calcNodesNoSelf), target)
199199
}
200200
}
201201
}
@@ -363,12 +363,12 @@ sampler_RW <- nimbleFunction(
363363
}
364364
}
365365
model[[target]] <<- propValue
366-
logMHR <- checkLogProb(model$calculateDiff(target))
366+
logMHR <- checkLogProb(model$calculateDiff(target), target)
367367
if(logMHR == -Inf) {
368368
jump <- FALSE
369369
nimCopy(from = mvSaved, to = model, row = 1, nodes = target, logProb = TRUE)
370370
} else {
371-
logMHR <- logMHR + checkLogProb(model$calculateDiff(calcNodesNoSelf)) + propLogScale
371+
logMHR <- logMHR + checkLogProb(model$calculateDiff(calcNodesNoSelf), target) + propLogScale
372372
jump <- decide(logMHR)
373373
if(jump) {
374374
##model$calculate(calcNodesPPomitted)
@@ -520,14 +520,14 @@ sampler_RW_noncentered <- nimbleFunction(
520520
}
521521
model[[target]] <<- propValue
522522

523-
logMHR <- checkLogProb(model$calculateDiff(target))
523+
logMHR <- checkLogProb(model$calculateDiff(target), target)
524524
if(logMHR == -Inf) {
525525
jump <- FALSE
526526
nimCopy(from = mvSaved, to = model, row = 1, nodes = target, logProb = TRUE)
527527
} else {
528528
## Shift effects and add log-determinant of Jacobian of transformation.
529529
logMHR <- logMHR + updateNoncentered(propValue, currentValue)
530-
logMHR <- logMHR + checkLogProb(model$calculateDiff(calcNodesNoSelf)) + propLogScale
530+
logMHR <- logMHR + checkLogProb(model$calculateDiff(calcNodesNoSelf), target) + propLogScale
531531
jump <- decide(logMHR)
532532
if(jump) {
533533
##model$calculate(calcNodesPPomitted)
@@ -684,17 +684,18 @@ sampler_RW_block <- nimbleFunction(
684684
if(!all(dim(propCov) == d)) stop('propCov matrix must have dimension ', d, 'x', d, '\n')
685685
if(!isSymmetric(propCov)) stop('propCov matrix must be symmetric')
686686
if(adaptInterval < 2) stop('sampler_RW_block: `adaptInterval` must be at least 2')
687+
targetNames <- createNamesString(target)
687688
},
688689
run = function() {
689690
for(i in 1:tries) {
690691
propValueVector <- generateProposalVector()
691692
values(model, targetAsScalar) <<- propValueVector
692-
lpD <- checkLogProb(model$calculateDiff(calcNodesProposalStage))
693+
lpD <- checkLogProb(model$calculateDiff(calcNodesProposalStage), targetNames)
693694
if(lpD == -Inf) {
694695
jump <- FALSE
695696
nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodesProposalStage, logProb = TRUE)
696697
} else {
697-
lpD <- lpD + checkLogProb(model$calculateDiff(calcNodesDepStage))
698+
lpD <- lpD + checkLogProb(model$calculateDiff(calcNodesDepStage), targetNames)
698699
jump <- decide(lpD)
699700
if(jump) {
700701
##model$calculate(calcNodesPPomitted)
@@ -941,9 +942,9 @@ sampler_slice <- nimbleFunction(
941942
setAndCalculateTarget = function(value = double()) {
942943
if(discrete) value <- floor(value)
943944
model[[target]] <<- value
944-
lp <- checkLogProb(model$calculate(target))
945+
lp <- checkLogProb(model$calculate(target), target)
945946
if(lp == -Inf) return(-Inf)
946-
lp <- lp + checkLogProb(model$calculate(calcNodesNoSelf))
947+
lp <- lp + checkLogProb(model$calculate(calcNodesNoSelf), target)
947948
returnType(double())
948949
return(lp)
949950
},
@@ -1082,10 +1083,10 @@ sampler_slice_noncentered <- nimbleFunction(
10821083
setAndCalculateTarget = function(value = double()) {
10831084
if(discrete) value <- floor(value)
10841085
model[[target]] <<- value
1085-
lp <- checkLogProb(model$calculate(target))
1086+
lp <- checkLogProb(model$calculate(target), target)
10861087
if(lp == -Inf) return(-Inf)
10871088
lp <- lp + updateNoncentered(value)
1088-
lp <- lp + checkLogProb(model$calculate(calcNodesNoSelf))
1089+
lp <- lp + checkLogProb(model$calculate(calcNodesNoSelf), target)
10891090
returnType(double())
10901091
return(lp)
10911092
},
@@ -1204,6 +1205,7 @@ sampler_ess <- nimbleFunction(
12041205
## checks
12051206
if(length(target) > 1) stop('elliptical slice sampler only applies to one target node')
12061207
if(!(model$getDistribution(target) %in% c('dnorm', 'dmnorm'))) stop('elliptical slice sampler only applies to normal distributions')
1208+
targetNames <- createNamesString(target)
12071209
},
12081210
run = function() {
12091211
u <- model$getLogProb(calcNodesNoSelf) - rexp(1, 1)
@@ -1215,7 +1217,7 @@ sampler_ess <- nimbleFunction(
12151217
theta_min <- theta - 2*Pi
12161218
theta_max <- theta
12171219
values(model, target) <<- f[1:d]*cos(theta) + nu[1:d]*sin(theta) + target_mean[1:d]
1218-
lp <- checkLogProb(model$calculate(calcNodesNoSelf))
1220+
lp <- checkLogProb(model$calculate(calcNodesNoSelf), targetNames)
12191221
numContractions <- 0
12201222
while((is.nan(lp) | lp < u) & theta_max - theta_min > eps & numContractions < maxContractions) { # must be is.nan()
12211223
## The checks for theta_max - theta_min small and max number of contractions are
@@ -1224,7 +1226,7 @@ sampler_ess <- nimbleFunction(
12241226
if(theta < 0) theta_min <- theta else theta_max <- theta
12251227
theta <- runif(1, theta_min, theta_max)
12261228
values(model, target) <<- f[1:d]*cos(theta) + nu[1:d]*sin(theta) + target_mean[1:d]
1227-
lp <- checkLogProb(model$calculate(calcNodesNoSelf))
1229+
lp <- checkLogProb(model$calculate(calcNodesNoSelf), targetNames)
12281230
numContractions <- numContractions + 1
12291231
}
12301232
if(theta_max - theta_min <= eps | numContractions == maxContractions) {
@@ -1300,6 +1302,7 @@ sampler_AF_slice <- nimbleFunction(
13001302
stop('sliceWidths must be a numeric vector')
13011303
if(length(widthVec) != d) stop('sliceWidths must have length = ', d)
13021304
if(adaptFactorInterval < 2) stop('sampler_AF_slice: `adaptFactorInterval` must be at least 2')
1305+
targetNames <- createNamesString(target)
13031306
},
13041307
run = function() {
13051308
maxContractionsReached <- FALSE
@@ -1370,9 +1373,9 @@ sampler_AF_slice <- nimbleFunction(
13701373
for(i in 1:d)
13711374
if(discrete[i] == 1) targetValues[i] <- floor(targetValues[i])
13721375
values(model, target) <<- targetValues
1373-
lp <- checkLogProb(model$calculate(calcNodesProposalStage))
1376+
lp <- checkLogProb(model$calculate(calcNodesProposalStage), targetNames)
13741377
if(lp == -Inf) return(lp)
1375-
lp <- lp + checkLogProb(model$calculate(calcNodesDepStage))
1378+
lp <- lp + checkLogProb(model$calculate(calcNodesDepStage), targetNames)
13761379
returnType(double())
13771380
return(lp)
13781381
},
@@ -1493,13 +1496,14 @@ sampler_crossLevel <- nimbleFunction(
14931496
lowConjugateGetLogDensityFunctions[[iLN]] <- getPosteriorDensityFromConjSampler(lowConjugateSamplerFunctions[[iLN]])
14941497
}
14951498
my_setAndCalculateTop <- setAndCalculate(model, target)
1499+
targetNames <- createNamesString(target)
14961500
},
14971501
run = function() {
14981502
modelLP0 <- model$getLogProb(calcNodes)
14991503
propLP0 <- 0
15001504
for(iSF in seq_along(lowConjugateGetLogDensityFunctions)) { propLP0 <- propLP0 + lowConjugateGetLogDensityFunctions[[iSF]]$run() }
15011505
propValueVector <- topRWblockSamplerFunction$generateProposalVector()
1502-
topLP <- checkLogProb(my_setAndCalculateTop$run(propValueVector))
1506+
topLP <- checkLogProb(my_setAndCalculateTop$run(propValueVector), targetNames)
15031507
if(topLP == -Inf) {
15041508
nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodes, logProb = TRUE)
15051509
}
@@ -1510,7 +1514,7 @@ sampler_crossLevel <- nimbleFunction(
15101514
propLP1 <- 0
15111515
for(iSF in seq_along(lowConjugateGetLogDensityFunctions))
15121516
propLP1 <- propLP1 + lowConjugateGetLogDensityFunctions[[iSF]]$run()
1513-
logMHR <- checkLogProb(modelLP1) - checkLogProb(modelLP0) - checkLogProb(propLP1) + checkLogProb(propLP0)
1517+
logMHR <- checkLogProb(modelLP1, targetNames) - checkLogProb(modelLP0, targetNames) - checkLogProb(propLP1, targetNames) + checkLogProb(propLP0, targetNames)
15141518
jump <- decide(logMHR)
15151519
if(jump) {
15161520
nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
@@ -1710,6 +1714,7 @@ sampler_RW_dirichlet <- nimbleFunction(
17101714
## checks
17111715
if(length(model$expandNodeNames(target)) > 1) stop('RW_dirichlet sampler only applies to one target node')
17121716
if(model$getDistribution(target) != 'ddirch') stop('can only use RW_dirichlet sampler for dirichlet distributions')
1717+
targetNames <- createNamesString(target)
17131718
},
17141719
run = function() {
17151720
if(thetaVec[1] == 0) thetaVec <<- values(model, target) ## initialization
@@ -1722,7 +1727,7 @@ sampler_RW_dirichlet <- nimbleFunction(
17221727
thetaVecProp <- thetaVec
17231728
thetaVecProp[i] <- propValue
17241729
values(model, target) <<- thetaVecProp / sum(thetaVecProp)
1725-
logMHR <- alphaVec[i]*propLogScale + currentValue - propValue + checkLogProb(model$calculateDiff(calcNodesNoSelf))
1730+
logMHR <- alphaVec[i]*propLogScale + currentValue - propValue + checkLogProb(model$calculateDiff(calcNodesNoSelf), targetNames)
17261731
jump <- decide(logMHR)
17271732
} else jump <- FALSE
17281733
if(adaptive & jump) timesAcceptedVec[i] <<- timesAcceptedVec[i] + 1
@@ -1815,6 +1820,7 @@ sampler_RW_wishart <- nimbleFunction(
18151820
if(!all(dim(propCov) == nTheta)) stop('propCov matrix must have dimension ', d, 'x', d)
18161821
if(!isSymmetric(propCov)) stop('propCov matrix must be symmetric')
18171822
if(adaptInterval < 2) stop('sampler_RW_wishart: `adaptInterval` must be at least 2')
1823+
targetNames <- createNamesString(target)
18181824
},
18191825
run = function() {
18201826
currentValue <<- model[[target]]
@@ -1842,7 +1848,7 @@ sampler_RW_wishart <- nimbleFunction(
18421848
## matrix multiply to get proposal value (matrix)
18431849
model[[target]] <<- t(propValue_chol) %*% propValue_chol
18441850
## decide and jump
1845-
logMHR <- checkLogProb(model$calculateDiff(calcNodes))
1851+
logMHR <- checkLogProb(model$calculateDiff(calcNodes), targetNames)
18461852
deltaDiag <- thetaVec_prop[1:d]-thetaVec[1:d]
18471853
for(i in 1:d) logMHR <- logMHR + (d+2-i)*deltaDiag[i] ## took me quite a while to derive this
18481854
jump <- decide(logMHR)
@@ -1940,6 +1946,7 @@ sampler_RW_lkj_corr_cholesky <- nimbleFunction(
19401946
if(dist != 'dlkj_corr_cholesky') stop('RW_lkj_corr_cholesky sampler can only be used with the dlkj_corr_cholesky distribution.')
19411947
if(d < 2) stop('RW_lkj_corr_cholesky sampler requires target node dimension to be at least 2x2.')
19421948
if(adaptFactorExponent < 0) stop('Cannot use RW_lkj_corr_cholesky sampler with adaptFactorExponent control parameter less than 0.')
1949+
targetNames <- createNamesString(target)
19431950
},
19441951
run = function() {
19451952
## calculate transformed values (in unconstrained space) and partial sums in each column
@@ -1965,8 +1972,8 @@ sampler_RW_lkj_corr_cholesky <- nimbleFunction(
19651972
propValue[jprime] <<- z[jprime, i] * sqrt(partialSumsProp[jprime])
19661973
}
19671974
model[[target]][j:i, i] <<- propValue[j:i]
1968-
logMHR <- checkLogProb(calculateDiff(model, calcNodesNoSelf)) +
1969-
checkLogProb(calculateDiff(model, target))
1975+
logMHR <- checkLogProb(calculateDiff(model, calcNodesNoSelf), targetNames) +
1976+
checkLogProb(calculateDiff(model, target), targetNames)
19701977
## Adjust MHR to account for non-symmetric proposal by adjusting prior on U to transformed scale (i.e., y).
19711978
## cosh component is for dz/dy and other component is for du/dz where 'u' is the corr matrix.
19721979
## This follows Stan reference manual Section 10.12 (for version 2.27).
@@ -2102,6 +2109,7 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction(
21022109
if(adaptFactorExponent < 0) stop('Cannot use RW_block_lkj_corr_cholesky sampler with adaptFactorExponent control parameter less than 0.')
21032110
if(adaptInterval < 2) stop('sampler_RW_block_lkj_corr_cholesky: `adaptInterval` must be at least 2')
21042111

2112+
targetNames <- createNamesString(target)
21052113
},
21062114
run = function() {
21072115
transform(model[[target]]) # compute z and partialSums
@@ -2138,12 +2146,12 @@ sampler_RW_block_lkj_corr_cholesky <- nimbleFunction(
21382146
## Adjust for log determinant term from initial values
21392147
logMHR <- logMHR - logDetJac
21402148

2141-
lpD <- checkLogProb(calculateDiff(model, calcNodesProposalStage))
2149+
lpD <- checkLogProb(calculateDiff(model, calcNodesProposalStage), targetNames)
21422150
if(lpD == -Inf) {
21432151
nimCopy(from = mvSaved, to = model, row = 1, nodes = calcNodesProposalStage, logProb = TRUE)
21442152
jump <- FALSE
21452153
} else {
2146-
logMHR <- logMHR + lpD + checkLogProb(calculateDiff(model, calcNodesDepStage))
2154+
logMHR <- logMHR + lpD + checkLogProb(calculateDiff(model, calcNodesDepStage), targetNames)
21472155
jump <- decide(logMHR)
21482156
if(jump) {
21492157
nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)
@@ -2418,7 +2426,7 @@ CAR_scalar_RW <- nimbleFunction(
24182426
propValue <- rnorm(1, mean = model[[targetScalar]], sd = scale)
24192427
model[[targetScalar]] <<- propValue
24202428
lp1 <- dcarList[[1]]$run() + model$calculate(depNodes)
2421-
logMHR <- checkLogProb(lp1) - checkLogProb(lp0)
2429+
logMHR <- checkLogProb(lp1, targetScalar) - checkLogProb(lp0, targetScalar)
24222430
jump <- decide(logMHR)
24232431
if(jump) {
24242432
model$calculate(targetScalar)
@@ -3325,6 +3333,7 @@ sampler_barker <- nimbleFunction(
33253333
## checks
33263334
if(!isTRUE(nimbleOptions('enableDerivs'))) stop('must enable NIMBLE derivatives, set nimbleOptions(enableDerivs = TRUE)', call. = FALSE)
33273335
if(!isTRUE(model$modelDef[['buildDerivs']])) stop('must set buildDerivs = TRUE when building model', call. = FALSE)
3336+
targetNames <- createNamesString(target)
33283337
},
33293338
run = function() {
33303339
current <- my_parameterTransform$transform(values(model, targetNodes))
@@ -3359,8 +3368,8 @@ sampler_barker <- nimbleFunction(
33593368

33603369
gradProposed <<- gradient(proposal)
33613370
newLogDetJacobian <- my_parameterTransform$logDetJacobian(proposal)
3362-
lpD <- checkLogProb(model$calculateDiff(calcNodes)) + checkLogProb(newLogDetJacobian) -
3363-
checkLogProb(oldLogDetJacobian) + checkLogProb(calculateLogHastingsRatio(diff))
3371+
lpD <- checkLogProb(model$calculateDiff(calcNodes), targetNames) + checkLogProb(newLogDetJacobian, targetNames) -
3372+
checkLogProb(oldLogDetJacobian, targetNames) + checkLogProb(calculateLogHastingsRatio(diff), targetNames)
33643373
jump <- decide(lpD)
33653374

33663375
if(jump) nimCopy(from = model, to = mvSaved, row = 1, nodes = calcNodes, logProb = TRUE)

‎packages/nimble/R/MCMC_utils.R‎

Lines changed: 10 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -53,10 +53,11 @@ decideAndJump <- nimbleFunction(
5353
setup = function(model, mvSaved, target, UNUSED) { ## should remove UNUSED argument, after next release of nimbleSMC -DT July 2024
5454
ccList <- mcmc_determineCalcAndCopyNodes(model, target)
5555
copyNodesDeterm <- ccList$copyNodesDeterm; copyNodesStoch <- ccList$copyNodesStoch # not used: calcNodes, calcNodesNoSelf
56+
targetNames <- createNamesString(target)
5657
},
5758
run = function(modelLP1 = double(), modelLP0 = double(), propLP1 = double(), propLP0 = double()) {
5859
## Check each one individually to catch case like `3 - Inf`.
59-
logMHR <- checkLogProb(modelLP1) - checkLogProb(modelLP0) - checkLogProb(propLP1) + checkLogProb(propLP0)
60+
logMHR <- checkLogProb(modelLP1, targetNames) - checkLogProb(modelLP0, targetNames) - checkLogProb(propLP1, targetNames) + checkLogProb(propLP0, targetNames)
6061
jump <- decide(logMHR)
6162
if(jump) {
6263
nimCopy(from = model, to = mvSaved, row = 1, nodes = target, logProb = TRUE)
@@ -72,11 +73,11 @@ decideAndJump <- nimbleFunction(
7273
}
7374
)
7475

75-
checkLogProb <- function(logProb) {
76+
checkLogProb <- function(logProb, target) {
7677
if(is.na(logProb))
7778
return(-Inf)
7879
if(logProb == Inf)
79-
print("MCMC sampling encountered a log probability density value of infinity. Results of sampling may not be valid.")
80+
cat("MCMC sampling of ", target, " encountered a log probability density value of infinity. Results of sampling may not be valid.\n")
8081
return(logProb)
8182
}
8283

@@ -551,6 +552,11 @@ mcmc_checkTargetAD <- function(model, targetNodes, samplerType) {
551552
paste0(dists[!ADok], collapse = ', '))
552553
}
553554

555+
createNamesString <- function(target) {
556+
if(length(target) == 1) return(target)
557+
if(length(target) > 4) target <- c(target[1:4], "...")
558+
return(paste0("{", paste(target, collapse = ', '), "}"))
559+
}
554560

555561
# This is function which builds a new MCMCconf from an old MCMCconf
556562
# This is required to be able to a new C-based MCMC without recompiling
@@ -587,3 +593,4 @@ getBaseClassName <- function(nf) {
587593

588594

589595

596+

‎packages/nimble/inst/CppCode/Utils.cpp‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -63,8 +63,8 @@ bool decide(double lMHr) { // simple function accept or reject based on log Metr
6363
return(false);
6464
}
6565

66-
void checkLogProbWarn() {
67-
_nimble_global_output<<"MCMC sampling encountered a log probability density value of infinity. Results of sampling may not be valid.\n";
66+
void checkLogProbWarn(std::string target) {
67+
_nimble_global_output<<"MCMC sampling of " << target << " encountered a log probability density value of infinity. Results of sampling may not be valid.\n";
6868
nimble_print_to_R(_nimble_global_output);
6969
}
7070

0 commit comments

Comments
 (0)