MCMC Algorithm for Multinomial Logit Model
rmnlIndepMetrop implements Independence Metropolis algorithm for the multinomial logit (MNL) model.
rmnlIndepMetrop(Data, Prior, Mcmc)
Data |
list(y, X, p) |
Prior |
list(A, betabar) |
Mcmc |
list(R, keep, nprint, nu) |
y ~ MNL(X, β)
Pr(y=j) = exp(x_j'β)/∑_k{e^{x_k'β}}
β ~ N(betabar, A^{-1})
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) |
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 |
Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.
For further discussion, see Chapter 3, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.
http://www.perossi.org/home/bsm-1
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)}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.