set.seed(1)
y1000<-rnorm(n=1000,mean=600,sd=30)
ModelData <-list(mass = y1000)
ModelConsts <- list(nobs = length(y1000))
ModelCode<-nimbleCode(
{
# Priors
population.mean ~ dunif(0,5000)
#population.mean ~ dnorm(0,sd=100)
population.sd ~ dunif(0,100)
# Normal distribution parameterized by precision = 1/variance in Nimble
population.variance <-
population.sd *
population.sd precision <- 1 / population.variance
# Likelihood
for(i in 1:nobs){
mass[i] ~ dnorm(population.mean, precision)
}
})
ModelInits <- function()
{list (population.mean = rnorm(1,600,90),
population.sd = runif(1, 1, 30))}
Nchains <- 3
set.seed(1)
Inits<-lapply(1:Nchains,function(x){ModelInits()})
#specifying the names of parameters to analyse and save:
params <- c("population.mean", "
population.sd")
calculations.for.modalDIC<-expression(
{
Model1.EC<-Model[[1]]
## first preparing the sampled parameters in a matrix format; uses the as.matrix function specific to mcmc.list objects; also stocking names of the parameters
samples.List.matrix.EC<-as.matrix(samplesList.temp)
names.samples.EC<-dimnames(
samples.List.matrix.EC)[[2]]
## second preparing the names of variables to calculate on:
varNames.EC<-Model1.EC$getVarNames()
DatavarNames.EC<-names(data)
notDatavarNames.EC<-setdiff(varNames.EC,DatavarNames.EC)
## third writing and compiling the nimbleFunction we will use:
logProbCalc.EC <- nimbleFunction(
setup = function(model,names.ref.list,notDatavarNames,DatavarNames) {
},
run = function(P = double(1)) { ##NB: double(1) means this if of double type and has one dimension
values(model,names.ref.list) <<- P
model$calculate(notDatavarNames)
return(model$calculate(DatavarNames))
returnType(double(0))
})
logProbCalcPrepared.EC <- logProbCalc.EC(Model1.EC,
names.samples.EC, notDatavarNames.EC, DatavarNames.EC)
ClogProbCalcPrepared.EC <- compileNimble(Model1.EC, logProbCalcPrepared.EC)
## fourth, running through all the samples in a lapply function:
logLiks.EC<-sapply(1:(dim(
samples.List.matrix.EC)[1]),function(toto)
{
ClogProbCalcPrepared.EC$logProbCalcPrepared.EC$run(
samples.List.matrix.EC[toto,])
})
#mode type DIC; cf. Celeux et al. 2006 Bayesian analysis
DIC.mode.EC<--4*mean(logLiks.EC)+2*max(logLiks.EC)
p.DIC.mode.EC<--2*mean(logLiks.EC)+2*max(logLiks.EC)
#calculation of classical DIC; cf. Celeux et al. 2006 Bayesian analysis
logLiks.meanparams.EC<-ClogProbCalcPrepared.EC$logProbCalcPrepared.EC$run(colMeans(
samples.List.matrix.EC))
DIC.EC<--4*mean(logLiks.EC)+2*
logLiks.meanparams.EC p.DIC.EC<--2*mean(logLiks.EC)+2*
logLiks.meanparams.EC list(DIC.mode=
DIC.mode.EC,p.DIC.mode=
p.DIC.mode.EC,DIC=
DIC.EC,p.DIC=
p.DIC.EC)
}
)
out.mcmc<-runMCMC_btadjust(code=ModelCode, constants = ModelConsts, data = ModelData, MCMC_language="Nimble",
Nchains=Nchains, params=params, inits=Inits,
niter.min=1000, niter.max=300000,
nburnin.min=100, nburnin.max=200000,
thin.min=1, thin.max=1000,
conv.max=1.05, neff.min=1000,
control=list(print.diagnostics=FALSE),
control.MCMC=list(includeParentNodes=TRUE,extraCalculations=calculations.for.modalDIC))
attributes(out.mcmc)<-c(attributes(out.mcmc),DIC=NA)
The run ends up - in Linux but not in Windows- with the following messages that are unclear to me:
[1] "###################################################################################"
[1] "Performing extra calculations"
[1] "###################################################################################"
Compiling
[Note] This may take a minute.
[Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
Warning: 8 objects were found from a DLL
Called from: cppProjects[[i]]$unloadSO(check = TRUE, force = FALSE)
Warning message:
In cppProjects[[i]]$unloadSO(check = TRUE, force = FALSE) :
A DLL to be unloaded has non-zero (population_dot_sd_L2_UID_5__1, population_dot_mean_L1_UID_4__1, Nimble__1, population_dot_variance_L3_UID_6__1, mass_L6_UID_9__1, lifted_d1_over_sqrt_oPprecision_cP_L6_UID_8__1, precision_L4_UID_7__1, numberedObjects) objects that need to be finalized first. It was not unloaded.
Browse[1]>
debug: if (!force) return(NULL)
Browse[2]> attributes(out.mcmc)<-c(attributes(out.mcmc),DIC=NA)
Error: object 'out.mcmc' not found
Browse[2]>
debug: return(NULL)
I have tried several ways of installing nimble - except from github which I am unsure I can install on this station - but the result was the same with all of them.
Thank you in advance for your help.