bayesm (version 3.1-7)

rmnpGibbs: Gibbs Sampler for Multinomial Probit

Description

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

Usage

rmnpGibbs(Data, Prior, Mcmc)

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

Arguments

Data

list(y, X, p)

Prior

list(betabar, A, nu, V)

Mcmc

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

Author

Peter Rossi, Anderson School, UCLA, [email protected].

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 \(j\)th.
\(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)

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