diff --git a/packages/nimble/R/MCEM_build.R b/packages/nimble/R/MCEM_build.R index acd039d54..336622a98 100644 --- a/packages/nimble/R/MCEM_build.R +++ b/packages/nimble/R/MCEM_build.R @@ -1,6 +1,6 @@ build_MCEM_expectation_noAD <- nimbleFunction( setup = function(model, mvSamples, burnInDefault, paramNodes, - latentNodes, calcNodes, calcNodesOther, + latentNodes, calcNodes, calcNodesOther, paramDeps, useTransform, transformer, lengthParams, lengthParamsTrans, derivsDelta = 0.0001) { @@ -9,19 +9,6 @@ build_MCEM_expectation_noAD <- nimbleFunction( if(lengthParamsTrans == 1) { lastParamsTrans <- c(-Inf, -1) # reset will be called first anyway } - paramDeps <- model$getDependencies(paramNodes, determOnly = TRUE, self=FALSE) - if(length(paramDeps) > 0) { - keep_paramDeps <- logical(length(paramDeps)) - for(i in seq_along(paramDeps)) { - if(any(paramDeps[i] == calcNodes)) keep_paramDeps[i] <- FALSE - else { - nextDeps <- model$getDependencies(paramDeps[i]) - keep_paramDeps[i] <- any(nextDeps %in% calcNodes) - } - } - paramDeps <- paramDeps[keep_paramDeps] - } - allCalcNodes <- model$topologicallySortNodes(c(paramDeps, calcNodes)) useOther <- length(calcNodesOther) > 0 @@ -31,6 +18,7 @@ build_MCEM_expectation_noAD <- nimbleFunction( reset = function() { lastParamsTrans <<- rep(-Inf, lengthParamsTrans) }, + resetDerivs = function() {}, set_burnIn = function(new_burnIn = integer()) { burnIn <<- new_burnIn }, @@ -52,7 +40,7 @@ build_MCEM_expectation_noAD <- nimbleFunction( nimCopy(from = mvSamples, to = model, nodes = latentNodes, row = i) sample_LL <- model$calculate(allCalcNodes) if(is.na(sample_LL) | is.nan(sample_LL) | sample_LL == -Inf | sample_LL == Inf) - stop("Non-finite log-likelihood occurred; the MCEM optimization cannot continue. Please check the state of the compiled model and determine which parameter values are causing the invalid log-likelihood by calling 'calculate' with subsets of the model parameters (e.g., 'name_of_compiled_model$calculate('y[3]')' to see if node 'y[3]' is the cause of the problem). Note that if your model is maximizing over parameters whose bounds are not constant (i.e., depend on other parameters), this is one possible cause of such problems; in that case you might try running the MCEM without bounds, by setting 'forceNoConstraints = TRUE'.") + stop("Non-finite log-likelihood occurred; the MCEM optimization cannot continue. Please check the state of the compiled model (e.g., checking the parameter values and by calling calculate on subsets of the model parameters). Note that if your model is maximizing over parameters whose bounds are not constant (i.e., depend on other parameters), this is one possible cause of such problems; in that case you might try running the MCEM without bounds, by setting the control list element 'forceNoConstraints = TRUE'.") sum_LL <- sum_LL + sample_LL } logLik <<- sum_LL / (nSamples - burnIn) @@ -241,7 +229,7 @@ build_MCEM_expectation_noAD <- nimbleFunction( build_MCEM_expectation <- nimbleFunction( setup = function(model, mvSamples, burnInDefault, paramNodes, - latentNodes, calcNodes, calcNodesOther, + latentNodes, calcNodes, calcNodesOther, paramDeps, useTransform, transformer, lengthParams, lengthParamsTrans) { logLik <- as.numeric(0) @@ -253,30 +241,27 @@ build_MCEM_expectation <- nimbleFunction( } transform_wrt <- as.integer(1:lengthParamsTrans) - paramDeps <- model$getDependencies(paramNodes, determOnly = TRUE, self=FALSE) - if(length(paramDeps) > 0) { - keep_paramDeps <- logical(length(paramDeps)) - for(i in seq_along(paramDeps)) { - if(any(paramDeps[i] == calcNodes)) keep_paramDeps[i] <- FALSE - else { - nextDeps <- model$getDependencies(paramDeps[i]) - keep_paramDeps[i] <- any(nextDeps %in% calcNodes) - } - } - paramDeps <- paramDeps[keep_paramDeps] - } - + # We do not calculate paramDeps separately allCalcNodes <- model$topologicallySortNodes(c(paramDeps, calcNodes)) useOther <- length(calcNodesOther) > 0 burnIn <- burnInDefault + resetDerivs_ <- TRUE + resetDerivsOther_ <- TRUE }, methods = list( reset = function() { lastParamsTrans <<- rep(-Inf, lengthParamsTrans) gradientLik <<- rep(-Inf, lengthParamsTrans) # perhaps unnecessary }, + resetDerivs = function() { + resetDerivs_ <<- TRUE + resetDerivsOther_ <<- TRUE + }, set_burnIn = function(new_burnIn = integer()) { + # burnin here should be provided as burnin/thin from the original MCMC, + # because this determines the number of samples from the mvSamples, + # which is really the original niter/thin. burnIn <<- new_burnIn }, ## @@ -298,7 +283,8 @@ build_MCEM_expectation <- nimbleFunction( values(model, paramNodes) <<- params for(i in (burnIn+1):nSamples){ nimCopy(from = mvSamples, to = model, nodes = latentNodes, row = i) - sample_derivs <- nimDerivs(model$calculate(allCalcNodes), wrt = paramNodes, order = c(0, 1)) + sample_derivs <- nimDerivs(model$calculate(allCalcNodes), wrt = paramNodes, order = c(0, 1), reset=resetDerivs_) + resetDerivs_ <<- FALSE sample_LL <- sample_derivs$value[1] if(is.na(sample_LL) | is.nan(sample_LL) | sample_LL == -Inf | sample_LL == Inf) stop("Non-finite log-likelihood occurred; the MCEM optimization cannot continue. Please check the state of the compiled model (accessible as 'name_of_model$CobjectInterface') and determine which parameter values are causing the invalid log-likelihood by calling 'calculate' with subsets of the model parameters (e.g., 'name_of_model$CobjectInterface$calculate('y[3]')' to see if node 'y[3]' is the cause of the problem). Note that if your model is maximizing over parameters whose bounds are not constant (i.e., depend on other parameters), this is one possible cause of such problems; in that case you might try running the MCEM without bounds, by setting 'forceNoConstraints = TRUE'.") @@ -311,7 +297,8 @@ build_MCEM_expectation <- nimbleFunction( sum_jacobian <- sum_jacobian / (nSamples - burnIn) # now this is mean of jacobian if(useOther) { - other_llh_derivs <- nimDerivs(model$calculate(calcNodesOther), wrt = paramNodes, order = c(0,1)) + other_llh_derivs <- nimDerivs(model$calculate(calcNodesOther), wrt = paramNodes, order = c(0,1), reset=resetDerivsOther_) + resetDerivsOther_ <<- FALSE logLik <<- logLik + other_llh_derivs$value[1] sum_jacobian <- sum_jacobian + other_llh_derivs$jacobian } @@ -611,7 +598,10 @@ R_MCEM_mcse <- nimbleRcall(function(samples = double(1), m = integer()) {}, #' that the initial states of one MCMC will be the last states from the previous #' MCMC, so they will often be good initial values after multiple iterations. Default=500. #' -#' \item \code{thin} Thinning interval for the MCMC in step 1. Default=1. +#' \item \code{thin} Thinning interval for the MCMC in step 1. Default=1. Note that +#' the computational cost of the maximization step depends on the size of the MCMC sample. +#' If chains are highly autocorrelated, thinning should be a good way to reduce the +#' maximization cost while maintaining most of the statistical information in each sample. #' #' \item \code{alpha} Type I error rate for determining when step 2 has moved #' uphill. See above. Default=0.25. @@ -762,7 +752,9 @@ R_MCEM_mcse <- nimbleRcall(function(samples = double(1), m = integer()) {}, #' last sample from \code{findMLE} will be used. If MLE convergence was #' reasonable, this sample can be used. However, if the last MCEM step made #' a big move in parameter space (e.g. if convergence was not achieved), -#' the last MCMC sample may not be accurate for obtaining \code{vcov}. Default=FALSE. +#' the last MCMC sample may not be accurate for obtaining \code{vcov}. Note that +#' \code{thin} and \code{burnin} will be used from the most recent call to `findMLE`, +#' or their defaults set in the `control` list provided to `buildMCEM`. Default=FALSE. #' #' \item \code{atMLE}. Logical indicating whether you believe the #' \code{params} represents the MLE. If TRUE, one part of the computation @@ -804,6 +796,11 @@ R_MCEM_mcse <- nimbleRcall(function(samples = double(1), m = integer()) {}, #' provided in the call to \code{buildMCEM}. The user does not normally need to #' call this. #' +#' \item \code{getParamNodes}. Return a vector of the parameter node +#' names. This facilitates being sure that the numeric vector of the MLE +#' parameters can be properly interpreted. If there is only one parameter, +#' an extra "_EXTRA_" will appear in the returned vector and should be ignored. +#' #' } #' #' @author Perry de Valpine, Clifford Anderson-Bergman and Nicholas Michaud @@ -858,7 +855,9 @@ buildMCEM <- nimbleFunction( "buffer", "alpha", "beta", "gamma", "C", "numReps", "forceNoConstraints", "verbose") %in% names(dotsList))) stop("From the arguments provided, it looks like you are trying to use the old version of buildMCEM.", - " buildMCEM was rewritten for nimble 1.2.0.") + " buildMCEM was rewritten for nimble 1.2.0. Some arguments were moved to the control list. See help(buildMCEM).") + if(!requireNamespace('mcmcse')) stop("MCEM requires that package `mcmcse` be installed.") + # Extract control elements or set defaults initMuse <- initMdefault <- extractControlElement(control, 'initM', 1000) MfactorUse <- MfactorDefault <- extractControlElement(control, 'Mfactor', 1/3) @@ -1069,19 +1068,29 @@ buildMCEM <- nimbleFunction( thin = thinDefault, control = mcmcControl, print = FALSE) } + warnOption <- getNimbleOption('MCMCwarnUnsampledStochasticNodes') + nimbleOptions('MCMCwarnUnsampledStochasticNodes' = FALSE) + on.exit(nimbleOptions('MCMCwarnUnsampledStochasticNodes' = warnOption)) mcmc_Latent <- buildMCMC(mcmc_Latent_Conf) mvSamples <- mcmc_Latent$mvSamples setupOutputs(mvSamples) nimbleOptions(verbose = nimbleVerbose) + ## New approach to setting up paramDeps borrowed from Laplace/AGHQ + paramDeps <- model$getDependencies(paramNodes, determOnly = TRUE, self=FALSE) + if(length(paramDeps)) { + calcNodesParents <- model$getParents(calcNodes, determOnly = TRUE) + paramDeps <- paramDeps[!paramDeps %in% calcNodes & paramDeps %in% calcNodesParents] + } + if(useDerivs) { E <- build_MCEM_expectation(model, mvSamples, burnInDefault, paramNodes, - latentNodes, calcNodes, calcNodesOther, + latentNodes, calcNodes, calcNodesOther, paramDeps, useTransform, transformer, lengthParams, lengthParamsTrans) } else { E <- build_MCEM_expectation_noAD(model, mvSamples, burnInDefault, paramNodes, - latentNodes, calcNodes, calcNodesOther, + latentNodes, calcNodes, calcNodesOther, paramDeps, useTransform, transformer, lengthParams, lengthParamsTrans, derivsDeltaUse) @@ -1093,8 +1102,14 @@ buildMCEM <- nimbleFunction( if(length(lastParamsTrans)==1) lastParamsTrans <- c(lastParamsTrans, -1) finalM <- 0L lastDiffQ <- 0 + paramNodesForUser <- paramNodesScalar + if(length(paramNodesForUser)==1) paramNodesForUser <- c(paramNodesForUser, "_EXTRA_") }, methods = list( + getParamNodes = function() { + return(paramNodesForUser) + returnType(character(1)) + }, fix_one_vec = function(x = double(1)) { if(length(x) == 2) { if(x[2] == -1) { @@ -1144,7 +1159,8 @@ buildMCEM <- nimbleFunction( if(M == -1) M <- finalM lastParamsTrans <<- paramsTrans values(model, paramNodes) <<- paramsActual - doMCMC(M, TRUE) + model$calculate(paramDeps) + doMCMC(M, thin = thinUse, reset = TRUE) } FI <- E$fisherInfrmtn(params, trans, returnTrans, atMLE) vc <- inverse(FI) @@ -1152,7 +1168,11 @@ buildMCEM <- nimbleFunction( returnType(double(2)) }, doMCMC = function(M = integer(), thin = integer(), reset = logical(0, default = TRUE)) { - mcmc_Latent$run(M, thin = thin, reset = reset, progressBar = MCMCprogressBarUse) + # MCMC$run errors out if reset=FALSE and thin is not the default of -1 + if(reset) + mcmc_Latent$run(M, thin = thin, reset = reset, progressBar = MCMCprogressBarUse) + else + mcmc_Latent$run(M, thin = -1, reset = reset, progressBar = MCMCprogressBarUse) }, transform = function(params = double(1)) { if(useTransform) @@ -1221,7 +1241,7 @@ buildMCEM <- nimbleFunction( } ## In case parameter nodes are not properly initialized if(any_na(pStart) | any_nan(pStart) | any(abs(pStart)==Inf)) - stop(" [Warning] For findMLE, pStart has some invalid values.") + stop(" [Warning] For findMLE, pStart has some invalid values. If you did not provide the pStart argument, current values in the model were used.") if(!continue) resetControls() else @@ -1247,8 +1267,20 @@ buildMCEM <- nimbleFunction( if(burnInUse >= initMuse) stop('mcem quitting: burnIn > initM value') - E$set_burnIn(burnInUse) - if(!continue) doMCMC(0, TRUE) # To get valid initial values + + # To be consistent with runMCMC: + # samples returned are floor((M-burnin)/thin) + # With automated increases in M, the floor may really matter + # The "E" object looks at the samples to get the post-thinning M, without burnin applied + # And we then use a burnin there that is effectively post-thinning. + # To get that right, we calculate here the burnin_post_thin as + # the number of recorded samples - the number of samples consistent with runMCMC use of burnin + # Actually I think it is redundant to do these right here as they are done below too. + M_post_thin <- floor(initMuse / thinUse) + burnin_post_thin <- M_post_thin - floor((initMuse-burnInUse) / thinUse) + E$set_burnIn(burnin_post_thin) + # no iterations, so thin doesn't matter here + if(!continue) doMCMC(0, thin = 1, reset = TRUE) # To get valid initial values params <- pStart if(!continue) { @@ -1270,8 +1302,9 @@ buildMCEM <- nimbleFunction( } } } - values(model, paramNodes) <<- params } + values(model, paramNodes) <<- params + model$calculate(paramDeps) } m <- initMuse finalM <<- m @@ -1288,6 +1321,7 @@ buildMCEM <- nimbleFunction( zGamma <<- qnorm(gammaUse, 0, 1, lower.tail=FALSE) hitMaxM <- FALSE paramsTrans <- transform(params) + E$resetDerivs() while(((numberOfTimesFlat < Cuse) & (itNum < maxIterUse)) | (itNum < minIterUse)) { # || itNum <= 1) { #endCrit > C){ # We forced at least 2 iterations. Why? checkInterrupt() @@ -1301,13 +1335,16 @@ buildMCEM <- nimbleFunction( hitMaxM <- TRUE } lastParamsTrans <<- paramsTrans # records params used for the latent state MCMC sample - doMCMC(m, thin = thin, reset = TRUE) # initial mcmc run of size m + M_post_thin <- floor(m / thinUse) + burnin_post_thin <- M_post_thin - floor((m-burnInUse) / thinUse) + E$set_burnIn(burnin_post_thin) + # First time through: params were set above + # Later times through: params were set below: + # at E$llh_sample(paramTrans) + # In the event of ascent with increasing M, we only get back here after a uphill move + # is determined, in which case E$llh_sample(paramTrans) will have been last use of model. + doMCMC(m, thin = thinUse, reset = TRUE) # initial mcmc run of size m # paramsPrev <- params #store previous params value - paramsTransPrev <- paramsTrans - ## if(itNum == 1) { - ## sample_logLiks_new <- E$llh_sample(paramsTransPrev) - ## } - ## sample_logLiks_prev <- sample_logLiks_new while(!acceptThisIter){ E$reset() if(useDerivs) { @@ -1330,12 +1367,14 @@ buildMCEM <- nimbleFunction( } } paramsTrans = optimOutput$par - # In the future, we could revise this to append llh_sample values when M is being extended. + # In the future, we could revise this to calculate only from appended samples # For now, we simply recalculate them all. This should be minor compared to cost of optim. - sample_logLiks_prev <- E$llh_sample(paramsTransPrev) - sample_logLiks_new <- E$llh_sample(paramsTrans) # Side effect: this puts the latest params into the model + # We compare to the parameters of the MCMC sampling. + sample_logLiks_prev <- E$llh_sample(lastParamsTrans) + ## PARAMETERS ARE UPDATED IN THE MODEL HERE: + sample_logLiks_new <- E$llh_sample(paramsTrans) # Side effect: this puts the latest params into the model and calculates # make nimbleRcall or other solution. - sdDeltaQ <- R_MCEM_mcse(sample_logLiks_new - sample_logLiks_prev, m-burnInUse) + sdDeltaQ <- R_MCEM_mcse(sample_logLiks_new - sample_logLiks_prev, m-burnInUse) # second argument is actually not used varDeltaQ <- sdDeltaQ*sdDeltaQ ## varDeltaQ <- cvarCalc$run(m, params, paramsPrev) ## sdDeltaQ <- sqrt(varDeltaQ) #asymptotic std. error @@ -1358,7 +1397,17 @@ buildMCEM <- nimbleFunction( hitMaxM <- TRUE } m <- m + mAdd - doMCMC(mAdd, thin = thin, reset = FALSE) + M_post_thin <- floor(m / thinUse) + burnin_post_thin <- M_post_thin - floor((m-burnInUse) / thinUse) + E$set_burnIn(burnin_post_thin) + # We need to revert to parameters used during sampling to continue sampling with the same parameters + tempParams <- inverseTransform(lastParamsTrans) # use same params as from doMCMC above + values(model, paramNodes) <<- tempParams + model$calculate(paramDeps) + + doMCMC(mAdd, thin = thinUse, reset = FALSE) + # Now we are guaranteed to back to the optim step, so we do not need to + # reset the model parameters right here (post-MCMC) } else { acceptThisIter <- TRUE stronglyFlat <- diffQ + zGamma*sdDeltaQ <= tolUse # We reject H0:DeltaQ==tol (one-sided to H0:DeltaQ < tol) with (low) Type I error rate gamma. i.e. we believe true diffQ < tol, so we did not move uphill more than tol, so we stayed flat. Using low Type I error rate means we have good evidence that true diffQ < tol. @@ -1381,10 +1430,10 @@ buildMCEM <- nimbleFunction( cat(" [Note] Convergence Criterion: ", diffQ+zGamma*sdDeltaQ, ".\n", sep = "") } } - if(itNum == maxIterUse) - cat(" [Note] Stopping MCEM because maximum number of MCEM iterations\n", - "(maxIter=", maxIterUse, ") was reached\n.") } + if(itNum == maxIterUse) + cat(" [Note] Stopping MCEM because maximum number of MCEM iterations\n", + "(maxIter=", maxIterUse, ") was reached.\n") } finalM <<- m if(!returnTrans) { diff --git a/packages/nimble/R/MCMC_configuration.R b/packages/nimble/R/MCMC_configuration.R index a627b26bc..3a91dc2d6 100644 --- a/packages/nimble/R/MCMC_configuration.R +++ b/packages/nimble/R/MCMC_configuration.R @@ -1,4 +1,3 @@ - samplerConf <- setRefClass( Class = 'samplerConf', fields = list( @@ -115,9 +114,9 @@ derivedConf <- setRefClass( #' conf$printMonitors() #' conf$printSamplers() MCMCconf <- setRefClass( - - Class = 'MCMCconf', - + + Class = 'MCMCconf', + fields = list( model = 'ANY', monitors = 'ANY', @@ -137,9 +136,9 @@ MCMCconf <- setRefClass( mvSamples1Conf = 'ANY', mvSamples2Conf = 'ANY' ), - + methods = list( - + initialize = function(model, nodes, control = list(), ##rules, monitors, thin = 1, monitors2 = character(), thin2 = 1, @@ -199,7 +198,7 @@ enableWAIC: A logical argument, specifying whether to enable WAIC calculations f controlWAIC A named list of inputs that control the behavior of the WAIC calculation, passed as the \'control\' input to \'buildWAIC\'. See \'help(waic)\`. -print: A logical argument specifying whether to print the montiors and samplers. Default is TRUE. +print: A logical argument specifying whether to print the monitors and samplers. Default is TRUE. ...: Additional named control list elements for default samplers, or additional arguments to be passed to the autoBlock function when autoBlock = TRUE. ' @@ -242,7 +241,7 @@ print: A logical argument specifying whether to print the montiors and samplers. } } } - + if(missing(nodes)) { nodes <- model$getNodeNames(stochOnly = TRUE, includeData = FALSE, includePredictive = samplePredictiveNodes) # Check of all(model$isStoch(nodes)) is not needed in this case @@ -250,21 +249,21 @@ print: A logical argument specifying whether to print the montiors and samplers. nodes <- character(0) } else nodes <- filterOutDataNodes(nodes) ## configureMCMC *never* assigns samplers to data nodes allowData_global <<- FALSE - + addDefaultSampler(nodes = nodes, useConjugacy = useConjugacy, onlyRW = onlyRW, onlySlice = onlySlice, multivariateNodesAsScalars = multivariateNodesAsScalars, print = FALSE) - + if(isTRUE (mean )) addDerivedQuantity('mean' , control = list(nodes = monitors)) if(is.character(mean )) addDerivedQuantity('mean' , control = list(nodes = mean )) if(isTRUE (variance)) addDerivedQuantity('variance', control = list(nodes = monitors)) if(is.character(variance)) addDerivedQuantity('variance', control = list(nodes = variance)) if(isTRUE (logProb )) addDerivedQuantity('logProb' , control = list(nodes = '.all' )) if(!is.logical (logProb )) addDerivedQuantity('logProb' , control = list(nodes = logProb )) - + if(print) show() ##printSamplers() }, @@ -280,7 +279,7 @@ print: A logical argument specifying whether to print the montiors and samplers. For internal use. Adds default MCMC samplers to the specified nodes. ' useNewConfigureMCMC <- isTRUE(getNimbleOption("useNewConfigureMCMC")) - + controlDefaultsArg <- list(...) for(i in seq_along(control)) controlDefaultsArg[[names(control)[i]]] <- control[[i]] @@ -349,11 +348,11 @@ For internal use. Adds default MCMC samplers to the specified nodes. nodes <- model$modelDef$maps$graphID_2_nodeName[nodeIDs] nodesForPosteriorPredictiveSampler <- model$modelDef$maps$graphID_2_nodeName[predictiveNodeIDsToSample] } - + if(useConjugacy) conjugacyResultsAll <- nimble:::conjugacyRelationshipsObject$checkConjugacy(model, nodeIDs) ## Later, this can go through model$checkConjugacy if we make it check whether nodes are already nodeIDs. To isolate changes, I am doing it directly here. nodeDeclIDs <- model$modelDef$maps$graphID_2_declID[nodeIDs] ## Below, nodeDeclIDs[i] gives the nodeDeclID. We could add an interface to get this. nodeDeclID_2_nodes <- split(nodes, nodeDeclIDs) - + uniqueNodeDeclIDs <- unique(nodeDeclIDs) nodeTraits <- lapply(uniqueNodeDeclIDs, function(x) { @@ -361,7 +360,7 @@ For internal use. Adds default MCMC samplers to the specified nodes. dist <- declInfo$distributionName distInfo <- getDistributionInfo(dist) discrete <- distInfo$discrete - ## Following can be replaced by an efficiency version model$isBinary + ## Following can be replaced by an efficiency version model$isBinary binary <- dist == 'dbern' ## If dist == 'dbin', then binary-ness will be checked for each node, below ## This could be improved to see if they all have a literal "1", for example. @@ -369,7 +368,7 @@ For internal use. Adds default MCMC samplers to the specified nodes. ## For nodeScalarComponents, we will check a single node ## We could use returnType = 'ids', but we have a warning generated in that case, ## for future investigation. - + ## Determining nodeLength is a bit tricky. ## The only purpose is to determine scalar vs. non-scalar. ## In the future, we may want to make this available in distributionInfo. @@ -393,11 +392,11 @@ For internal use. Adds default MCMC samplers to the specified nodes. } ) names(nodeTraits) <- as.character(uniqueNodeDeclIDs) - + allDists <- unlist(lapply(model$modelDef$declInfo, `[[`, 'distributionName')) allDists <- allDists[!is.na(allDists)] check_dCRP <- any(allDists == "dCRP") - + clusterNodeInfo <- NULL; dcrpNode <- NULL; numCRPnodes <- 0; clusterNodeParams <- NULL for(i in seq_along(nodes)) { @@ -429,7 +428,7 @@ For internal use. Adds default MCMC samplers to the specified nodes. ## if node is the root of a posterior predictive (entirely non-data) network of nodes, assign 'posterior_predictive' sampler if(node %in% nodesForPosteriorPredictiveSampler) { addSampler(target = node, type = 'posterior_predictive', control = controlDefaultsArg); next } - + ## for multivariate nodes, either add a conjugate sampler, RW_multinomial, or RW_block sampler if(nodeLength > 1) { if(useConjugacy) { @@ -472,10 +471,10 @@ For internal use. Adds default MCMC samplers to the specified nodes. addSampler(target = node, type = 'barker', silent = TRUE, control = controlDefaultsArg); next } else { addSampler(target = node, type = 'RW_block', silent = TRUE, control = controlDefaultsArg); next } } - + if(onlyRW && !discrete) { addSampler(target = node, type = 'RW', control = controlDefaultsArg); next } if(onlySlice) { addSampler(target = node, type = 'slice', control = controlDefaultsArg); next } - + ## if node passes checkConjugacy(), assign 'conjugate_dxxx' sampler if(useConjugacy) { conjugacyResult <- conjugacyResultsAll[[node]] @@ -483,16 +482,16 @@ For internal use. Adds default MCMC samplers to the specified nodes. addConjugateSampler(conjugacyResult = conjugacyResult, dynamicallyIndexed = model$modelDef$varInfo[[model$getVarNames(nodes=node)]]$anyDynamicallyIndexed); next } } - + ## if node is discrete 0/1 (binary), assign 'binary' sampler if(binary) { addSampler(target = node, type = 'binary', control = controlDefaultsArg); next } - + ## for categorical nodes, assign a 'categorical' sampler if(nodeDist == 'dcat') { addSampler(target = node, type = 'categorical', control = controlDefaultsArg); next } - + ## if node distribution is discrete, assign 'slice' sampler if(discrete) { addSampler(target = node, type = 'slice', control = controlDefaultsArg); next } - + ## if node distribution is dgamma and its dependency is dCRP, assign 'CRP_concentration' sampler if(check_dCRP) { if(nodeDist == 'dgamma'){ @@ -506,7 +505,7 @@ For internal use. Adds default MCMC samplers to the specified nodes. } } } - + ## default: 'RW' sampler addSampler(target = node, type = 'RW', control = controlDefaultsArg); next } @@ -517,7 +516,7 @@ For internal use. Adds default MCMC samplers to the specified nodes. ## If anything contraindicates wrapping, we avoid it. E.g., a dangerous case is if a single hyperparameter is involved ## in multiple sets of cluster parameters, but another hyperparameter is involved in one of those sets. ## In that case, the processing of second hyperparameter could turn on the wrapping for some cluster parameters - ## which would mess up sampling of the first hyperparameter. + ## which would mess up sampling of the first hyperparameter. wrap <- TRUE if(!is.null(clusterNodeInfo)) { allClusterNodes <- lapply(clusterNodeInfo, function(x) x$clusterNodes) @@ -536,7 +535,7 @@ For internal use. Adds default MCMC samplers to the specified nodes. if(!all(clusterNodeDeps %in% model$getDependencies(dcrpNode[[k]], stochOnly = TRUE, self = FALSE))) wrap <- FALSE - + ## For now avoid wrapper if any overlap of clusterNodes, as hard to determine if cluster is occupied. ## We'll need to come back to this to handle the mu[xi[i],eta[j]] case if we want to ## avoid sampling empty clusters in that case. @@ -589,7 +588,7 @@ For internal use. Adds default MCMC samplers to the specified nodes. } } } - + setUnsampledNodes() if(print) printSamplers(byType = TRUE) ##show() ##printSamplers() }, @@ -616,7 +615,7 @@ For internal use. Adds default MCMC samplers to the specified nodes. nameToPrint <- gsub('^sampler_', '', conjSamplerName) addSampler(target = conjugacyResult$target, type = conjSamplerFunction, control = conjugacyResult$control, print = print, name = nameToPrint) }, - + addSampler = function(target = character(), ## target argument is *not* expanded (unless targetByNode = TRUE) type = 'RW', control = list(), @@ -698,7 +697,7 @@ The second usage of \'multivariateNodesAsScalars\' occurs when \'default\' is TR if(!missing(expandTarget)) { messageIfVerbose(' [Warning] `expandTarget` argument name has been deprecated; please use `targetByNode.'); targetByNode <- expandTarget } if(!missing(scalarComponents)) { messageIfVerbose(' [Warning] `scalarComponents` argument name has been deprecated; please use `multivariateNodesAsScalars.'); multivariateNodesAsScalars <- scalarComponents } - + nameProvided <- !missing(name) if(is.character(type)) { if(type == 'conjugate') { @@ -764,7 +763,7 @@ The second usage of \'multivariateNodesAsScalars\' occurs when \'default\' is TR ##requiredControlNames <- controlNamesLibrary[[libraryTag]] controlArgs <- c(control, list(...)) thisControlList <- mcmc_generateControlListArgument(control=controlArgs, controlDefaults=controlDefaults) ## should name arguments - + if(!targetByNode) { ## when targetByNode = FALSE, ## no node expansion takes place @@ -778,7 +777,7 @@ The second usage of \'multivariateNodesAsScalars\' occurs when \'default\' is TR addOneSampler(thisSamplerName, samplerFunction, targetExpanded[i], thisControlList, allowData, print) } } - + return(invisible(samplerConfs)) }, @@ -810,7 +809,7 @@ For internal use only nodes <- model$expandNodeNames(nodes) return(nodes[!model$isData(nodes)]) }, - + removeSamplers = function(..., ind, print = FALSE) { ' Removes one or more samplers from an MCMCconf object. @@ -844,7 +843,7 @@ Alias for removeSamplers method ' removeSamplers(...) }, - + ## Note: function prototype is identical to addSampler replaceSamplers = function(...) { ' @@ -877,7 +876,7 @@ Alias for replaceSamplers method ' replaceSamplers(...) }, - + setSamplers = function(..., ind, print = FALSE) { ' Sets the ordering of the list of MCMC samplers. @@ -896,7 +895,7 @@ Alternatively, a character vector may be used to specify a set of model nodes an As another alternative, a list of samplerConf objects may be used as the argument, in which case this ordered list of samplerConf objects will define the samplers in this MCMC configuration object, completely over-writing the current list of samplers. No checking is done to ensure the validity of the contents of these samplerConf objects; only that all elements of the list argument are, in fact, samplerConf objects. print: A logical argument specifying whether to print the new list of samplers (default FALSE). -' +' if(missing(ind)) { ind <- list(...) ind <- unname(unlist(ind)) @@ -921,7 +920,7 @@ Alias for setSamplers method ' setSamplers(...) }, - + printSamplers = function(..., ind, type, displayControlDefaults = FALSE, displayNonScalars = FALSE, displayConjugateDependencies = FALSE, executionOrder = FALSE, byType = FALSE) { ' Prints details of the MCMC samplers. @@ -1185,11 +1184,11 @@ Invisibly returns a list of the current derived quantity function configurations if(length(interval) != 1) stop('derived quantity interval should be a single number') if(interval < 0) stop('derived quantity interval must be at least 1') ## interval = 0 corresponds to matching thin interval if(floor(interval) != interval) stop('derived quantity interval must be an integer') - + thisControlList <- c(control, list(...)) - + addOneDerivedQuantity(thisDerivedName, derivedFunction, interval, thisControlList, print) - + return(invisible(derivedConfs)) }, @@ -1327,7 +1326,7 @@ Details: See the initialize() function ' - + if(isMvSamplesReady(ind)){ messageIfVerbose(' [Note] Changing monitors, even though an MCMC has been built already. When compiling the MCMC, use resetFunctions = TRUE option.') if(ind == 1) @@ -1406,7 +1405,7 @@ See the initialize() function ' setMonitors(..., ind = 2, print = print) }, - + resetMonitors = function() { ' Resets the current monitors and monitors2 lists to nothing. @@ -1417,17 +1416,17 @@ See the initialize() function ' monitors <<- character() monitors2 <<- character() - + if(isMvSamplesReady(1) || isMvSamplesReady(2)){ message('Changing monitors, even though an MCMC has been built already. When compiling the MCMC, use resetFunctions = TRUE option.') mvSamples1Conf <<- NULL mvSamples2Conf <<- NULL } - + return(invisible(NULL)) }, - + printMonitors = function() { ' Prints all current monitors and monitors2 @@ -1501,7 +1500,7 @@ See the initialize() function }, getMvSamplesConf = function(ind = 1){ - + if(isMvSamplesReady(ind) == TRUE) { if(ind == 1) return(mvSamples1Conf) return(mvSamples2Conf) @@ -1515,13 +1514,13 @@ See the initialize() function return(output) } }, - + isMvSamplesReady = function(ind){ if(ind == 1) return(is(mvSamples1Conf, 'function')) #Probably really want to give mvConfs there own class... if(ind == 2) return(is(mvSamples2Conf, 'function')) stop('invalid indicator for isMvSsamplesReady') }, - + makeMvSamplesConf = function(ind){ modelSymbolObjects = model$getSymbolTable()$getSymbolObjects() if(ind == 1) monitorNames = monitors @@ -1529,9 +1528,9 @@ See the initialize() function if(!all(monitorNames %in% names(modelSymbolObjects))) stop('some monitor names are not in the model symbol table; this should never occur') thisModelValuesConf = modelValuesConf(symbolTable(symbols = modelSymbolObjects[monitorNames])) if(ind == 1) mvSamples1Conf <<- thisModelValuesConf - if(ind == 2) mvSamples2Conf <<- thisModelValuesConf + if(ind == 2) mvSamples2Conf <<- thisModelValuesConf }, - + setUnsampledNodes = function() { samplerTargetNodes <- model$expandNodeNames(unlist(lapply(samplerConfs, `[[`, 'target'))) additionalNodesBeingSampled <- character() @@ -1550,12 +1549,12 @@ See the initialize() function allNodesBeingSampled <- unique(c(samplerTargetNodes, additionalNodesBeingSampled)) unsampledNodes <<- setdiff(model$getNodeNames(stochOnly = TRUE, includeData = FALSE), allNodesBeingSampled) }, - + getUnsampledNodes = function() { setUnsampledNodes() return(unsampledNodes) }, - + warnUnsampledNodes = function(includeConfGetUnsampledNodes = TRUE) { if(length(unsampledNodes)) { numUnsampled <- length(unsampledNodes) @@ -1635,10 +1634,10 @@ See the initialize() function #'@param oldConf An optional MCMCconf object to modify rather than creating a new MCMCconf from scratch #'@param ... Additional named control list elements for default samplers, or additional arguments to be passed to the \code{\link{autoBlock}} function when \code{autoBlock = TRUE} #'@author Daniel Turek -#'@export +#'@export #'@details See \code{\link{MCMCconf}} for details on how to manipulate the \code{MCMCconf} object #'@seealso \code{\link{buildMCMC}} \code{\link{runMCMC}} \code{\link{nimbleMCMC}} -configureMCMC <- function(model, nodes, control = list(), +configureMCMC <- function(model, nodes, control = list(), monitors, thin = 1, monitors2 = character(), thin2 = 1, useConjugacy = getNimbleOption('MCMCuseConjugacy'), onlyRW = FALSE, onlySlice = FALSE, @@ -1651,16 +1650,16 @@ configureMCMC <- function(model, nodes, control = list(), ## samplerAssignmentRules system deprecated Nov 2020 -DT ##rules = getNimbleOption('MCMCdefaultSamplerAssignmentRules'), ...) { - + ## samplerAssignmentRules system deprecated Nov 2020 -DT ##if(!inherits(rules, 'samplerAssignmentRules')) stop('rules argument must be a samplerAssignmentRules object') if(!missing(oldConf)){ if(!is(oldConf, 'MCMCconf')) stop('oldConf must be an MCMCconf object, as built by the configureMCMC function') - return(makeNewConfFromOldConf(oldConf)) + return(makeNewConfFromOldConf(oldConf)) } - + if(missing(model)) stop('Either oldConf or model must be supplied') if(missing(monitors)) monitors <- NULL @@ -1677,6 +1676,3 @@ configureMCMC <- function(model, nodes, control = list(), print = print, ...) return(invisible(thisConf)) } - - - diff --git a/packages/nimble/tests/testthat/mcemTestLog_Correct.Rout b/packages/nimble/tests/testthat/mcemTestLog_Correct.Rout index 11913499d..d88d781fb 100644 --- a/packages/nimble/tests/testthat/mcemTestLog_Correct.Rout +++ b/packages/nimble/tests/testthat/mcemTestLog_Correct.Rout @@ -1,3 +1,5 @@ +MCEM error-trapping +MCEM pump with all defaults (roxygen example) [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -53,9 +55,10 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 8. [Note] Current number of MCMC iterations: 9400. [Note] Parameter Estimates: -0.823025 - 1.26278 - [Note] Convergence Criterion: 0.000165076. +0.822968 + 1.26253 + [Note] Convergence Criterion: 0.000157632. +MCEM with pump transform=FALSE [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -74,9 +77,9 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 2188. [Note] Parameter Estimates: -0.818181 - 1.25068 - [Note] Convergence Criterion: 0.0049061. +0.817874 + 1.24692 + [Note] Convergence Criterion: 0.0041643. [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 500. [Note] Parameter Estimates: @@ -84,21 +87,22 @@ Monte Carlo error too big: increasing MCMC sample size. 1.12139 [Note] Convergence Criterion: 0.662289. [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -. [Note] Iteration Number: 1. +(maxIter=1) was reached. + [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: 0.808784 1.20671 [Note] Convergence Criterion: 0.0356889. Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 2. - [Note] Current number of MCMC iterations: 1546. + [Note] Current number of MCMC iterations: 1962. [Note] Parameter Estimates: -0.817103 - 1.2491 - [Note] Convergence Criterion: 0.00535501. +0.817835 + 1.24893 + [Note] Convergence Criterion: 0.00483164. [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -106,67 +110,49 @@ Monte Carlo error too big: increasing MCMC sample size. 1.12136 [Note] Convergence Criterion: 0.670036. [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -. [Note] Iteration Number: 1. +(maxIter=1) was reached. + [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: 0.821669 1.2273 [Note] Convergence Criterion: 0.0335632. [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -.Monte Carlo error too big: increasing MCMC sample size. - [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -.Monte Carlo error too big: increasing MCMC sample size. - [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -.Monte Carlo error too big: increasing MCMC sample size. - [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -. [Note] Iteration Number: 1. +(maxIter=1) was reached. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 2663. [Note] Parameter Estimates: -0.819348 - 1.25075 - [Note] Convergence Criterion: 0.00341703. - [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -.Monte Carlo error too big: increasing MCMC sample size. - [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -.Monte Carlo error too big: increasing MCMC sample size. +0.819054 + 1.24727 + [Note] Convergence Criterion: 0.00284864. [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -.Monte Carlo error too big: increasing MCMC sample size. - [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -. [Note] Iteration Number: 1. +(maxIter=1) was reached. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 8277. [Note] Parameter Estimates: -0.825168 - 1.27138 - [Note] Convergence Criterion: 0.000994448. - [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -.Monte Carlo error too big: increasing MCMC sample size. - [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -.Monte Carlo error too big: increasing MCMC sample size. +0.82409 +1.26675 + [Note] Convergence Criterion: 0.000942662. [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -.Monte Carlo error too big: increasing MCMC sample size. - [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -. [Note] Iteration Number: 1. +(maxIter=1) was reached. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 20000. [Note] Parameter Estimates: -0.824445 - 1.2667 - [Note] Convergence Criterion: 0.000112564. +0.824414 + 1.26661 + [Note] Convergence Criterion: 1.3166e-05. [Note] Stopping MCEM because maximum number of MCEM iterations -(maxIter=1) was reached -. [Note] Iteration Number: 1. +(maxIter=1) was reached. + [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: 0.810924 @@ -226,15 +212,15 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: -0.819486 - 1.24728 - [Note] Convergence Criterion: 0.00296793. +0.819425 + 1.24522 + [Note] Convergence Criterion: 0.00260915. [Note] Iteration Number: 4. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: -0.819365 - 1.25064 - [Note] Convergence Criterion: 0.000266769. +0.819229 + 1.2499 + [Note] Convergence Criterion: 0.000392739. [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -253,57 +239,58 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: -0.819486 - 1.24728 - [Note] Convergence Criterion: 0.00296793. +0.819425 + 1.24522 + [Note] Convergence Criterion: 0.00260915. [Note] Iteration Number: 4. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: -0.819365 - 1.25064 - [Note] Convergence Criterion: 0.000266769. +0.819229 + 1.2499 + [Note] Convergence Criterion: 0.000392739. [Note] Iteration Number: 5. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: -0.826071 - 1.27092 - [Note] Convergence Criterion: 0.00144698. +0.826009 + 1.27067 + [Note] Convergence Criterion: 0.00149861. [Note] Iteration Number: 6. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: -0.824685 - 1.27635 - [Note] Convergence Criterion: 0.000616001. +0.824662 + 1.27626 + [Note] Convergence Criterion: 0.000627415. [Note] Iteration Number: 7. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: -0.825622 - 1.27394 - [Note] Convergence Criterion: 0.000277843. +0.825615 + 1.27392 + [Note] Convergence Criterion: 0.000275094. [Note] Iteration Number: 8. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: -0.833176 +0.833173 1.29178 - [Note] Convergence Criterion: 0.00114422. + [Note] Convergence Criterion: 0.00114573. [Note] Iteration Number: 9. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: 0.824275 1.28081 - [Note] Convergence Criterion: 0.00104714. + [Note] Convergence Criterion: 0.0010469. [Note] Iteration Number: 10. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: 0.824481 1.26263 - [Note] Convergence Criterion: 0.00215988. + [Note] Convergence Criterion: 0.0021597. [Note] Iteration Number: 11. [Note] Current number of MCMC iterations: 2000. [Note] Parameter Estimates: 0.826653 1.26951 - [Note] Convergence Criterion: 0.000385994. + [Note] Convergence Criterion: 0.000386032. +MCEM pump with transform=TRUE [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -323,9 +310,9 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 2516. [Note] Parameter Estimates: -0.819206 - 1.24929 - [Note] Convergence Criterion: 0.00320976. +0.818707 + 1.24532 + [Note] Convergence Criterion: 0.00263509. [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -345,9 +332,30 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 2516. [Note] Parameter Estimates: -0.819206 - 1.24929 - [Note] Convergence Criterion: 0.00320976. +0.818707 + 1.24532 + [Note] Convergence Criterion: 0.00263509. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 4. + [Note] Current number of MCMC iterations: 9642. + [Note] Parameter Estimates: +0.823945 + 1.26618 + [Note] Convergence Criterion: 0.00101863. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 5. + [Note] Current number of MCMC iterations: 20000. + [Note] Parameter Estimates: +0.824171 + 1.26615 + [Note] Convergence Criterion: 7.78593e-06. +MCEM pump without AD, with transform=FALSE [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -366,9 +374,10 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 2188. [Note] Parameter Estimates: -0.818182 - 1.25068 - [Note] Convergence Criterion: 0.0049061. +0.817876 + 1.24692 + [Note] Convergence Criterion: 0.0041643. +MCEM pump without AD, with transform=TRUE [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -387,9 +396,9 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 2188. [Note] Parameter Estimates: -0.818185 - 1.25068 - [Note] Convergence Criterion: 0.00490643. +0.817883 + 1.24693 + [Note] Convergence Criterion: 0.0041649. [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -408,9 +417,10 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 2188. [Note] Parameter Estimates: -0.818185 - 1.25068 - [Note] Convergence Criterion: 0.00490643. +0.817883 + 1.24693 + [Note] Convergence Criterion: 0.0041649. +MCEM with pump transform=FALSE and some calcNodesOther [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -429,18 +439,18 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 2188. [Note] Parameter Estimates: -1.06877 -2.21092 - [Note] Convergence Criterion: 0.00617455. -Monte Carlo error too big: increasing MCMC sample size. +1.06593 +2.19633 + [Note] Convergence Criterion: 0.00470096. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 4. - [Note] Current number of MCMC iterations: 6197. + [Note] Current number of MCMC iterations: 4298. [Note] Parameter Estimates: -1.07893 -2.25326 - [Note] Convergence Criterion: 0.00137602. + 1.0764 +2.24303 + [Note] Convergence Criterion: 0.00181612. +MCEM with pump transform=FALSE and some calcNodesOther, noAD [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -459,18 +469,18 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 2188. [Note] Parameter Estimates: -1.06877 -2.21093 - [Note] Convergence Criterion: 0.00617451. -Monte Carlo error too big: increasing MCMC sample size. +1.06593 +2.19633 + [Note] Convergence Criterion: 0.00470093. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 4. - [Note] Current number of MCMC iterations: 6197. + [Note] Current number of MCMC iterations: 4298. [Note] Parameter Estimates: -1.07893 -2.25326 - [Note] Convergence Criterion: 0.00137595. + 1.0764 +2.24303 + [Note] Convergence Criterion: 0.00181607. +MCEM simplest 1D works [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -502,13 +512,9 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 4. [Note] Current number of MCMC iterations: 20000. [Note] Parameter Estimates: -4.00374 - [Note] Convergence Criterion: 5.52622e-05. - [Note] Iteration Number: 5. - [Note] Current number of MCMC iterations: 20000. - [Note] Parameter Estimates: -3.99762 - [Note] Convergence Criterion: 1.55888e-05. +4.00925 + [Note] Convergence Criterion: 3.39403e-05. +MCEM simplest 1D with a constrained parameter works [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -525,33 +531,28 @@ Monte Carlo error too big: increasing MCMC sample size. 4.03535 [Note] Convergence Criterion: 0.0125879. Monte Carlo error too big: increasing MCMC sample size. - [Note] Iteration Number: 4. - [Note] Current number of MCMC iterations: 1167. - [Note] Parameter Estimates: -3.94797 - [Note] Convergence Criterion: 0.00145563. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. - [Note] Iteration Number: 5. - [Note] Current number of MCMC iterations: 2611. - [Note] Parameter Estimates: -4.00733 - [Note] Convergence Criterion: 0.000591767. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 4. + [Note] Current number of MCMC iterations: 9400. + [Note] Parameter Estimates: +4.01069 + [Note] Convergence Criterion: 0.000113579. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. - [Note] Iteration Number: 6. + [Note] Iteration Number: 5. [Note] Current number of MCMC iterations: 20000. [Note] Parameter Estimates: -3.99967 - [Note] Convergence Criterion: 1.94222e-05. +4.01258 + [Note] Convergence Criterion: 4.30593e-06. [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -568,33 +569,29 @@ Monte Carlo error too big: increasing MCMC sample size. 4.03535 [Note] Convergence Criterion: 0.0125879. Monte Carlo error too big: increasing MCMC sample size. - [Note] Iteration Number: 4. - [Note] Current number of MCMC iterations: 1167. - [Note] Parameter Estimates: -3.94797 - [Note] Convergence Criterion: 0.00145563. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. - [Note] Iteration Number: 5. - [Note] Current number of MCMC iterations: 2611. - [Note] Parameter Estimates: -4.00733 - [Note] Convergence Criterion: 0.000591767. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 4. + [Note] Current number of MCMC iterations: 9400. + [Note] Parameter Estimates: +4.01069 + [Note] Convergence Criterion: 0.000113579. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. - [Note] Iteration Number: 6. + [Note] Iteration Number: 5. [Note] Current number of MCMC iterations: 20000. [Note] Parameter Estimates: -3.99967 - [Note] Convergence Criterion: 1.94222e-05. +4.01258 + [Note] Convergence Criterion: 4.30593e-06. +MCEM simplest 1D with deterministic intermediates and multiple data works [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -608,91 +605,87 @@ Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. -Monte Carlo error too big: increasing MCMC sample size. -Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 2. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: -1.39417 - [Note] Convergence Criterion: 0.00496179. +1.3734 + [Note] Convergence Criterion: 0.00753366. [Note] Iteration Number: 3. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: -1.10611 - [Note] Convergence Criterion: 0.0182309. +1.2447 + [Note] Convergence Criterion: 0.00713057. [Note] Iteration Number: 4. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: -0.673105 - [Note] Convergence Criterion: 0.038268. +0.931467 + [Note] Convergence Criterion: 0.0232516. [Note] Iteration Number: 5. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: -0.26249 - [Note] Convergence Criterion: 0.0340481. +0.514863 + [Note] Convergence Criterion: 0.0396897. [Note] Iteration Number: 6. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: --0.0389393 - [Note] Convergence Criterion: 0.0201176. +0.0685438 + [Note] Convergence Criterion: 0.045638. [Note] Iteration Number: 7. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: --0.552011 - [Note] Convergence Criterion: 0.0506355. +-0.435173 + [Note] Convergence Criterion: 0.0555813. [Note] Iteration Number: 8. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: --1.0231 - [Note] Convergence Criterion: 0.0413861. +-1.01771 + [Note] Convergence Criterion: 0.0656141. [Note] Iteration Number: 9. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: --1.49505 - [Note] Convergence Criterion: 0.0410266. +-1.33294 + [Note] Convergence Criterion: 0.0212339. [Note] Iteration Number: 10. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: --1.95698 - [Note] Convergence Criterion: 0.0380804. +-1.81544 + [Note] Convergence Criterion: 0.0459213. [Note] Iteration Number: 11. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: --2.30157 - [Note] Convergence Criterion: 0.021926. +-2.22244 + [Note] Convergence Criterion: 0.0324886. [Note] Iteration Number: 12. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: --2.53196 - [Note] Convergence Criterion: 0.0106752. +-2.57693 + [Note] Convergence Criterion: 0.0274394. [Note] Iteration Number: 13. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: --2.71536 - [Note] Convergence Criterion: 0.0072198. +-2.68666 + [Note] Convergence Criterion: 0.00336188. [Note] Iteration Number: 14. - [Note] Current number of MCMC iterations: 12367. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: --2.81595 - [Note] Convergence Criterion: 0.00263783. +-2.79567 + [Note] Convergence Criterion: 0.00352499. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 15. [Note] Current number of MCMC iterations: 20000. [Note] Parameter Estimates: --2.8345 - [Note] Convergence Criterion: 0.000251499. +-2.80731 + [Note] Convergence Criterion: 0.000137197. [Note] Iteration Number: 16. [Note] Current number of MCMC iterations: 20000. [Note] Parameter Estimates: --2.87312 - [Note] Convergence Criterion: 0.000558247. - [Note] Iteration Number: 17. - [Note] Current number of MCMC iterations: 20000. - [Note] Parameter Estimates: --2.88834 - [Note] Convergence Criterion: 0.000177332. +-2.8166 + [Note] Convergence Criterion: 0.000108485. +MCEM 1D with a constrained parameter and deterministic intermediates works [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -808,47 +801,47 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 23. [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -35.0047 - [Note] Convergence Criterion: 0.00513322. +34.9424 + [Note] Convergence Criterion: 0.00398613. [Note] Iteration Number: 24. [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -35.6171 - [Note] Convergence Criterion: 0.0105147. +35.5599 + [Note] Convergence Criterion: 0.0106472. [Note] Iteration Number: 25. [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -36.0888 - [Note] Convergence Criterion: 0.00738943. +36.0364 + [Note] Convergence Criterion: 0.00749461. [Note] Iteration Number: 26. [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -36.4474 - [Note] Convergence Criterion: 0.004954. +36.3993 + [Note] Convergence Criterion: 0.00503566. [Note] Iteration Number: 27. [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -37.4108 - [Note] Convergence Criterion: 0.0215292. +37.3667 + [Note] Convergence Criterion: 0.0216714. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 28. [Note] Current number of MCMC iterations: 1687. [Note] Parameter Estimates: -37.7451 - [Note] Convergence Criterion: 0.0040169. +37.6673 + [Note] Convergence Criterion: 0.0034711. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 29. [Note] Current number of MCMC iterations: 2611. [Note] Parameter Estimates: -38.0069 - [Note] Convergence Criterion: 0.00245493. +37.8566 + [Note] Convergence Criterion: 0.00158479. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 30. [Note] Current number of MCMC iterations: 3315. [Note] Parameter Estimates: -38.2099 - [Note] Convergence Criterion: 0.00159944. +38.0445 + [Note] Convergence Criterion: 0.00144046. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. @@ -857,268 +850,264 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 31. [Note] Current number of MCMC iterations: 12367. [Note] Parameter Estimates: -38.3511 - [Note] Convergence Criterion: 0.000616037. +38.1405 + [Note] Convergence Criterion: 0.00035821. [Note] Iteration Number: 32. [Note] Current number of MCMC iterations: 12367. [Note] Parameter Estimates: -38.4734 - [Note] Convergence Criterion: 0.000497953. +38.2801 + [Note] Convergence Criterion: 0.000602561. [Note] Iteration Number: 33. [Note] Current number of MCMC iterations: 12367. [Note] Parameter Estimates: -38.6567 - [Note] Convergence Criterion: 0.00091512. +38.4794 + [Note] Convergence Criterion: 0.00103894. [Note] Iteration Number: 34. [Note] Current number of MCMC iterations: 12367. [Note] Parameter Estimates: -38.7284 - [Note] Convergence Criterion: 0.000243766. +38.5658 + [Note] Convergence Criterion: 0.000311096. [Note] Iteration Number: 35. [Note] Current number of MCMC iterations: 12367. [Note] Parameter Estimates: -38.901 - [Note] Convergence Criterion: 0.000831227. +38.7518 + [Note] Convergence Criterion: 0.000930625. [Note] Iteration Number: 36. [Note] Current number of MCMC iterations: 12367. [Note] Parameter Estimates: -38.9998 - [Note] Convergence Criterion: 0.000372776. +38.8629 + [Note] Convergence Criterion: 0.000438304. [Note] Iteration Number: 37. [Note] Current number of MCMC iterations: 12367. [Note] Parameter Estimates: -39.0785 - [Note] Convergence Criterion: 0.000274747. +38.9529 + [Note] Convergence Criterion: 0.00032834. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 38. [Note] Current number of MCMC iterations: 20000. [Note] Parameter Estimates: -39.1399 - [Note] Convergence Criterion: 0.000166908. - [Note] Iteration Number: 39. - [Note] Current number of MCMC iterations: 20000. - [Note] Parameter Estimates: -39.1567 - [Note] Convergence Criterion: 3.56707e-05. +39.0072 + [Note] Convergence Criterion: 0.000142358. [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -3.88796 - [Note] Convergence Criterion: 0.199104. +3.61819 + [Note] Convergence Criterion: 0.172391. [Note] Iteration Number: 2. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -6.70707 - [Note] Convergence Criterion: 0.144179. +6.56944 + [Note] Convergence Criterion: 0.156342. [Note] Iteration Number: 3. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -8.93779 - [Note] Convergence Criterion: 0.095655. +9.05379 + [Note] Convergence Criterion: 0.115103. [Note] Iteration Number: 4. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -11.0122 - [Note] Convergence Criterion: 0.0853815. +11.5595 + [Note] Convergence Criterion: 0.11775. [Note] Iteration Number: 5. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -13.5612 - [Note] Convergence Criterion: 0.121944. +14.1227 + [Note] Convergence Criterion: 0.122365. [Note] Iteration Number: 6. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -15.7403 - [Note] Convergence Criterion: 0.0919451. +16.3636 + [Note] Convergence Criterion: 0.0961715. [Note] Iteration Number: 7. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -17.4694 - [Note] Convergence Criterion: 0.0610221. +18.6672 + [Note] Convergence Criterion: 0.101793. [Note] Iteration Number: 8. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -19.1703 - [Note] Convergence Criterion: 0.0596965. +20.5049 + [Note] Convergence Criterion: 0.0673011. [Note] Iteration Number: 9. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -20.5289 - [Note] Convergence Criterion: 0.0404434. +22.0339 + [Note] Convergence Criterion: 0.0509787. [Note] Iteration Number: 10. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -22.3118 - [Note] Convergence Criterion: 0.0623274. +23.1647 + [Note] Convergence Criterion: 0.0304122. [Note] Iteration Number: 11. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -23.4629 - [Note] Convergence Criterion: 0.0311416. +24.1035 + [Note] Convergence Criterion: 0.0233678. [Note] Iteration Number: 12. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -24.8283 - [Note] Convergence Criterion: 0.0416213. +25.4908 + [Note] Convergence Criterion: 0.043811. [Note] Iteration Number: 13. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -26.3874 - [Note] Convergence Criterion: 0.0516532. +26.4272 + [Note] Convergence Criterion: 0.023357. [Note] Iteration Number: 14. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -27.7639 - [Note] Convergence Criterion: 0.0429361. +27.1034 + [Note] Convergence Criterion: 0.0142114. [Note] Iteration Number: 15. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -28.7151 - [Note] Convergence Criterion: 0.0237817. +28.1416 + [Note] Convergence Criterion: 0.0265473. [Note] Iteration Number: 16. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -29.3797 - [Note] Convergence Criterion: 0.0137453. +29.2719 + [Note] Convergence Criterion: 0.0305335. [Note] Iteration Number: 17. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -30.173 - [Note] Convergence Criterion: 0.0178812. +29.8921 + [Note] Convergence Criterion: 0.0123142. [Note] Iteration Number: 18. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: -30.6031 - [Note] Convergence Criterion: 0.00746577. +30.6213 + [Note] Convergence Criterion: 0.0161238. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 19. - [Note] Current number of MCMC iterations: 1000. + [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -31.1897 - [Note] Convergence Criterion: 0.0115362. +30.9228 + [Note] Convergence Criterion: 0.00394432. [Note] Iteration Number: 20. - [Note] Current number of MCMC iterations: 1000. + [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -31.821 - [Note] Convergence Criterion: 0.01316. +31.3233 + [Note] Convergence Criterion: 0.00542803. [Note] Iteration Number: 21. - [Note] Current number of MCMC iterations: 1000. + [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -32.5598 - [Note] Convergence Criterion: 0.0157753. +31.8851 + [Note] Convergence Criterion: 0.0092428. [Note] Iteration Number: 22. - [Note] Current number of MCMC iterations: 1000. + [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -33.2559 - [Note] Convergence Criterion: 0.0149313. +32.4052 + [Note] Convergence Criterion: 0.00841461. [Note] Iteration Number: 23. - [Note] Current number of MCMC iterations: 1000. + [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -33.8701 - [Note] Convergence Criterion: 0.0122046. +32.6848 + [Note] Convergence Criterion: 0.00356841. [Note] Iteration Number: 24. - [Note] Current number of MCMC iterations: 1000. + [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -34.7837 - [Note] Convergence Criterion: 0.0217216. +33.3674 + [Note] Convergence Criterion: 0.0125512. [Note] Iteration Number: 25. - [Note] Current number of MCMC iterations: 1000. + [Note] Current number of MCMC iterations: 1390. [Note] Parameter Estimates: -35.3102 - [Note] Convergence Criterion: 0.00996641. +34.1309 + [Note] Convergence Criterion: 0.0150872. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 26. - [Note] Current number of MCMC iterations: 1167. + [Note] Current number of MCMC iterations: 1687. [Note] Parameter Estimates: -35.665 - [Note] Convergence Criterion: 0.00533183. +34.412 + [Note] Convergence Criterion: 0.00318845. +Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 27. - [Note] Current number of MCMC iterations: 1167. + [Note] Current number of MCMC iterations: 2083. [Note] Parameter Estimates: -36.1096 - [Note] Convergence Criterion: 0.00711728. -Monte Carlo error too big: increasing MCMC sample size. +34.6108 + [Note] Convergence Criterion: 0.00184824. [Note] Iteration Number: 28. - [Note] Current number of MCMC iterations: 1390. + [Note] Current number of MCMC iterations: 2083. [Note] Parameter Estimates: -36.3998 - [Note] Convergence Criterion: 0.00374219. +35.0625 + [Note] Convergence Criterion: 0.00582845. [Note] Iteration Number: 29. - [Note] Current number of MCMC iterations: 1390. + [Note] Current number of MCMC iterations: 2083. [Note] Parameter Estimates: -36.962 - [Note] Convergence Criterion: 0.00938013. -Monte Carlo error too big: increasing MCMC sample size. +35.4323 + [Note] Convergence Criterion: 0.00435585. [Note] Iteration Number: 30. - [Note] Current number of MCMC iterations: 1687. + [Note] Current number of MCMC iterations: 2083. [Note] Parameter Estimates: -37.3095 - [Note] Convergence Criterion: 0.00434664. +35.9355 + [Note] Convergence Criterion: 0.00684756. +Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 31. - [Note] Current number of MCMC iterations: 1687. + [Note] Current number of MCMC iterations: 2611. [Note] Parameter Estimates: -37.7594 - [Note] Convergence Criterion: 0.00620774. +36.2229 + [Note] Convergence Criterion: 0.00278768. [Note] Iteration Number: 32. - [Note] Current number of MCMC iterations: 1687. + [Note] Current number of MCMC iterations: 2611. [Note] Parameter Estimates: -37.9889 - [Note] Convergence Criterion: 0.00245218. -Monte Carlo error too big: increasing MCMC sample size. +36.5298 + [Note] Convergence Criterion: 0.00302671. [Note] Iteration Number: 33. - [Note] Current number of MCMC iterations: 2083. + [Note] Current number of MCMC iterations: 2611. [Note] Parameter Estimates: -38.2551 - [Note] Convergence Criterion: 0.0027239. -Monte Carlo error too big: increasing MCMC sample size. -Monte Carlo error too big: increasing MCMC sample size. +36.959 + [Note] Convergence Criterion: 0.0049841. [Note] Iteration Number: 34. - [Note] Current number of MCMC iterations: 3315. + [Note] Current number of MCMC iterations: 2611. [Note] Parameter Estimates: -38.4409 - [Note] Convergence Criterion: 0.00139352. +37.2769 + [Note] Convergence Criterion: 0.00318619. [Note] Iteration Number: 35. - [Note] Current number of MCMC iterations: 3315. + [Note] Current number of MCMC iterations: 2611. [Note] Parameter Estimates: -38.6373 - [Note] Convergence Criterion: 0.00150371. +37.7139 + [Note] Convergence Criterion: 0.00520789. [Note] Iteration Number: 36. - [Note] Current number of MCMC iterations: 3315. + [Note] Current number of MCMC iterations: 2611. [Note] Parameter Estimates: -38.8325 - [Note] Convergence Criterion: 0.00149478. -Monte Carlo error too big: increasing MCMC sample size. -Monte Carlo error too big: increasing MCMC sample size. -Monte Carlo error too big: increasing MCMC sample size. -Monte Carlo error too big: increasing MCMC sample size. -Monte Carlo error too big: increasing MCMC sample size. -Monte Carlo error too big: increasing MCMC sample size. -Monte Carlo error too big: increasing MCMC sample size. +38.0998 + [Note] Convergence Criterion: 0.00426572. [Note] Iteration Number: 37. - [Note] Current number of MCMC iterations: 20000. + [Note] Current number of MCMC iterations: 2611. [Note] Parameter Estimates: -38.9006 - [Note] Convergence Criterion: 0.000192296. +38.2964 + [Note] Convergence Criterion: 0.00165787. [Note] Iteration Number: 38. - [Note] Current number of MCMC iterations: 20000. + [Note] Current number of MCMC iterations: 2611. [Note] Parameter Estimates: -39.0068 - [Note] Convergence Criterion: 0.00035607. +38.5742 + [Note] Convergence Criterion: 0.00268214. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 39. - [Note] Current number of MCMC iterations: 20000. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: -39.0837 - [Note] Convergence Criterion: 0.000225167. +38.6901 + [Note] Convergence Criterion: 0.000555406. [Note] Iteration Number: 40. - [Note] Current number of MCMC iterations: 20000. + [Note] Current number of MCMC iterations: 7175. [Note] Parameter Estimates: -39.2204 - [Note] Convergence Criterion: 0.00051648. +38.7954 + [Note] Convergence Criterion: 0.000493818. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 41. [Note] Current number of MCMC iterations: 20000. [Note] Parameter Estimates: -39.2688 - [Note] Convergence Criterion: 0.000124143. +38.802 + [Note] Convergence Criterion: 1.29801e-05. +MCEM simplest 2x1D works, with multiple data for each [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -1146,8 +1135,9 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 5. [Note] Current number of MCMC iterations: 2083. [Note] Parameter Estimates: -9.65659 - [Note] Convergence Criterion: 0.00057379. +9.65931 + [Note] Convergence Criterion: 0.000549532. +MCMC for simple LME case works [Note] Iteration Number: 1. [Note] Current number of MCMC iterations: 1000. [Note] Parameter Estimates: @@ -1396,95 +1386,216 @@ Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 28. [Note] Current number of MCMC iterations: 1546. [Note] Parameter Estimates: - 2.33427 -0.135925 - 9.71458 -0.288106 - 0.18677 - [Note] Convergence Criterion: 0.00610788. + 2.33425 +0.136228 + 9.71434 +0.288558 +0.186775 + [Note] Convergence Criterion: 0.00466801. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 29. [Note] Current number of MCMC iterations: 1962. [Note] Parameter Estimates: - 2.33724 -0.134264 - 9.71506 -0.288378 -0.187002 - [Note] Convergence Criterion: 0.00439994. + 2.33717 +0.134718 + 9.71501 +0.288453 +0.186971 + [Note] Convergence Criterion: 0.00387192. [Note] Iteration Number: 30. [Note] Current number of MCMC iterations: 1962. [Note] Parameter Estimates: - 2.33678 -0.132272 - 9.71363 -0.289139 -0.186554 - [Note] Convergence Criterion: 0.00660849. + 2.33672 +0.132695 + 9.7136 +0.289208 +0.186539 + [Note] Convergence Criterion: 0.00667758. [Note] Iteration Number: 31. [Note] Current number of MCMC iterations: 1962. [Note] Parameter Estimates: - 2.33568 -0.130075 - 9.71343 -0.290403 -0.187144 - [Note] Convergence Criterion: 0.00787739. + 2.33563 +0.130473 + 9.7134 +0.290466 +0.187132 + [Note] Convergence Criterion: 0.00794951. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 32. [Note] Current number of MCMC iterations: 3255. [Note] Parameter Estimates: - 2.33573 - 0.129 - 9.71202 -0.291863 -0.187177 - [Note] Convergence Criterion: 0.00342782. + 2.33567 +0.129721 + 9.7121 +0.291513 +0.187156 + [Note] Convergence Criterion: 0.00211147. [Note] Iteration Number: 33. [Note] Current number of MCMC iterations: 3255. [Note] Parameter Estimates: - 2.33624 -0.127576 - 9.71232 -0.292143 -0.186954 - [Note] Convergence Criterion: 0.00305783. + 2.33617 +0.128259 + 9.71249 +0.291841 +0.186934 + [Note] Convergence Criterion: 0.00316764. [Note] Iteration Number: 34. [Note] Current number of MCMC iterations: 3255. [Note] Parameter Estimates: - 2.33712 -0.125168 - 9.71289 -0.292336 -0.187388 - [Note] Convergence Criterion: 0.00699576. + 2.33704 + 0.12581 + 9.71303 +0.292071 + 0.18737 + [Note] Convergence Criterion: 0.00713368. [Note] Iteration Number: 35. [Note] Current number of MCMC iterations: 3255. [Note] Parameter Estimates: - 2.3366 -0.124069 - 9.71384 -0.290924 -0.187017 - [Note] Convergence Criterion: 0.00385993. + 2.33652 +0.124687 + 9.71398 +0.290678 +0.187001 + [Note] Convergence Criterion: 0.00386159. [Note] Iteration Number: 36. [Note] Current number of MCMC iterations: 3255. [Note] Parameter Estimates: - 2.33804 -0.124058 - 9.71217 -0.293027 -0.187268 - [Note] Convergence Criterion: 0.00386721. + 2.33798 +0.124657 + 9.71228 +0.292821 +0.187255 + [Note] Convergence Criterion: 0.00395005. +Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. Monte Carlo error too big: increasing MCMC sample size. [Note] Iteration Number: 37. - [Note] Current number of MCMC iterations: 5554. - [Note] Parameter Estimates: - 2.33613 -0.123456 - 9.71184 -0.292731 -0.187387 - [Note] Convergence Criterion: 0.000937941. + [Note] Current number of MCMC iterations: 7306. + [Note] Parameter Estimates: + 2.33632 +0.124165 + 9.71177 +0.293147 +0.187335 + [Note] Convergence Criterion: 0.000695707. +MCEM thin control entry works (simplest 1D case) + [Note] Iteration Number: 1. + [Note] Current number of MCMC iterations: 1000. + [Note] Parameter Estimates: +1.33036 + [Note] Convergence Criterion: 0.146406. + [Note] Stopping MCEM because maximum number of MCEM iterations +(maxIter=1) was reached. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 1. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +1.18128 + [Note] Convergence Criterion: 0.00836943. + [Note] Iteration Number: 2. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +0.978408 + [Note] Convergence Criterion: 0.0134633. + [Note] Iteration Number: 3. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +0.678763 + [Note] Convergence Criterion: 0.0208968. + [Note] Iteration Number: 4. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +0.23017 + [Note] Convergence Criterion: 0.0460741. + [Note] Iteration Number: 5. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +-0.155733 + [Note] Convergence Criterion: 0.0331282. + [Note] Iteration Number: 6. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +-0.709289 + [Note] Convergence Criterion: 0.0655549. + [Note] Iteration Number: 7. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +-1.24358 + [Note] Convergence Criterion: 0.0564287. + [Note] Iteration Number: 8. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +-1.67465 + [Note] Convergence Criterion: 0.0372498. + [Note] Iteration Number: 9. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +-2.04859 + [Note] Convergence Criterion: 0.0297808. + [Note] Iteration Number: 10. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +-2.34762 + [Note] Convergence Criterion: 0.0187838. + [Note] Iteration Number: 11. + [Note] Current number of MCMC iterations: 6823. + [Note] Parameter Estimates: +-2.74534 + [Note] Convergence Criterion: 0.0311541. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 12. + [Note] Current number of MCMC iterations: 11742. + [Note] Parameter Estimates: +-2.80128 + [Note] Convergence Criterion: 0.00117352. +Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 13. + [Note] Current number of MCMC iterations: 15490. + [Note] Parameter Estimates: +-2.85236 + [Note] Convergence Criterion: 0.000958239. +Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 14. + [Note] Current number of MCMC iterations: 20000. + [Note] Parameter Estimates: +-2.87815 + [Note] Convergence Criterion: 0.000367584. +MCEM thin control and simplest 1D with deterministic intermediates and multiple data works + [Note] Iteration Number: 1. + [Note] Current number of MCMC iterations: 1000. + [Note] Parameter Estimates: +1.33442 + [Note] Convergence Criterion: 0.145608. + [Note] Stopping MCEM because maximum number of MCEM iterations +(maxIter=1) was reached. + [Note] Iteration Number: 1. + [Note] Current number of MCMC iterations: 1000. + [Note] Parameter Estimates: +1.33442 + [Note] Convergence Criterion: 0.145608. + [Note] Stopping MCEM because maximum number of MCEM iterations +(maxIter=1) was reached. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. +Monte Carlo error too big: increasing MCMC sample size. + [Note] Iteration Number: 1. + [Note] Current number of MCMC iterations: 535. + [Note] Parameter Estimates: +1.47333 + [Note] Convergence Criterion: 0.11241. + [Note] Stopping MCEM because maximum number of MCEM iterations +(maxIter=1) was reached. + [Note] Iteration Number: 1. + [Note] Current number of MCMC iterations: 1000. + [Note] Parameter Estimates: +1.33036 + [Note] Convergence Criterion: 0.146406. + [Note] Stopping MCEM because maximum number of MCEM iterations +(maxIter=1) was reached. diff --git a/packages/nimble/tests/testthat/test-mcem.R b/packages/nimble/tests/testthat/test-mcem.R index a439aa5d3..44be92f1c 100644 --- a/packages/nimble/tests/testthat/test-mcem.R +++ b/packages/nimble/tests/testthat/test-mcem.R @@ -9,8 +9,7 @@ nimbleOptions(verbose = TRUE) nimbleProgressBarSetting <- nimbleOptions('MCMCprogressBar') nimbleOptions(MCMCprogressBar = FALSE) -context("Testing of MCEM") - +# tests begin with a cat() to help organize the gold file output goldFileName <- 'mcemTestLog_Correct.Rout' tempFileName <- 'mcemTestLog.Rout' generatingGoldFile <- !is.null(nimbleOptions('generateGoldFileForMCEMtesting')) @@ -65,6 +64,7 @@ pump_mle_trans <- pump_optim_trans$par pump_vcov_trans <- inverse(-pump_optim_trans$hessian) test_that("MCEM error-trapping", { + cat("MCEM error-trapping\n") pumpCode <- nimbleCode({ for (i in 1:N){ theta[i] ~ dgamma(alpha,beta); @@ -101,6 +101,7 @@ test_that("MCEM error-trapping", { }) test_that("MCEM pump with all defaults (roxygen example)", { + cat("MCEM pump with all defaults (roxygen example)\n") pumpCode <- nimbleCode({ for (i in 1:N){ theta[i] ~ dgamma(alpha,beta); @@ -142,6 +143,7 @@ test_that("MCEM pump with all defaults (roxygen example)", { }) test_that("MCEM with pump transform=FALSE",{ + cat("MCEM with pump transform=FALSE\n") pumpCode <- nimbleCode({ for (i in 1:N){ theta[i] ~ dgamma(alpha,beta) @@ -239,6 +241,7 @@ test_that("MCEM with pump transform=FALSE",{ }) test_that("MCEM pump with transform=TRUE", { + cat("MCEM pump with transform=TRUE\n") # transformed version: Use parameter transformation ## build an MCEM algorithm with Ascent-based convergence criterion pumpCode <- nimbleCode({ @@ -286,7 +289,7 @@ test_that("MCEM pump with transform=TRUE", { set.seed(0) cpump$alpha <- pump$alpha cpump$beta <- pump$beta - out <- cpumpMCEM$findMLE(initM = 1000, returnTrans=TRUE) + out <- cpumpMCEM$findMLE(initM = 1000, returnTrans=TRUE, tol = 0.001) expect_equal(out$par, pump_mle_trans, tol=0.01) vc <- cpumpMCEM$vcov(out$par, trans = TRUE) expect_equal(vc, pump_vcov_trans, tol=0.02) @@ -295,6 +298,7 @@ test_that("MCEM pump with transform=TRUE", { }) test_that("MCEM pump without AD, with transform=FALSE", { + cat("MCEM pump without AD, with transform=FALSE\n") #### # version without AD pumpCode <- nimbleCode({ @@ -343,6 +347,7 @@ test_that("MCEM pump without AD, with transform=FALSE", { test_that("MCEM pump without AD, with transform=TRUE", { + cat("MCEM pump without AD, with transform=TRUE\n") #### # version without AD pumpCode <- nimbleCode({ @@ -398,6 +403,7 @@ test_that("MCEM pump without AD, with transform=TRUE", { ## are not related to missing data or latent states test_that("MCEM with pump transform=FALSE and some calcNodesOther",{ + cat("MCEM with pump transform=FALSE and some calcNodesOther\n") pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) @@ -486,6 +492,7 @@ test_that("MCEM with pump transform=FALSE and some calcNodesOther",{ test_that("MCEM with pump transform=FALSE and some calcNodesOther, noAD",{ + cat("MCEM with pump transform=FALSE and some calcNodesOther, noAD\n") pumpConsts <- list(N = 10, t = c(94.3, 15.7, 62.9, 126, 5.24, 31.4, 1.05, 1.05, 2.1, 10.5)) @@ -578,6 +585,7 @@ test_that("MCEM with pump transform=FALSE and some calcNodesOther, noAD",{ # Adopting some (far from all) tests from test-ADlaplace test_that("MCEM simplest 1D works", { + cat("MCEM simplest 1D works\n") m <- nimbleModel( nimbleCode({ y ~ dnorm(a, sd = 2) @@ -598,6 +606,7 @@ test_that("MCEM simplest 1D works", { }) test_that("MCEM simplest 1D with a constrained parameter works", { + cat("MCEM simplest 1D with a constrained parameter works\n") m <- nimbleModel( nimbleCode({ y ~ dnorm(a, sd = 2) @@ -652,6 +661,7 @@ test_that("MCEM simplest 1D with a constrained parameter works", { # that we don't need to test here. test_that("MCEM simplest 1D with deterministic intermediates and multiple data works", { + cat("MCEM simplest 1D with deterministic intermediates and multiple data works\n") set.seed(1) m <- nimbleModel( nimbleCode({ @@ -706,6 +716,7 @@ test_that("MCEM simplest 1D with deterministic intermediates and multiple data w # Case "1D with deterministic intermediates works" from Laplace seems redundant test_that("MCEM 1D with a constrained parameter and deterministic intermediates works", { + cat("MCEM 1D with a constrained parameter and deterministic intermediates works\n") m <- nimbleModel( nimbleCode({ y ~ dnorm(0.2 * a, sd = 2) @@ -773,6 +784,7 @@ test_that("MCEM 1D with a constrained parameter and deterministic intermediates # Skipping "with deterministic intermediates and multiple data works" test_that("MCEM simplest 2x1D works, with multiple data for each", { + cat("MCEM simplest 2x1D works, with multiple data for each\n") set.seed(1) y <- matrix(rnorm(6, 4, 5), nrow = 2) m <- nimbleModel( @@ -876,6 +888,7 @@ test_that("MCEM simplest 2x1D works, with multiple data for each", { ## }) test_that("MCMC for simple LME case works", { + cat("MCMC for simple LME case works\n") set.seed(1) g <- rep(1:10, each = 5) n <- length(g) @@ -942,6 +955,155 @@ test_that("MCMC for simple LME case works", { } }) +## The next two tests are added as regression tests for fixes made in 1.4.0. +## Prior to 1.4.0, there were several small problems: +## - control list "thin" element was not used during sampling +## - when increasing M, new samples were not appended by rather a full new sample was generated +## - in vcov with resetSamples=TRUE, thin was not used. +## - burnin was applied post-thinning instead of pre-thinning (which would be consistent with runMCMC etc.) +## - reset vs resetMV was not used, although I don't have a regression test for that here. +## - nimDerivs did not have reset=TRUE for the first all in a findMLe run, although I don't have a regression test for that here. + +test_that("MCEM thin control entry works (simplest 1D case)",{ + cat("MCEM thin control entry works (simplest 1D case)\n") + nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions=TRUE) + on.exit(nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions=FALSE)) + set.seed(1) + m <- nimbleModel( + nimbleCode({ + mu ~ dnorm(0, sd = 5) + a ~ dexp(rate = exp(0.5 * mu)) + for (i in 1:5){ + y[i] ~ dnorm(0.2 * a, sd = 2) + } + }), data = list(y = rnorm(5, 1, 2)), inits = list(mu = 2, a = 1), + buildDerivs = TRUE + ) + MCEM <- buildMCEM(model = m, control=list(thin=2)) + cm <- compileNimble(m) + cMCEM <- compileNimble(MCEM, project = m) + set.seed(1) + expect_no_error(opt <- cMCEM$findMLE(maxIter = 1, ascent=FALSE)) # PASSES only as of 1.4.0 + # We should have "effective burnin" = burnin/thin = 250; this is E$burnIn + expect_equal(cMCEM$E$burnIn, 250) # PASSES only as of 1.4.0 + opt <- cMCEM$findMLE(C = 2, initM = 2000, burnin = 500, thin = 2) + expect_equal(opt$par, -2.856, .04) + vcov <- cMCEM$vcov(opt$par, trans=FALSE, resetSamples = TRUE, M = 2000) + expect_equal(dim(as.matrix(cMCEM$mcmc_Latent$mvSamples))[1], 2000/2) # PASSES only as of 1.4.0 +}) + +test_that("MCEM thin control and simplest 1D with deterministic intermediates and multiple data works", { + cat("MCEM thin control and simplest 1D with deterministic intermediates and multiple data works\n") + nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions=TRUE) + on.exit(nimbleOptions(buildInterfacesForCompiledNestedNimbleFunctions=FALSE)) + set.seed(1) + m <- nimbleModel( + nimbleCode({ + mu ~ dnorm(0, sd = 5) + a ~ dexp(rate = exp(0.5 * mu)) + for (i in 1:5){ + y[i] ~ dnorm(0.2 * a, sd = 2) + } + }), data = list(y = rnorm(5, 1, 2)), inits = list(mu = 2, a = 1), + buildDerivs = TRUE + ) + MCEM <- buildMCEM(model = m) + cm <- compileNimble(m) + cMCEM <- compileNimble(MCEM, project = m) + + expect_identical(cMCEM$getParamNodes()[1], "mu") # NEW FEATURE in 1.4.0 + + ## First, check that MCMC is replicable. + ## This should work and is being checked to be cautious. + for(v in m$getVarNames()) + cm[[v]] <- m[[v]] + cm$calculate() + cm$mu <- 2 + set.seed(1) + cMCEM$doMCMC(1000, 1, TRUE) + firstSamples <- as.matrix(cMCEM$mcmc_Latent$mvSamples) + + for(v in m$getVarNames()) + cm[[v]] <- m[[v]] + cm$calculate() + cm$mu <- 2 + set.seed(1) + cMCEM$doMCMC(1000, 1, TRUE) + secondSamples <- as.matrix(cMCEM$mcmc_Latent$mvSamples) + + expect_identical(firstSamples, secondSamples) + + ## Next, check that with reset=FALSE, samples are appended, so the + ## first 1000 are the same as above. This is a regression test for + ## nimble 1.4.0. Prior to this version, samples were not appended, + ## but rather a full new sample was obtained. This was inefficient but + ## not wrong. + cMCEM$doMCMC(500, 1, FALSE) + thirdSamples <- as.matrix(cMCEM$mcmc_Latent$mvSamples) + expect_identical(firstSamples[1:1000], thirdSamples[1:1000, 1]) + + ## Second, test that samples are reproducible and + ## additional samples are appended when running fully + ## in the algorithm. + for(v in m$getVarNames()) + cm[[v]] <- m[[v]] + cm$calculate() + cm$mu <- 2 + set.seed(1) + opt <- cMCEM$findMLE(maxIter = 1, ascent=FALSE) + firstSamples <- as.matrix(cMCEM$mcmc_Latent$mvSamples) + + for(v in m$getVarNames()) + cm[[v]] <- m[[v]] + cm$calculate() + cm$mu <- 2 + set.seed(1) + opt <- cMCEM$findMLE(maxIter = 1, ascent=FALSE) + secondSamples <- as.matrix(cMCEM$mcmc_Latent$mvSamples) + + expect_identical(firstSamples, secondSamples) + + for(v in m$getVarNames()) + cm[[v]] <- m[[v]] + cm$calculate() + cm$mu <- 2 + set.seed(1) + # by choosing only 10 more samples than burnin, we are + # forcing the adaptive increase of samples and can check that + # it appends to rather than replaces previous samples. + opt <- cMCEM$findMLE(initM = 510, maxIter = 1, ascent=TRUE ) + thirdSamples <- as.matrix(cMCEM$mcmc_Latent$mvSamples) + nThird <- dim(thirdSamples)[1] + expect_identical(firstSamples[1:nThird], thirdSamples[1:nThird, 1]) # PASSES only as of 1.4.0 + + ## Next, check that thinning is respected + for(v in m$getVarNames()) + cm[[v]] <- m[[v]] + cm$mu <- 2 + cm$calculate() + set.seed(1) + opt <- cMCEM$findMLE(maxIter = 1, ascent=FALSE, thin = 2) ## PASSES only as of 1.4.0 + firstSamplesT <- as.matrix(cMCEM$mcmc_Latent$mvSamples) + expect_equal(firstSamples[seq(2, 1000, by=2)], firstSamplesT[1:500]) + expect_equal(dim(firstSamplesT)[1], 500) + expect_equal(cMCEM$E$burnIn, 250) # PASSES only as of 1.4.0 + + # and that appending samples works correctly with thinning + # Nope: This test is disabled because in fact when niter is not a + # multiple of thin, even with reset=FALSE the MCMC samples will + # not match but will be offset by niter %% thin at the junction + # of the first and second runs. + ## for(v in m$getVarNames()) + ## cm[[v]] <- m[[v]] + ## cm$mu <- 2 + ## cm$calculate() + ## set.seed(1) + ## opt <- cMCEM$findMLE(initM = 510, maxIter = 1, ascent=TRUE, thin = 2) + ## thirdSamplesT <- as.matrix(cMCEM$mcmc_Latent$mvSamples) + ## nThird <- dim(thirdSamplesT)[1] + ## expect_identical(firstSamplesT[1:nThird], thirdSamplesT[1:nThird, 1]) +}) + sink(NULL) if(!generatingGoldFile) {