bayesm (version 3.1-6)

rmnlIndepMetrop: MCMC Algorithm for Multinomial Logit Model

Description

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

Usage

rmnlIndepMetrop(Data, Prior, Mcmc)

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

Arguments

Data

list(y, X, p)

Prior

list(A, betabar)

Mcmc

list(R, keep, nprint, nu)

Author

Peter Rossi, Anderson School, UCLA, .

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)

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