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 ~ MNL(X, β)
Pr(y=j) = exp(x_j'β)/∑_k{e^{x_k'β}}

β ~ 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.
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)

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-4 Index]