# this function calls on other functions in order to return the sampled estimates
# and the credible intervals together with the posterior distribution objects
# to be passed on for forther analysis
gibbsFun <- function(data, estimates, n.iter, n.burnin, thin, n.chains, interval, item.dropped, pairwise,
callback = function(){}, k0, df0, a0, b0, m0) {
p <- ncol(data)
res <- list()
if ("alpha" %in% estimates || "lambda2" %in% estimates || "lambda4" %in% estimates || "lambda6" %in% estimates ||
"glb" %in% estimates) {
tmp_out <- covSamp(data, n.iter, n.burnin, thin, n.chains, pairwise, callback, k0, df0)
C <- tmp_out$cov_mat
res$covsamp <- C
res$data_mis_samp_cov <- tmp_out$dat_mis_samp_cov
if (item.dropped) {
Ctmp <- array(0, c(n.chains, length(seq(1, n.iter - n.burnin, thin)), p, p - 1, p - 1))
for (i in 1:p) {
Ctmp[, , i, , ] <- C[, , -i, -i]
}
}
}
if ("alpha" %in% estimates) {
res$samp$Bayes_alpha <- coda::mcmc(apply(C, MARGIN = c(1, 2), applyalpha, callback))
int <- coda::HPDinterval(coda::mcmc(as.vector(res$samp$Bayes_alpha)), prob = interval)
res$cred$low$Bayes_alpha <- int[1]
res$cred$up$Bayes_alpha <- int[2]
res$est$Bayes_alpha <- mean(res$samp$Bayes_alpha)
if (item.dropped) {
res$ifitem$samp$alpha <- (apply(Ctmp, c(1, 2, 3), applyalpha, callback))
res$ifitem$est$alpha <- apply(res$ifitem$samp$alpha, 3, mean)
res$ifitem$cred$alpha <- coda::HPDinterval(coda::mcmc(apply(res$ifitem$samp$alpha, 3, as.vector)),
prob = interval)
}
}
if ("lambda2" %in% estimates) {
res$samp$Bayes_lambda2 <- coda::mcmc(apply(C, MARGIN = c(1, 2), applylambda2, callback))
int <- coda::HPDinterval(coda::mcmc(as.vector(res$samp$Bayes_lambda2)), prob = interval)
res$cred$low$Bayes_lambda2 <- int[1]
res$cred$up$Bayes_lambda2 <- int[2]
res$est$Bayes_lambda2 <- mean(res$samp$Bayes_lambda2)
if (item.dropped) {
res$ifitem$samp$lambda2 <- apply(Ctmp, c(1, 2, 3), applylambda2, callback)
res$ifitem$est$lambda2 <- apply(res$ifitem$samp$lambda2, 3, mean)
res$ifitem$cred$lambda2 <- coda::HPDinterval(coda::mcmc(apply(res$ifitem$samp$lambda2, 3, as.vector)),
prob = interval)
}
}
if ("lambda4" %in% estimates) {
res$samp$Bayes_lambda4 <- coda::mcmc(apply(C, MARGIN = c(1, 2), applylambda4NoCpp, callback))
int <- coda::HPDinterval(coda::mcmc(as.vector(res$samp$Bayes_lambda4)), prob = interval)
res$cred$low$Bayes_lambda4 <- int[1]
res$cred$up$Bayes_lambda4 <- int[2]
res$est$Bayes_lambda4 <- mean(res$samp$Bayes_lambda4)
if (item.dropped) {
res$ifitem$samp$lambda4 <- (apply(Ctmp, c(1, 2, 3), applylambda4NoCpp, callback))
res$ifitem$est$lambda4 <- apply(res$ifitem$samp$lambda4, 3, mean)
res$ifitem$cred$lambda4 <- coda::HPDinterval(coda::mcmc(apply(res$ifitem$samp$lambda4, 3, as.vector)),
prob = interval)
}
}
if ("lambda6" %in% estimates) {
res$samp$Bayes_lambda6 <- coda::mcmc(apply(C, MARGIN = c(1, 2), applylambda6, callback))
int <- coda::HPDinterval(coda::mcmc(as.vector(res$samp$Bayes_lambda6)), prob = interval)
res$cred$low$Bayes_lambda6 <- int[1]
res$cred$up$Bayes_lambda6 <- int[2]
res$est$Bayes_lambda6 <- mean(res$samp$Bayes_lambda6)
if (item.dropped) {
res$ifitem$samp$lambda6 <- (apply(Ctmp, c(1, 2, 3), applylambda6, callback))
res$ifitem$est$lambda6 <- apply(res$ifitem$samp$lambda6, 3, mean)
res$ifitem$cred$lambda6 <- coda::HPDinterval(coda::mcmc(apply(res$ifitem$samp$lambda6, 3, as.vector)),
prob = interval)
}
}
if ("glb" %in% estimates) {
res$samp$Bayes_glb <- coda::mcmc(t(apply(C, c(1), glbOnArrayCustom, callback)))
if (sum(is.na(res$samp$Bayes_glb) > 0)) {
int <- c(NA, NA)
} else {
int <- coda::HPDinterval(coda::mcmc(as.vector(res$samp$Bayes_glb)), prob = interval)
}
res$cred$low$Bayes_glb <- int[1]
res$cred$up$Bayes_glb <- int[2]
res$est$Bayes_glb <- mean(res$samp$Bayes_glb)
if (item.dropped){
res$ifitem$samp$glb <- aperm(apply(Ctmp, c(1, 3), glbOnArrayCustom, callback), c(2, 1, 3))
res$ifitem$est$glb <- apply(res$ifitem$samp$glb, 3, mean)
res$ifitem$cred$glb <- coda::HPDinterval(coda::mcmc(apply(res$ifitem$samp$glb, 3, as.vector)),
prob = interval)
}
}
# special case omega -----------------------------------------------------------------
if ("omega" %in% estimates) {
om_samp <- omegaSampler(data, n.iter, n.burnin, thin, n.chains, pairwise, callback, a0, b0, m0)
res$samp$Bayes_omega <- om_samp$omega
res$data_mis_samp_fm <- om_samp$dat_mis_samp_fm
res$loadings <- apply(om_samp$lambda, 3, as.vector)
res$resid_var <- apply(om_samp$psi, 3, as.vector)
res$f_var <- c(om_samp$phi)
int <- coda::HPDinterval(coda::mcmc(as.vector(res$samp$Bayes_omega)), prob = interval)
res$cred$low$Bayes_omega <- int[1]
res$cred$up$Bayes_omega <- int[2]
res$est$Bayes_omega <- mean(res$samp$Bayes_omega)
if (item.dropped) {
om_samp_ifitem <- array(0, c(n.chains, length(seq(1, n.iter - n.burnin, thin)), p))
for (i in 1:p) {
tmp <- data[, -i]
om_samp_ifitem[, , i] <- omegaSampler(tmp, n.iter, n.burnin, thin, n.chains,
pairwise, callback, a0, b0, m0)$omega
}
res$ifitem$samp$omega <- om_samp_ifitem
res$ifitem$est$omega <- apply(om_samp_ifitem, 3, mean)
res$ifitem$cred$omega <- coda::HPDinterval(coda::mcmc(apply(res$ifitem$samp$omega, 3, as.vector)),
prob = interval)
}
}
return(res)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.