bayesm (version 3.0-2)

rmnlIndepMetrop: MCMC Algorithm for Multinomial Logit Model

Description

rmnIndepMetrop implements Independence Metropolis for the MNL.

Usage

rmnlIndepMetrop(Data, Prior, Mcmc)

Arguments

Data

list(p,y,X)

Prior

list(A,betabar) optional

Mcmc

list(R,keep,nprint,nu)

Value

a list containing:

betadraw

R/keep x nvar array of beta draws

loglike

R/keep vector of loglike values for each draw

acceptr

acceptance rate of Metropolis draws

Details

Model: y \(\sim\) MNL(X,\(\beta\)). \(\Pr(y=j) = exp(x_j'\beta)/\sum_k{e^{x_k'\beta}}\).

Prior: \(\beta\) \(\sim\) \(N(betabar,A^{-1})\)

list arguments contain:

• p number of alternatives

• y nobs vector of multinomial outcomes (1,…, p)

• X nobs*p x nvar matrix

• A nvar x nvar pds prior prec matrix (def: .01I)

• betabar nvar x 1 prior mean (def: 0)

• R number of MCMC draws

• keep MCMC thinning parm: keep every keepth draw (def: 1)

• nprint print the estimated time remaining for every nprint'th draw (def: 100)

• nu degrees of freedom parameter for independence t density (def: 6)

References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch, Chapter 5. http://www.perossi.org/home/bsm-1

See Also

rhierMnlRwMixture

Examples

##

if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}

set.seed(66)
n=200; p=3; beta=c(1,-1,1.5,.5)

simmnl= function(p,n,beta) {
  #   note: create X array with 2 alt.spec vars
  k=length(beta)
  X1=matrix(runif(n*p,min=-1,max=1),ncol=p)
  X2=matrix(runif(n*p,min=-1,max=1),ncol=p)
  X=createX(p,na=2,nd=NULL,Xd=NULL,Xa=cbind(X1,X2),base=1)
  Xbeta=X%*%beta # now do probs
  p=nrow(Xbeta)/n
  Xbeta=matrix(Xbeta,byrow=TRUE,ncol=p)
  Prob=exp(Xbeta)
  iota=c(rep(1,p))
  denom=Prob%*%iota
  Prob=Prob/as.vector(denom)
  # draw y
  y=vector("double",n)
  ind=1:p
  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))
}

simout=simmnl(p,n,beta)

Data1=list(y=simout$y,X=simout$X,p=p); Mcmc1=list(R=R,keep=1)
out=rmnlIndepMetrop(Data=Data1,Mcmc=Mcmc1)

cat("Summary of beta draws",fill=TRUE)
summary(out$betadraw,tvalues=beta)

if(0){
## plotting examples
plot(out$betadraw)
}