I’m sure I’ve done this before, and I’m equally sure it was a kind of ugly procedure, and dependent on the constrain and unconstrain_pars functions that don’t seem to be intended for the purpose and pose their own problems. Is there a nice way to achieve this? I ask today because I’d like to compute the quantiles of a partial correlation matrix, but I only sampled the correlation matrix.
I am not sure that I understand. You did as.array
or something to get the raw draws, and you then calculated a partial correlation matrix at each iteration, and now you want to put the partial correlations into the stanfit
object? If so, I don’t see how that would accomplish anything more than calling monitor
on the array of partial correlations.
I had a correlation matrix as one part of the model that was fit, and I wanted to compute the quantiles of the partial correlation matrix. I’ve done it using the relist function based on a list skeleton that I realized was available in the @inits subobject. It was even fiddlier than expected, but perhaps I’m missing something!
#' Compute functions of matrices from samples of a stanfit object
#'
#' @param stanfit object of class stanfit.
#' @param object name of stan sub object from stanfit to use for calculations.
#' @param objectindices matrix of indices, with the number of columns matching
#' the number of dimensions of the object. 'all' computes \code{which( array(1,objdims)==1,arr.ind=TRUE)},
#' where objdims is what would be returned by dim(object) if the object existed in the R environment.
#' @param calc string containing R calculation to evaluate, with the string 'object' in place of the actual object name.
#' @param summary if FALSE, a iterations * parameters matrix is returned, if TRUE,
#' rstan::monitor is first run on the output.
#'
#' @return matrix of values of the specified interactions at each iteration.
#' @export
#'
#' @examples
#' temp<-stan_postcalc(stanfit=ctstantestfit$stanfit,
#' object='DRIFT', objectindices='all', calc='exp(object)')
stan_postcalc <-function(stanfit,object,calc='object', objectindices='all', summary=TRUE){
mc=As.mcmc.list(stanfit)
m=do.call(rbind,mc)
outdims = dim(stanfit@inits[[1]][[object]]) #complete
if(objectindices=='all') objectindices <- which( array(1,outdims)==1,arr.ind=TRUE) else {
if(ncol(as.matrix(objectindices))!= length(outdims)) stop('Number of columns of object indices must match number of dimensions in object')
}
objectindices = objectindices[,outdims!=1] #subset
outdims = outdims[outdims!=1] #subset
out=array(apply(m,1,function(x) {
object = array(array(relist(x,skeleton = stanfit@inits[[1]])[[object]],dim=outdims)[objectindices],outdims)
ret = eval(parse(text=calc))
} ),dim = c(outdims,nrow(m)) )
if(summary) {
parnames <- array(1,dim=outdims)
parnames <- which(parnames==1,arr.ind = TRUE)
parnames <- apply(parnames,1,function(x) paste0(x, collapse=', '))
parnames <- paste0(object,'[',parnames,']')
out <- array(out,dim=c(prod(outdims), nrow(mc[[1]]), length(mc)))
out <- aperm(out,c(2,3,1))
dimnames(out)[[3]] <- parnames
out <- monitor(out,warmup=0,digits_summary = 2)
}
return(out)
}