R/rheteroLinearIndepMetrop.R

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

#' @title Distributed Independence Metropolis-Hastings Algorithm for Draws From Multivariate Normal Distribution
#' 
#' @description
#' \code{rheteroLinearIndepMetrop} implements an Independence Metropolis-Hastings algorithm with a Gibbs sampler.
#' 
#' @param Data A list of data partitions where each partition includes: `regdata` - A \eqn{nreg} size list of multinomial regdata, and optional `Z`- \eqn{nreg \times nz} matrix of unit characteristics.
#' @param betadraws A list of betadraws returned from either \code{rhierLinearMixtureParallel} or \code{rhierLinearDPParallel}
#' @param Mcmc A list with one required parameter: `R`-number of iterations, and optional parameters: `keep` and `nprint`.
#' @param Prior A list with one required parameter: `ncomp`, and optional parameters: `deltabar`, `Ad`, `mubar`, `Amu`, `nu`, `V`, `nu.e`, and `ssq`.
#' 
#' @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.
#' }
#' 
#' @details
#' \subsection{Model and Priors}{
#'     \code{nreg} regression equations with \code{nvar} as the number of \eqn{X} vars in each equation \cr
#'     \eqn{y_i = X_i\beta_i + e_i} with \eqn{e_i} \eqn{\sim}{~} \eqn{N(0, \tau_i)}
#'
#'     \eqn{\tau_i} \eqn{\sim}{~} \eqn{nu.e*ssq_i/\chi^2_{nu.e}} where  \eqn{\tau_i} is the variance of \eqn{e_i}\cr
#'     \eqn{B = Z\Delta + U} or \eqn{\beta_i = \Delta' Z[i,]' + u_i} \cr
#'     \eqn{\Delta} is an \eqn{nz \times nvar} matrix \cr
#'
#'     \eqn{u_i} \eqn{\sim}{~} \eqn{N(\mu_{ind}, \Sigma_{ind})}\cr
#'     
#'     \eqn{delta = vec(\Delta)} \eqn{\sim}{~} \eqn{N(deltabar, A_d^{-1})}\cr
#'     \eqn{\mu_j} \eqn{\sim}{~} \eqn{N(mubar, \Sigma_j(x) Amu^{-1})}\cr
#'     \eqn{\Sigma_j} \eqn{\sim}{~} \eqn{IW(nu, V)} \cr
#'    
#'     \eqn{\theta_i = (\mu_i, \Sigma_i)} \eqn{\sim}{~} \eqn{DP(G_0(\lambda), alpha)}\cr
#'     
#'     \eqn{G_0(\lambda):}\cr
#'     \eqn{\mu_i | \Sigma_i} \eqn{\sim}{~} \eqn{N(0, \Sigma_i (x) a^{-1})}\cr
#'     \eqn{\Sigma_i} \eqn{\sim}{~} \eqn{IW(nu, nu*v*I)}\cr
#'     \eqn{delta = vec(\Delta)} \eqn{\sim}{~} \eqn{N(deltabar, A_d^{-1})}\cr
#'     
#'      \eqn{\lambda(a, nu, v):}\cr
#'      \eqn{a} \eqn{\sim}{~} uniform[alim[1], alimb[2]]\cr
#'      \eqn{nu} \eqn{\sim}{~}  dim(data)-1 + exp(z) \cr
#'      \eqn{z} \eqn{\sim}{~}  uniform[dim(data)-1+nulim[1], nulim[2]]\cr
#'      \eqn{v} \eqn{\sim}{~} uniform[vlim[1], vlim[2]]
#'     
#'     \eqn{Z} should \emph{not} include an intercept and should be centered for ease of interpretation. 
#'     The mean of each of the \code{nreg} \eqn{\beta}s is the mean of the normal mixture. 
#'     Use \code{summary()} to compute this mean from the \code{compdraw} output.
#'     
#'     The prior on \eqn{\Sigma_i} is parameterized such that \eqn{mode(\Sigma) = nu/(nu+2) vI}. The support of nu enforces a non-degenerate IW density; \eqn{nulim[1] > 0}
#'     
#'     The default choices of alim, nulim, and vlim determine the location and approximate size of candidate "atoms" or possible normal components. The defaults are sensible given a reasonable scaling of the X variables. You want to insure that alim is set for a wide enough range of values (remember a is a precision parameter) and the v is big enough to propose Sigma matrices wide enough to cover the data range.
#'
#'     A careful analyst should look at the posterior distribution of a, nu, v to make sure that the support is set correctly in alim, nulim, vlim. In other words, if we see the posterior bunched up at one end of these support ranges, we should widen the range and rerun.
#'
#'     If you want to force the procedure to use many small atoms, then set nulim to consider only large values and set vlim to consider only small scaling constants. Set alphamax to a large number. This will create a very "lumpy" density estimate somewhat like the classical Kernel density estimates. Of course, this is not advised if you have a prior belief that densities are relatively smooth.
#' }
#' 
#' \subsection{Argument Details}{
#' \emph{\code{Data  = list(regdata, Z)} [\code{Z} optional]}
#' \tabular{ll}{
#'     \code{regdata:        } \tab A \eqn{nreg/s_shard} size list of regdata \cr
#'     \code{regdata[[i]]$X: } \tab \eqn{n_i \times nvar} design matrix for equation \eqn{i}  \cr
#'     \code{regdata[[i]]$y: } \tab \eqn{n_i \times 1} vector of observations for equation \eqn{i} \cr
#'     \code{Z:              } \tab A list of s partitions where each partition include \eqn{(nreg/s_shard) \times nz} matrix of unit characteristics 
#'     }
#'     
#'  
#' \emph{\code{betadraw:}} A matrix with \eqn{R} rows and \eqn{nvar} columns of beta draws.
#' 
#'     
#' \emph{\code{Prior = list(deltabar, Ad, Prioralphalist, lambda_hyper, nu, V, nu_e, mubar, Amu, ssq, ncomp)} [all but \code{ncomp} are optional]}
#' \tabular{ll}{
#'     \code{deltabar:       } \tab \eqn{(nz \times nvar) \times 1} vector of prior means (def: 0) \cr
#'     \code{Ad:             } \tab prior precision matrix for vec(D) (def: 0.01*I) \cr
#'     \code{mubar:          } \tab \eqn{nvar \times 1} prior mean vector for normal component mean (def: 0) \cr
#'     \code{Amu:            } \tab prior precision for normal component mean (def: 0.01) \cr
#'     \code{nu.e:             } \tab d.f. parameter for regression error variance prior (def: 3) \cr
#'     \code{V:              } \tab PDS location parameter for IW prior on normal component Sigma (def: nu*I) \cr
#'     \code{ssq:            } \tab scale parameter for regression error variance prior (def: \code{var(y_i)}) \cr
#'     \code{ncomp:          } \tab number of components used in normal mixture \cr
#'     }
#'  
#' \emph{\code{Mcmc  = list(R, keep, nprint)} [only \code{R} required]}
#' \tabular{ll}{
#'     \code{R:              } \tab number of MCMC draws \cr
#'     \code{keep:           } \tab MCMC thinning parameter -- keep every \code{keep}th draw (def: 1) \cr
#'     \code{nprint:         } \tab print the estimated time remaining for every \code{nprint}'th draw (def: 100, set to 0 for no print) \cr
#'     } 
#' }
#' 
#' @seealso 
#' \code{\link{rhierLinearMixtureParallel}},
#' \code{\link{rheteroMnlIndepMetrop}}
#' 
#' @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
#' 
#'######### Single Component with rhierLinearMixtureParallel########
#'R = 500
#'
#'set.seed(66)
#'nreg=1000
#'nobs=5 #number of observations
#'nvar=3 #columns
#'nz=2
#'
#'Z=matrix(runif(nreg*nz),ncol=nz) 
#'Z=t(t(Z)-apply(Z,2,mean))
#'
#'Delta=matrix(c(1,-1,2,0,1,0),ncol=nz)
#'tau0=.1
#'iota=c(rep(1,nobs)) 
#'
#'#Default
#'tcomps=NULL
#'a=matrix(c(1,0,0,0.5773503,1.1547005,0,-0.4082483,0.4082483,1.2247449),ncol=3) 
#'tcomps[[1]]=list(mu=c(0,-1,-2),rooti=a) 
#'tpvec=c(1)                               
#'
#'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)
#'}
#'
#'
#'Data1=list(list(regdata=regdata,Z=Z))
#'s = 1
#'Data2=scalablebayesm::partition_data(Data1,s=s)
#'
#'Prior1=list(ncomp=1)
#'Mcmc1=list(R=R,keep=1)
#'
#'set.seed(1)
#'out2 = parallel::mclapply(Data2, FUN = rhierLinearMixtureParallel, 
#'Prior = Prior1, Mcmc = Mcmc1,
#' mc.cores = s, mc.set.seed = FALSE)
#' 
#' betadraws = parallel::mclapply(out2,FUN=drawPosteriorParallel,Z=Z, 
#' Prior = Prior1, Mcmc = Mcmc1, mc.cores=s,mc.set.seed = FALSE)
#' betadraws = combine_draws(betadraws, R)
#' 
#' out_indep = parallel::mclapply(Data2, FUN=rheteroLinearIndepMetrop, 
#' betadraws = betadraws, Mcmc = Mcmc1, Prior = Prior1, mc.cores = s, mc.set.seed = FALSE)
#'
#'
#'######### Multiple Components with rhierLinearMixtureParallel########
#'R = 500
#'
#'set.seed(66)
#'nreg=1000
#'nobs=5 #number of observations
#'nvar=3 #columns
#'nz=2
#'
#'Z=matrix(runif(nreg*nz),ncol=nz) 
#'Z=t(t(Z)-apply(Z,2,mean))
#'
#'Delta=matrix(c(1,-1,2,0,1,0),ncol=nz)
#'tau0=.1
#'iota=c(rep(1,nobs)) 
#'
#'#Default
#'tcomps=NULL
#'a=matrix(c(1,0,0,0.5773503,1.1547005,0,-0.4082483,0.4082483,1.2247449),ncol=3) 
#'tcomps[[1]]=list(mu=c(0,-1,-2),rooti=a) 
#'tcomps[[2]]=list(mu=c(0,-1,-2)*2,rooti=a)
#'tcomps[[3]]=list(mu=c(0,-1,-2)*4,rooti=a)
#'tpvec=c(.4,.2,.4)                                   
#'
#'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)
#'}
#'
#'
#'Data1=list(list(regdata=regdata,Z=Z))
#'s = 1
#'Data2=scalablebayesm::partition_data(Data1, s=s)
#'
#'Prior1=list(ncomp=3)
#'Mcmc1=list(R=R,keep=1)
#'
#'set.seed(1)
#'out2 = parallel::mclapply(Data2, FUN = rhierLinearMixtureParallel, Prior = Prior1, Mcmc = Mcmc1,
#' mc.cores = s, mc.set.seed = FALSE)
#' 
#' betadraws = parallel::mclapply(out2,FUN=drawPosteriorParallel,Z=Z, 
#' Prior = Prior1, Mcmc = Mcmc1, mc.cores=s,mc.set.seed = TRUE)
#' betadraws = combine_draws(betadraws, R)
#' 
#' out_indep = parallel::mclapply(Data2, FUN=rheteroLinearIndepMetrop, 
#' betadraws = betadraws, Mcmc = Mcmc1, Prior = Prior1, mc.cores = s, mc.set.seed = TRUE)
#'
#'
#'
#' @export
#' 
rheteroLinearIndepMetrop=function(Data, betadraws, Mcmc, Prior){
  
  nvar = ncol(Data$regdata[[1]]$X)
  
  getvar=function(l) { 
    v=var(l$y)
    if(is.na(v)) return(1)
    if(v>0) return (v) else return (1)}
  
  if(missing(Prior)) 
  {stop("Requires Prior list argument (at least ncomp)", .call=FALSE)} 
  
  if(is.null(Prior$nu)) {
    nu=nvar+ScalableBayesmConstant.nuInc
  } else {
    nu=Prior$nu
  }
  
  if(is.null(Prior$ssq)) {
    ssq=sapply(Data$regdata,getvar)
  } else {
    ssq=Prior$ssq
  }
  
  
  drawsindep=rheteroLinearIndepMetrop_rcpp_loop(Data,betadraws,Mcmc, nu, ssq)
  return(drawsindep)
  #return(list(betadraw = drawsindep))
  
}




#attr(betadraw,"dimnames")<-NULL
#return(list(betadraw=betadraw))