Rstan output structure -- coerce samples back to 'rich' structure (arrays, matrices, etc)

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)
}