rmnlIndepMetrop {bayesm} R Documentation

## MCMC Algorithm for Multinomial Logit Model

### Description

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

### Usage

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 \sim MNL(X, \beta)
Pr(y=j) = exp(x_j'\beta)/\sum_k{e^{x_k'\beta}}

\beta \sim 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 keepth 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)

### Value

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.

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


[Package bayesm version 3.1-6 Index]