Skip to content

Commit 79bd822

Browse files
danielturekpaciorek
danielturek
and
paciorek
authored
updates for correct dimension (#1605)
* updates for correct dimension * added a test * now the conj offset and coeff sizes are actually correct --------- Co-authored-by: Christopher Paciorek <paciorek@stat.berkeley.edu>
1 parent 8c0a6b2 commit 79bd822

2 files changed

Lines changed: 75 additions & 23 deletions

File tree

‎packages/nimble/R/MCMC_conjugacy.R‎

Lines changed: 40 additions & 23 deletions
Original file line numberDiff line numberDiff line change
@@ -570,15 +570,27 @@ conjugacyClass <- setRefClass(
570570
}, list(DEP_NAMES = as.name(paste0( 'dep_', distLinkName, '_nodeNames')),
571571
N_DEP = as.name(paste0('N_dep_', distLinkName)),
572572
DEP_CONTROL_NAME = as.name(paste0( 'dep_', distLinkName))))
573+
forSizeDetermination <- character()
574+
## need to include dependent node 'value', if multivariate or if doing dependentScreen:
575+
## Dec 2025
576+
if(distDim > 0 || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
577+
forSizeDetermination <- c(forSizeDetermination, 'value')
578+
}
579+
## if we'll be calculating coeff and offset later, then we also (always) need to extract the sizes
580+
## of the relevant parameter of the dependent node distributions
581+
## (in order to set the correct sizes (and max size) for the coeff and offset terms)
582+
## Dec 2025
583+
if(currentLink %in% c('additive', 'multiplicative', 'multiplicativeScalar', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
584+
forSizeDetermination <- c(forSizeDetermination, dependents[[distName]]$param)
585+
}
573586
## revamp of the code below for size determination,
574587
## to now separate the sizes (and maximum size) of any multivariate
575588
## dependent nodes, and the parameters of dependent nodes, separately.
576589
## this was inititiated by the addition of the gamma-dcar_normal conjugacy,
577590
## where the parameters of dcar_normal often have longer length than the dcar_normal node itself.
578591
## Nov 2025
579-
mvParams <- c('value', neededParams)
580-
mvParams <- mvParams[distDimParams[mvParams] > 0]
581-
for(param in mvParams) {
592+
forSizeDetermination <- c(forSizeDetermination, neededParams[distDimParams[neededParams] > 0])
593+
for(param in forSizeDetermination) {
582594
if(param == 'value') {
583595
## particular call to extract sizes of the *node itself*
584596
functionBody$addCode({
@@ -617,7 +629,7 @@ conjugacyClass <- setRefClass(
617629
## July 2017
618630
targetNdim <- as.numeric(getDimension(prior))
619631
targetCoeffNdim <- switch(as.character(targetNdim), `0`=0, `1`=2, `2`=2, stop())
620-
if(targetCoeffNdim == 2 && link == 'multiplicativeScalar') ## Handles wish/invwish. There are no cases where we allow non-scalar 'coeff'.
632+
if(targetCoeffNdim == 2 && link == 'multiplicativeScalar') ## handles wish/invwish, there are no cases where we allow non-scalar 'coeff'
621633
targetCoeffNdim <- 0
622634
for(iDepCount in seq_along(dependentCounts)) {
623635
distLinkName <- distLinkNameList[[iDepCount]]
@@ -627,13 +639,14 @@ conjugacyClass <- setRefClass(
627639
## the 2's here are *only* to prevent warnings about assigning into member
628640
## variable names using local assignment '<-', so changed the names to "...2"
629641
## so it doesn't recognize the ref class field name
642+
paramSizeMaxName <- as.name(paste0('dep_', distLinkName, '_', dependents[[distName]]$param, '_sizeMax'))
630643
inputList <- list(DEP_OFFSET_VAR2 = as.name(paste0('dep_', distLinkName, '_offset')),
631644
DEP_COEFF_VAR2 = as.name(paste0('dep_', distLinkName, '_coeff')),
632-
DECLARE_SIZE_OFFSET = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_value_sizeMax')), as.name(paste0('dep_', distLinkName, '_value_sizeMax')), targetNdim),
633-
DECLARE_SIZE_COEFF = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), as.name(paste0('dep_', distLinkName, '_value_sizeMax')), quote(d), targetCoeffNdim))
645+
DECLARE_SIZE_OFFSET = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), paramSizeMaxName, paramSizeMaxName, targetNdim ),
646+
DECLARE_SIZE_COEFF = makeDeclareSizeField(as.name(paste0('N_dep_', distLinkName)), paramSizeMaxName, quote(d), targetCoeffNdim))
634647
if(currentLink == 'additive')
635648
functionBody$addCode(
636-
DEP_OFFSET_VAR2 <- array(0, dim = DECLARE_SIZE_OFFSET),
649+
DEP_OFFSET_VAR2 <- array(0, dim = DECLARE_SIZE_OFFSET),
637650
inputList)
638651
if(currentLink %in% c('multiplicative', 'multiplicativeScalar'))
639652
functionBody$addCode(
@@ -881,18 +894,18 @@ conjugacyClass <- setRefClass(
881894
distName <- distNameList[[iDepCount]]
882895
currentLink <- currentLinkList[[iDepCount]]
883896

884-
if(currentLink %in% c('additive', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
897+
if(currentLink %in% c('additive', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))
885898
functionBody$addCode({
886899
for(iDep in 1:N_DEP) {
887-
thisSize <- DEP_SIZES[iDep]
900+
thisSize <- DEP_PARAM_SIZES[iDep]
888901
DEP_OFFSET_VAR[iDep, 1:thisSize] <<- model$getParam(DEP_NAMES[iDep], PARAM_NAME)
889902
}
890903
},
891-
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
892-
DEP_SIZES = as.name(paste0('dep_', distLinkName, '_value_sizes')),
893-
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')),
894-
DEP_NAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
895-
PARAM_NAME = dependents[[distName]]$param))
904+
list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
905+
DEP_PARAM_SIZES = as.name(paste0('dep_', distLinkName, '_', dependents[[distName]]$param, '_sizes')),
906+
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')),
907+
DEP_NAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
908+
PARAM_NAME = dependents[[distName]]$param))
896909
}
897910
}
898911
if(any(allCurrentLinks %in% c('multiplicative', 'linear')) || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
@@ -912,22 +925,22 @@ conjugacyClass <- setRefClass(
912925
currentLink <- currentLinkList[[iDepCount]]
913926

914927
if(currentLink %in% c('multiplicative', 'linear') || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
915-
inputList <- list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
916-
DEP_SIZES = as.name(paste0('dep_', distLinkName, '_value_sizes')),
917-
DEP_COEFF_VAR = as.name(paste0('dep_', distLinkName, '_coeff')),
918-
DEP_NAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
919-
PARAM_NAME = dependents[[distName]]$param,
920-
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')))
928+
inputList <- list(N_DEP = as.name(paste0('N_dep_', distLinkName)),
929+
DEP_PARAM_SIZES = as.name(paste0('dep_', distLinkName, '_', dependents[[distName]]$param, '_sizes')),
930+
DEP_COEFF_VAR = as.name(paste0('dep_', distLinkName, '_coeff')),
931+
DEP_NAMES = as.name(paste0('dep_', distLinkName, '_nodeNames')),
932+
PARAM_NAME = dependents[[distName]]$param,
933+
DEP_OFFSET_VAR = as.name(paste0('dep_', distLinkName, '_offset')))
921934
if(currentLink == 'linear' || (getNimbleOption('allowDynamicIndexing') && doDependentScreen)) {
922935
forLoopBody$addCode(
923936
for(iDep in 1:N_DEP) {
924-
thisSize <- DEP_SIZES[iDep]
937+
thisSize <- DEP_PARAM_SIZES[iDep]
925938
DEP_COEFF_VAR[iDep, 1:thisSize, sizeIndex] <<- model$getParam(DEP_NAMES[iDep], PARAM_NAME) - DEP_OFFSET_VAR[iDep, 1:thisSize]
926939
}, inputList)
927940
} else {
928941
forLoopBody$addCode(
929942
for(iDep in 1:N_DEP) {
930-
thisSize <- DEP_SIZES[iDep]
943+
thisSize <- DEP_PARAM_SIZES[iDep]
931944
DEP_COEFF_VAR[iDep, 1:thisSize, sizeIndex] <<- model$getParam(DEP_NAMES[iDep], PARAM_NAME)
932945
}, inputList)
933946
}
@@ -1042,7 +1055,11 @@ conjugacyClass <- setRefClass(
10421055
}
10431056

10441057
for(p in c('value', depParamsAvailable)) {
1045-
if(distDimParams[[p]] > 0) {
1058+
## need to include dependent node 'value', if multivariate or if doing dependentScreen:
1059+
if( (p == 'value' && (distDim > 0 || (getNimbleOption('allowDynamicIndexing') && doDependentScreen))) ||
1060+
## for other parameters, they need to be multivariate:
1061+
distDimParams[[p]] > 0
1062+
) {
10461063
forLoopBody$addCode(SIZE_NAME <- SIZE_VALUE[iDep],
10471064
list(SIZE_NAME = as.name(paste0(p, '_size')),
10481065
SIZE_VALUE = as.name(paste0('dep_', distLinkName, '_', p, '_sizes'))))

‎packages/nimble/tests/testthat/test-mcmc.R‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3058,6 +3058,41 @@ test_that('asymptotic correct results from conjugate gamma - CAR_normal sampler'
30583058
expect_true(all(abs(summary_RW - summary_conj) < 0.04))
30593059
})
30603060

3061+
test_that('conjugate correct dimensions for dim=1 target nodes', {
3062+
code <- nimbleCode({
3063+
for(k in 1:ncomponents) {
3064+
for(v in 1:Nn) {
3065+
Nprob[Ni[v]:Nf[v], k] ~ ddirch(alpha = Nalpha0[Ni[v]:Nf[v]])
3066+
}
3067+
}
3068+
for(d in 1:npoints) {
3069+
K[d] ~ dcat(prob = W[1:ncomponents])
3070+
for(v in 1:Nn) {
3071+
Ndata[d, v] ~ dcat(prob = Nprob[Ni[v]:Nf[v], K[d]])
3072+
}
3073+
}
3074+
})
3075+
##
3076+
constants <- list(
3077+
ncomponents = 32,
3078+
npoints = 5,
3079+
Nn = 1,
3080+
Ni = c(1, NA_integer_),
3081+
Nf = c(5, NA_integer_),
3082+
Nalpha0 = rep(0.2, 5))
3083+
data <- list(Ndata = matrix(1:5, nrow = 5, ncol = 1))
3084+
inits <- list(W = rep(1/32, 32), K = 1:5, Nprob = matrix(0.2, nrow = 5, ncol = 32))
3085+
##
3086+
Rmodel <- nimbleModel(code, constants, data, inits)
3087+
##
3088+
expect_no_error(conf <- configureMCMC(Rmodel))
3089+
expect_no_error(Rmcmc <- buildMCMC(conf))
3090+
expect_no_error(Cmodel <- compileNimble(Rmodel))
3091+
expect_no_error(Cmcmc <- compileNimble(Rmcmc, project = Rmodel))
3092+
expect_no_error(Rmcmc$run(1))
3093+
expect_no_error(Cmcmc$run(1))
3094+
})
3095+
30613096
test_that("correct handling of changing `thin` value with `reset=FALSE`", {
30623097
code <- nimbleCode({
30633098
y ~ dnorm(mu, sd = sigma)

0 commit comments

Comments
 (0)