rheteroMnlIndepMetrop: Independence Metropolis-Hastings Algorithm for Draws From...

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

Independence Metropolis-Hastings Algorithm for Draws From Multinomial Distribution

Description

rheteroMnlIndepMetrop implements an Independence Metropolis algorithm with a Gibbs sampler.

Usage

rheteroMnlIndepMetrop(Data, draws, Mcmc)

Arguments

Data	A list of s partitions where each partition includes: 'p'- number of choice alternatives, 'lgtdata' - An nlgt size list of multinomial logistic data, and optional 'Z'- matrix of unit characteristics.
draws	A list of draws returned from either rhierMnlRwMixtureParallel.
Mcmc	A list with one required parameter: 'R'-number of iterations, and optional parameters: 'keep' and 'nprint'.

Details

Model and Priors

y_i \sim MNL(X_i,\beta_i) with i = 1, \ldots, length(lgtdata) and where \beta_i is 1 \times nvar

\beta_i = Z\Delta[i,] + u_i
Note: Z\Delta is the matrix Z \times \Delta and [i,] refers to ith row of this product
Delta is an nz \times nvar array

u_i \sim N(\mu_{ind},\Sigma_{ind}) with ind \sim multinomial(pvec)

pvec \sim dirichlet(a)
delta = vec(\Delta) \sim N(deltabar, A_d^{-1})
\mu_j \sim N(mubar, \Sigma_j (x) Amu^{-1})
\Sigma_j \sim IW(nu, V)

Note: Z should NOT include an intercept and is centered for ease of interpretation. The mean of each of the nlgt \betas is the mean of the normal mixture. Use summary() to compute this mean from the compdraw output.

Be careful in assessing prior parameter Amu: 0.01 is too small for many applications. See chapter 5 of Rossi et al for full discussion.

Argument Details

Data = list(lgtdata, Z, p) [Z optional]

lgtdata:	A nlgt/shards size list of multinominal lgtdata
lgtdata[[i]]$y:	n_i \times 1 vector of multinomial outcomes (1, ..., p) for ith unit
lgtdata[[i]]$X:	(n_i \times p) \times nvar design matrix for ith unit
Z:	A list of s partitions where each partition include (nlgt/number of shards) \times nz matrix of unit characteristics
p:	number of choice alternatives

draws: A matrix with R rows and nlgt columns of beta draws.

Mcmc = list(R, keep, nprint, s, w) [only R required]

R:	number of MCMC draws
keep:	MCMC thinning parameter -- keep every keepth draw (def: 1)
nprint:	print the estimated time remaining for every nprint'th draw (def: 100, set to 0 for no print)
s:	scaling parameter for RW Metropolis (def: 2.93/sqrt(nvar))
w:	fractional likelihood weighting parameter (def: 0.1)

Value

A list containing:

• betadraw: A matrix of size R \times nvar containing the drawn beta values from the Gibbs sampling procedure.

Author(s)

Federico Bumbaca, Leeds School of Business, University of Colorado Boulder, [e-mail]

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.

See Also

rhierMnlRwMixtureParallel, rheteroLinearIndepMetrop

Examples

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)

scalablebayesm documentation built on April 3, 2025, 7:55 p.m.