R/get_coda.R In bfa: Bayesian Factor Analysis

Documented in get_coda

```#' Get coda object
#'
#' Returns an \code{mcmc} object for use in \code{coda} functions for
#' convergence diagnostics, HPD intervals, etc.
#'
#' @param model a \code{bfa} model
#' @param scores Return samples of the factor scores? (default is FALSE)
#' @param scale Return factor loadings on the correlation scale? (default is TRUE for mixed/copula models, FALSE for Gaussian models)
#' @param positive Post-process to enforce positivity constraints
#' @return An \code{mcmc} object
#' @export

k = model\$K; p = model\$P; n=model\$N
nsim = model\$nsim; thin = model\$thin; nburn = model\$nburn
samples = NULL
kn=c()
for (i in 1:k) kn = c(kn, rep(i, p))

if (scale) {
if (dim(pl)[2]>1) {
for (i in 1:dim(pl)[3]) {
pl[,,i] = pl[,,i]/sqrt(model\$post.sigma2[i,]+rowSums(pl[,,i]^2))
}
} else {
for (i in 1:dim(pl)[3]) {
pl[,,i] = pl[,,i]/sqrt(model\$post.sigma2[i,]+pl[,,i]^2)
}
}
if (positive) {
for (i in 1:k) {
}
}
}
dim(pl) = c(1,p*k, nsim/thin)
pl = pl[1,,]
samples = cbind(samples, t(pl))
}

kn=c()
for (i in 1:k) kn = c(kn, rep(i, n))
scorenames = paste(model\$obslabel, kn)
if(scores) {
ps = model\$post.scores
dim(ps) = c(1,n*k, nsim/thin)
ps = ps[1,,]
rownames(ps) = scorenames
samples = cbind(samples, t(ps))
}

codaobj = mcmc(samples, nburn+1, nsim+nburn, thin=thin)
}
```

Try the bfa package in your browser

Any scripts or data that you put into this service are public.

bfa documentation built on May 29, 2017, 1:38 p.m.