rmnpGibbs {bayesm}

Gibbs Sampler for Multinomial Probit

Description

rmnpGibbs implements the McCulloch/Rossi Gibbs Sampler for the multinomial probit model.

Usage

rmnpGibbs(Data, Prior, Mcmc)

Arguments

Data

list(y, X, p)

Prior

list(betabar, A, nu, V)

Mcmc

list(R, keep, nprint, beta0, sigma0)

Details

Model and Priors

w_i = X_i\beta + e with e \sim N(0, \Sigma). Note: w_i and e are (p-1) x 1.
y_i = j if w_{ij} > max(0,w_{i,-j}) for j=1, \ldots, p-1 and where w_{i,-j} means elements of w_i other than the jth.
y_i = p, if all w_i < 0

\beta is not identified. However, \beta/sqrt(\sigma_{11}) and \Sigma/\sigma_{11} are identified. See reference or example below for details.

\beta \sim N(betabar,A^{-1})
\Sigma \sim IW(nu,V)

Argument Details

Data = list(y, X, p)

y:	n x 1 vector of multinomial outcomes (1, ..., p)
X:	n*(p-1) x k design matrix. To make X matrix use createX with DIFF=TRUE
p:	number of alternatives

Prior = list(betabar, A, nu, V) [optional]

betabar:	k x 1 prior mean (def: 0)
A:	k x k prior precision matrix (def: 0.01*I)
nu:	d.f. parameter for Inverted Wishart prior (def: (p-1)+3)
V:	PDS location parameter for Inverted Wishart prior (def: nu*I)

Mcmc = list(R, keep, nprint, beta0, sigma0) [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)
beta0:	initial value for beta (def: 0)
sigma0:	initial value for sigma (def: I)

Value

A list containing:

betadraw

R/keep x k matrix of betadraws

sigmadraw

R/keep x (p-1)*(p-1) matrix of sigma draws – each row is the vector form of sigma

Author(s)

Peter Rossi, Anderson School, UCLA, [e-mail].

References

For further discussion, see Chapter 4, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

See Also

rmvpGibbs

Examples

if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10}
set.seed(66)

simmnp = function(X, p, n, beta, sigma) {
  indmax = function(x) {which(max(x)==x)}
  Xbeta = X%*%beta
  w = as.vector(crossprod(chol(sigma),matrix(rnorm((p-1)*n),ncol=n))) + Xbeta
  w = matrix(w, ncol=(p-1), byrow=TRUE)
  maxw = apply(w, 1, max)
  y = apply(w, 1, indmax)
  y = ifelse(maxw < 0, p, y)
  return(list(y=y, X=X, beta=beta, sigma=sigma))
}

p = 3
n = 500
beta = c(-1,1,1,2)
Sigma = matrix(c(1, 0.5, 0.5, 1), ncol=2)
k = length(beta)
X1 = matrix(runif(n*p,min=0,max=2),ncol=p)
X2 = matrix(runif(n*p,min=0,max=2),ncol=p)
X = createX(p, na=2, nd=NULL, Xa=cbind(X1,X2), Xd=NULL, DIFF=TRUE, base=p)

simout = simmnp(X,p,500,beta,Sigma)

Data1 = list(p=p, y=simout$y, X=simout$X)
Mcmc1 = list(R=R, keep=1)

out = rmnpGibbs(Data=Data1, Mcmc=Mcmc1)

cat(" Summary of Betadraws ", fill=TRUE)
betatilde = out$betadraw / sqrt(out$sigmadraw[,1])
attributes(betatilde)$class = "bayesm.mat"
summary(betatilde, tvalues=beta)

cat(" Summary of Sigmadraws ", fill=TRUE)
sigmadraw = out$sigmadraw / out$sigmadraw[,1]
attributes(sigmadraw)$class = "bayesm.var"
summary(sigmadraw, tvalues=as.vector(Sigma[upper.tri(Sigma,diag=TRUE)]))

## plotting examples
if(0){plot(betatilde,tvalues=beta)}