R/drawPosteriorParallel.R

In scalablebayesm: Distributed Markov Chain Monte Carlo for Bayesian Inference in Marketing

#' Draw from Posterior Parallel Distribution
#' 
#' @description
#' \code{drawPosteriorParallel} draws from a posterior predictive distribution.
#' 
#' @param draws (list) - a list of length s where each sublist contains compdraw, 
#' @param Z (matrix) - (optional) an \eqn{nreg/s \times nz} matrix of unit characteristics
#' @param Prior (list) - (optional) a list of optional parameters 'v' and 'nu'
#' @param Mcmc (list) - a list containing 'R' and optionally 'keep'
#' 
#' @return A list containing:
#' \itemize{
#'   \item \strong{betadraw}: A matrix of size \eqn{R \times nvar} containing the drawn \code{beta} values from the Gibbs sampling procedure.
#' }
#' 
#' @author Federico Bumbaca, Leeds School of Business, University of Colorado Boulder, \email{federico.bumbaca@colorado.edu}
#' @references Bumbaca, F. (Rico), Misra, S., & Rossi, P. E. (2020). Scalable Target Marketing: Distributed Markov Chain Monte Carlo for Bayesian Hierarchical Models. Journal of Marketing Research, 57(6), 999-1018.
#' @examples
#'
#'\donttest{
#'
#'s=1
#'R=2000
#'nreg = 2000
#'nobs=5 #number of observations
#'nvar=3 #columns
#'nz=2
#'
#'Z=NULL
#'Delta=matrix(c(1,0,1,0,1,2),ncol=nz) 
#'tau0=1
#'iota=c(rep(1,nobs)) 
#'
#'## create arguments for rmixture
#'
#'#Default
#'tcomps=NULL
#'a = diag(1, nrow=3)
#'tcomps[[1]] = list(mu=c(0,-1,-2),rooti=a) 
#'tpvec = 1                             
#'ncomp=length(tcomps)
#'
#'regdata=NULL
#'betas=matrix(double(nreg*nvar),ncol=nvar) 
#'tind=double(nreg) 
#'
#'for (reg in 1:nreg) { 
#'  tempout=bayesm::rmixture(1,tpvec,tcomps)
#'  if (is.null(Z)){
#'    betas[reg,]= as.vector(tempout$x)  
#'  }else{
#'    betas[reg,]=Delta%*%Z[reg,]+as.vector(tempout$x)} 
#'  tind[reg]=tempout$z
#'  X=cbind(iota,matrix(runif(nobs*(nvar-1)),ncol=(nvar-1))) 
#'  tau=tau0*runif(1,min=0.5,max=1) 
#'  y=X%*%betas[reg,]+sqrt(tau)*rnorm(nobs)
#'  regdata[[reg]]=list(y=y,X=X,beta=betas[reg,],tau=tau) 
#'}
#'
#'Prior1=list(ncomp=ncomp) 
#'keep=1
#'Mcmc1=list(R=R,keep=keep)
#'Data1=list(list(regdata=regdata,Z=Z))
#'
#'
#'Data2 = partition_data(Data1, s)
#'
#'draws = parallel::mclapply(Data2, FUN = rhierLinearMixtureParallel, Prior = Prior1, Mcmc = Mcmc1, 
#'mc.cores = s, mc.set.seed = TRUE)
#'
#'out = parallel::mclapply(draws,FUN=drawPosteriorParallel,
#'Z=Z, Prior = Prior1, Mcmc = Mcmc1, mc.cores=s,
#'mc.set.seed = TRUE)
#'
#'}
#'
#'@author Federico Bumbaca, \email{federico.bumbaca@@colorado.edu}
#'@export

drawPosteriorParallel = function(draws, Z, Prior,Mcmc) {
  
  s = length(draws)
  nvar = length(draws$compdraw[[1]][[1]]$mu)
  
  if(is.null(Prior$nu)) {
    nu=nvar+ScalableBayesmConstant.nuInc
  } else {
    nu=Prior$nu
  }
  if(is.null(Prior$v)) {
    v=nu*diag(nvar)
  } else {
    v=Prior$v
  }
  r = Mcmc$R
  post_burn_in = length(draws$compdraw)
  drawdelta=TRUE
  if(is.null(Z)){
    Z=matrix(0,nrow=1,ncol=1)
    draws$Deltadraw=matrix(0,nrow=1,ncol=1)
    drawdelta=FALSE
  }else{
    Z=Z
  }
  if(is.null(Mcmc$keep)) {
    keep=ScalableBayesmConstant.keep
  } else {
    keep=Mcmc$keep
  }
  betaDraws = drawPosteriorParallel_cpp(draws$compdraw, draws$probdraw, draws$Deltadraw, v, r, s, post_burn_in, keep, Z, drawdelta)
  return(betaDraws)
}