Skip to content

Optional Burn in for the models - mtc.run#79

Open
frostwix wants to merge 9 commits intogertvv:masterfrom
frostwix:master
Open

Optional Burn in for the models - mtc.run#79
frostwix wants to merge 9 commits intogertvv:masterfrom
frostwix:master

Conversation

@frostwix
Copy link

@frostwix frostwix commented Aug 2, 2024

We want to make a proposition to update the mtc.run function as it is not supporting a burn-in phase currently.
The burn-in phase is often required by national agencies for all mcmc runs, and even for conjugates.
We make the mtc.run update as minimalistic as posible and backward compatible.
That is why for backward compatiblity the NEW ARGUMENT n.burnin is added as the last one, so mtc.run will works even if somebody used it without argument names mtc.run(value, value) instead of mtc.run(arg1 = value, arg2 = value).
Please give us your thoughts on the update.

Adaptation and burnin are not the same and currently only the former can be set in mtc.run.
Please read more about it here https://stackoverflow.com/a/78807748/5442527 (example where adaptation is skipped as all models are conjugates) and here https://stackoverflow.com/a/38875637/5442527. Summing up it is important to recognize that burn-in and adaptation refer to two different things and are implemented independently in JAGS/rjags. The adaptation will be run only if at least one parameter in the model requires it, whereas burn-in is run always.

Our update fix the #78 issue.

We did not add any unit tests as the current tests structure will demand a lot of copy pasting. Please feel free to propose something.

Our update is backward compatible and results from mtc.run are the same.
The default value for n.burnin is 0.
Proof that the results from mtc.run are the same.

remotes::install_github("gertvv/gemtc/gemtc@master", force = TRUE)

library(gemtc)

set.seed(42)
model <- mtc.model(parkinson)
# By default, the model will have 4 chains - generate a seed for each
seeds <- sample.int(4, n = .Machine$integer.max)
# Apply JAGS RNG settings to each chain
model$inits <- mapply(c, model$inits, list(
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[1]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[2]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[3]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[4])), SIMPLIFY=FALSE)

results <- mtc.run(model, thin=1)

saveRDS(results, "mtc_run_before.RDS")

remotes::install_github("frostwix/gemtc/gemtc@master", force = TRUE)

library(gemtc)

set.seed(42)
model <- mtc.model(parkinson)
# By default, the model will have 4 chains - generate a seed for each
seeds <- sample.int(4, n = .Machine$integer.max)
# Apply JAGS RNG settings to each chain
model$inits <- mapply(c, model$inits, list(
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[1]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[2]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[3]),
  list(.RNG.name="base::Mersenne-Twister", .RNG.seed=seeds[4])), SIMPLIFY=FALSE)

results <- mtc.run(model, thin=1)

saveRDS(results, "mtc_run_after.RDS")

identical(readRDS("mtc_run_before.RDS")$samples, readRDS("mtc_run_after.RDS")$samples)
# TRUE
identical(readRDS("mtc_run_before.RDS")$deviance, readRDS("mtc_run_after.RDS")$deviance)
# TRUE

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants