rmnlIndepMetrop function

MCMC Algorithm for Multinomial Logit Model

rmnlIndepMetrop implements Independence Metropolis algorithm for the multinomial logit (MNL) model.

rmnlIndepMetrop(Data, Prior, Mcmc)

Arguments

• Data: list(y, X, p)
• Prior: list(A, betabar)
• Mcmc: list(R, keep, nprint, nu)

Details

Model and Priors

y $~$ MNL(X, $\beta$)

$Pr(y=j) = exp(x_j'\beta)/\sum_k{e^{x_k'\beta}}$

$\beta$ $~$ $N(betabar, A^{-1})$

Argument Details

Data = list(y, X, p)

y:	$n x 1$ vector of multinomial outcomes (1, , p)
X:	$n*p x k$ matrix
p:	number of alternatives

Prior = list(A, betabar) [optional]

A:	$k x k$ prior precision matrix (def: 0.01*I)
betabar:	$k x 1$ prior mean (def: 0)

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

R:	number of MCMC draws
keep:	MCMC thinning parameter -- keep every keep th draw (def: 1)
nprint:	print the estimated time remaining for every nprint 'th draw (def: 100, set to 0 for no print)
nu:	d.f. parameter for independent t density (def: 6)

Returns

A list containing: - betadraw: $R/keep x k$ matrix of beta draws

• loglike: $R/keep x 1$ vector of log-likelihood values evaluated at each draw

• acceptr: acceptance rate of Metropolis draws

Author(s)

Peter Rossi, Anderson School, UCLA, perossichi@gmail.com .

References

For further discussion, see Chapter 3, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

See Also

rhierMnlRwMixture

Examples

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

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))
}

n = 200
p = 3
beta = c(1, -1, 1.5, 0.5)

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)

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