1 | jags.fit(inData, inInits, pars.to.save, model, model.file, n.chains, niters, conv.limit, setsize, nruns, Covs)
|
inData |
|
inInits |
|
pars.to.save |
|
model |
|
model.file |
|
n.chains |
|
niters |
|
conv.limit |
|
setsize |
|
nruns |
|
Covs |
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 | ##---- Should be DIRECTLY executable !! ----
##-- ==> Define data, use random,
##-- or do help(data=index) for the standard data sets.
## The function is currently defined as
function (inData, inInits, pars.to.save, model, model.file, n.chains,
niters, conv.limit, setsize, nruns, Covs)
{
mod = jags.model(model.file, data = inData, inits = inInits,
n.chains = n.chains, n.adapt = 0)
done.adapt = adapt(mod, n.iter = 2 * setsize, end.adaptation = FALSE)
while (!done.adapt) done.adapt = adapt(mod, n.iter = 2 *
setsize, end.adaptation = FALSE)
samples <- coda.samples(model = mod, variable.names = pars.to.save,
n.iter = setsize, thin = 1)
nsamples = setsize
varnames <- dimnames(samples[[3]])[[2]]
alpha.vars <- grep("alpha", varnames)
beta.vars <- grep("beta", varnames)
Sd.vars <- grep("Sd", varnames)
if (!is.null(Covs))
slope.vars <- grep("slope", varnames)
Mean.vars <- grep("theta", varnames)
if (is.null(Covs)) {
temp.samples.array = array(NA, c(niters, n.chains, length(alpha.vars) +
length(beta.vars) + length(Sd.vars)))
temp.samples.array[seq(nsamples), 1, ] <- samples[[1]][,
c(alpha.vars, beta.vars, Sd.vars)]
temp.samples.array[seq(nsamples), 2, ] <- samples[[2]][,
c(alpha.vars, beta.vars, Sd.vars)]
temp.samples.array[seq(nsamples), 3, ] <- samples[[3]][,
c(alpha.vars, beta.vars, Sd.vars)]
dimnames(temp.samples.array) <- list(NULL, NULL, varnames[c(alpha.vars,
beta.vars, Sd.vars)])
}
else {
temp.samples.array = array(NA, c(niters, n.chains, length(alpha.vars) +
length(beta.vars) + length(Sd.vars) + length(slope.vars)))
temp.samples.array[seq(nsamples), 1, ] <- samples[[1]][,
c(alpha.vars, slope.vars, beta.vars, Sd.vars)]
temp.samples.array[seq(nsamples), 2, ] <- samples[[2]][,
c(alpha.vars, slope.vars, beta.vars, Sd.vars)]
temp.samples.array[seq(nsamples), 3, ] <- samples[[3]][,
c(alpha.vars, slope.vars, beta.vars, Sd.vars)]
dimnames(temp.samples.array) <- list(NULL, NULL, varnames[c(Sd.vars,
slope.vars, beta.vars, Sd.vars)])
}
max.bgrRatio = max(apply(temp.samples.array[seq(nsamples),
, , drop = F], 3, maxbgrRatio))
check <- max.bgrRatio > conv.limit
print(max.bgrRatio)
flush.console()
if (check) {
count <- 1
while (check & (count < niters/setsize)) {
samples <- coda.samples(mod, variable.names = pars.to.save,
n.iter = setsize, thin = 1)
count <- count + 1
if (is.null(Covs)) {
temp.samples.array[nsamples + seq(setsize), 1,
] = samples[[1]][, c(alpha.vars, beta.vars,
Sd.vars)]
temp.samples.array[nsamples + seq(setsize), 2,
] = samples[[2]][, c(alpha.vars, beta.vars,
Sd.vars)]
temp.samples.array[nsamples + seq(setsize), 3,
] = samples[[3]][, c(alpha.vars, beta.vars,
Sd.vars)]
}
else {
temp.samples.array[nsamples + seq(setsize), 1,
] = samples[[1]][, c(alpha.vars, slope.vars,
beta.vars, Sd.vars)]
temp.samples.array[nsamples + seq(setsize), 2,
] = samples[[2]][, c(alpha.vars, slope.vars,
beta.vars, Sd.vars)]
temp.samples.array[nsamples + seq(setsize), 3,
] = samples[[3]][, c(alpha.vars, slope.vars,
beta.vars, Sd.vars)]
}
nsamples = nsamples + setsize
max.bgrRatio = max(apply(temp.samples.array[seq(nsamples),
, , drop = F], 3, maxbgrRatio))
check <- max.bgrRatio > conv.limit
print(max.bgrRatio)
flush.console()
}
}
no.to.converge = nsamples
print(max.bgrRatio)
flush.console()
no.to.keep <- no.to.converge/2
if (nruns > no.to.keep)
samples = coda.samples(mod, variable.names = pars.to.save,
n.iter = nruns, thin = 1)
else {
thin <- floor(no.to.keep/nruns)
samples <- coda.samples(mod, variable.names = pars.to.save,
n.iter = no.to.keep, n.thin = thin)
}
samples.array = array(NA, c(dim(samples[[1]])[1], n.chains,
length(varnames)))
samples.array[, 1, ] = samples[[1]]
samples.array[, 2, ] = samples[[2]]
samples.array[, 3, ] = samples[[3]]
dimnames(samples.array) = list(NULL, NULL, varnames)
out <- list(no.to.converge, dim(samples.array)[1], samples.array)
names(out) <- c("BurnIn", "No. Runs Per Chain", "Samples")
return(out)
}
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.