rmnlIndepMetrop: MCMC Algorithm for Multinomial Logit Model

In bayesm: Bayesian Inference for Marketing/Micro-Econometrics

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, [e-mail].

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

bayesm documentation built on Sept. 24, 2023, 1:07 a.m.