R/rheteroMnlIndepMetrop.R

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

#' @title Independence Metropolis-Hastings Algorithm for Draws From Multinomial Distribution

#' @description
#' \code{rheteroMnlIndepMetrop} implements an Independence Metropolis algorithm with a Gibbs sampler.
#' 
#' @param Data A list of s partitions where each partition includes: `p`- number of choice alternatives, `lgtdata` - An \eqn{nlgt} size list of multinomial logistic data, and optional `Z`- matrix of unit characteristics.
#' @param draws A list of draws returned from either \code{rhierMnlRwMixtureParallel}.
#' @param Mcmc A list with one required parameter: `R`-number of iterations, and optional parameters: `keep` and `nprint`.
#' 
#' @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}{
#' \eqn{y_i} \eqn{\sim}{~} \eqn{MNL(X_i,\beta_i)} with \eqn{i = 1, \ldots,} length(lgtdata) 
#' and where \eqn{\beta_i} is \eqn{1 \times nvar}
#' 
#' \eqn{\beta_i} = \eqn{Z\Delta}[i,] + \eqn{u_i} \cr
#' Note:  Z\eqn{\Delta} is the matrix Z \eqn{ \times \Delta} and [i,] refers to \eqn{i}th row of this product \cr
#' Delta is an \eqn{nz \times nvar} array 
#' 
#' \eqn{u_i} \eqn{\sim}{~} \eqn{N(\mu_{ind},\Sigma_{ind})} with \eqn{ind} \eqn{\sim}{~} multinomial(pvec)
#' 
#' \eqn{pvec}                \eqn{\sim}{~} dirichlet(a) \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)}
#' 
#' Note: \eqn{Z} should NOT include an intercept and is centered for ease of interpretation. 
#' The mean of each of the \code{nlgt} \eqn{\beta}s is the mean of the normal mixture. 
#' Use \code{summary()} to compute this mean from the \code{compdraw} output.\cr
#' 
#' Be careful in assessing prior parameter \code{Amu}: 0.01 is too small for many applications. 
#' See chapter 5 of Rossi et al for full discussion.
#' }
#' \subsection{Argument Details}{
#' \emph{\code{Data  = list(lgtdata, Z, p)} [\code{Z} optional]}
#' \tabular{ll}{
#' \code{lgtdata:        } \tab A \eqn{nlgt/shards} size list of multinominal lgtdata \cr
#' \code{lgtdata[[i]]$y: } \tab \eqn{n_i \times 1} vector of multinomial outcomes (1, \ldots, p) for \eqn{i}th unit\cr
#' \code{lgtdata[[i]]$X: } \tab \eqn{(n_i \times p) \times nvar} design matrix for \eqn{i}th unit \cr
#' \code{Z:              } \tab A list of s partitions where each partition include \eqn{(nlgt/number of shards) \times nz} matrix of unit characteristics\cr 
#' \code{p:              } \tab number of choice alternatives 
#' }
#'  
#' \emph{\code{draws:}} A matrix with \eqn{R} rows and \eqn{nlgt} columns of beta draws.
#' 
#' \emph{\code{Mcmc  = list(R, keep, nprint, s, w)} [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
#' \code{s:              } \tab scaling parameter for RW Metropolis (def: 2.93/\code{sqrt(nvar)}) \cr
#' \code{w:              } \tab fractional likelihood weighting parameter (def: 0.1)
#' }
#' }
#' 
#' @seealso 
#' \code{\link{rhierMnlRwMixtureParallel}},
#' \code{\link{rheteroLinearIndepMetrop}}
#' 
#' @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{
#' R = 500
#' 
#' ######### Single Component with rhierMnlRwMixtureParallel########
#' ##parameters
#' p=3 # num of choice alterns                            
#' ncoef=3  
#' nlgt=1000                          
#' nz=2
#' Z=matrix(runif(nz*nlgt),ncol=nz)
#' Z=t(t(Z)-apply(Z,2,mean)) # demean Z
#' 
#' ncomp=1 # no of mixture components
#' Delta=matrix(c(1,0,1,0,1,2),ncol=2)
#' comps=NULL
#' comps[[1]]=list(mu=c(0,2,1),rooti=diag(rep(1,3)))
#' pvec=c(1)
#' 
#' simmnlwX= function(n,X,beta){
#'   k=length(beta)
#'   Xbeta=X %*% beta
#'   j=nrow(Xbeta)/n
#'   Xbeta=matrix(Xbeta,byrow=TRUE,ncol=j)
#'   Prob=exp(Xbeta)
#'   iota=c(rep(1,j))
#'   denom=Prob %*% iota
#'   Prob=Prob/as.vector(denom)
#'   y=vector("double",n)
#'   ind=1:j
#'   for (i in 1:n) { 
#'     yvec = rmultinom(1, 1, Prob[i,])
#'     y[i] = ind%*%yvec
#'   }
#'   return(list(y=y,X=X,beta=beta,prob=Prob))
#' }
#' 
#' ## simulate data
#' simlgtdata=NULL
#' ni=rep(5,nlgt) 
#' for (i in 1:nlgt) 
#' {
#'   if (is.null(Z))
#'   {
#'     betai=as.vector(bayesm::rmixture(1,pvec,comps)$x)
#'   } else {
#'     betai=Delta %*% Z[i,]+as.vector(bayesm::rmixture(1,pvec,comps)$x)
#'   }
#'   Xa=matrix(runif(ni[i]*p,min=-1.5,max=0),ncol=p)
#'   X=bayesm::createX(p,na=1,nd=NULL,Xa=Xa,Xd=NULL,base=1)
#'   outa=simmnlwX(ni[i],X,betai)
#'   simlgtdata[[i]]=list(y=outa$y,X=X,beta=betai)
#' }
#' 
#' ## set MCMC parameters
#' Prior1=list(ncomp=ncomp) 
#' keep=1
#' Mcmc1=list(R=R,keep=keep) 
#' Data1=list(list(p=p,lgtdata=simlgtdata,Z=Z))
#' s = 1
#' Data2 = partition_data(Data1, s=s)
#' 
#' out2 = parallel::mclapply(Data2, FUN = rhierMnlRwMixtureParallel, 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=rheteroMnlIndepMetrop, draws = betadraws, 
#' Mcmc = Mcmc1, mc.cores = s, mc.set.seed = FALSE)
#'
#' 
#' ######### Multiple Components with rhierMnlRwMixtureParallel########
#' ##parameters
#' R=500
#' p=3 # num of choice alterns                            
#' ncoef=3  
#' nlgt=1000  # num of cross sectional units                         
#' nz=2
#' Z=matrix(runif(nz*nlgt),ncol=nz)
#' Z=t(t(Z)-apply(Z,2,mean)) # demean Z
#' 
#' ncomp=3
#' Delta=matrix(c(1,0,1,0,1,2),ncol=2) 
#' 
#' comps=NULL 
#' comps[[1]]=list(mu=c(0,2,1),rooti=diag(rep(1,3)))
#' comps[[2]]=list(mu=c(1,0,2),rooti=diag(rep(1,3)))
#' comps[[3]]=list(mu=c(2,1,0),rooti=diag(rep(1,3)))
#' pvec=c(.4,.2,.4)
#' 
#' simmnlwX= function(n,X,beta) {
#'   k=length(beta)
#'   Xbeta=X %*% beta
#'   j=nrow(Xbeta)/n
#'   Xbeta=matrix(Xbeta,byrow=TRUE,ncol=j)
#'   Prob=exp(Xbeta)
#'   iota=c(rep(1,j))
#'   denom=Prob %*% iota
#'   Prob=Prob/as.vector(denom)
#'   y=vector("double",n)
#'   ind=1:j
#'   for (i in 1:n) 
#'   {yvec=rmultinom(1,1,Prob[i,]); y[i]=ind %*% yvec}
#'   return(list(y=y,X=X,beta=beta,prob=Prob))
#' }
#' 
#' ## simulate data
#' simlgtdata=NULL
#' ni=rep(5,nlgt) 
#' for (i in 1:nlgt) 
#' {
#'   if (is.null(Z))
#'   {
#'     betai=as.vector(bayesm::rmixture(1,pvec,comps)$x)
#'   } else {
#'     betai=Delta %*% Z[i,]+as.vector(bayesm::rmixture(1,pvec,comps)$x)
#'   }
#'   Xa=matrix(runif(ni[i]*p,min=-1.5,max=0),ncol=p)
#'   X=bayesm::createX(p,na=1,nd=NULL,Xa=Xa,Xd=NULL,base=1)
#'   outa=simmnlwX(ni[i],X,betai) 
#'   simlgtdata[[i]]=list(y=outa$y,X=X,beta=betai)
#' }
#' 
#' ## set parameters for priors and Z
#' Prior1=list(ncomp=ncomp) 
#' keep = 1
#' Mcmc1=list(R=R,keep=keep) 
#' Data1=list(list(p=p,lgtdata=simlgtdata,Z=Z))
#' s = 1
#' Data2 = partition_data(Data1,s)
#' 
#' out2 = parallel::mclapply(Data2, FUN = rhierMnlRwMixtureParallel, 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=rheteroMnlIndepMetrop, draws = betadraws, 
#' Mcmc = Mcmc1, mc.cores = s, mc.set.seed = FALSE)
#' }
#' @export
rheteroMnlIndepMetrop = function(Data, draws, Mcmc) {
  out = rheteroMnlIndepMetrop_rcpp_loop(Data = Data, draws = draws, Mcmc = Mcmc)
  return(out)
}