bayesm (version 1.1-2)

rmnpGibbs: 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(p, y, X)

Prior

list(betabar,A,nu,V) (optional)

Mcmc

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

Value

• a list containing:
• betadraw
  R/keep x k array of betadraws
• sigmadraw
  R/keep x (p-1)*(p-1) array of sigma draws -- each row is in vector form

concept

• bayes
• multinomial probit
• MCMC
• Gibbs Sampling

Details

model: $w_i = X_i\beta + e$. $e$ $\sim$ $N(0,Sigma)$. note: $w_i, e$ are (p-1) x 1. $y_i = j$, if $w_{ij} > max(0,w_{i,-j})$ j=1,...,p-1. $w_{i,-j}$ means elements of $w_i$ other than the jth. $y_i = p$, if all $w_i < 0$. priors: $beta$ $\sim$ $N(betabar,A^{-1})$ $Sigma$ $\sim$ IW(nu,V) to make up X matrix use createX with DIFF=TRUE. List arguments contain

• p

{number of choices or possible multinomial outcomes} y{n x 1 vector of multinomial outcomes} X{n*(p-1) x k Design Matrix} betabar{k x 1 prior mean (def: 0)} A{k x k prior precision matrix (def: .01I)} nu{ d.f. parm for IWishart prior (def: (p-1) + 3)} V{ pds location parm for IWishart prior (def: nu*I)} beta0{ initial value for beta} sigma0{ initial value for sigma } R{ number of MCMC draws } keep{ thinning parameter - keep every keepth draw (def: 1)}

References

For further discussion, see Bayesian Statistics and Marketing by Allenby, McCulloch, and Rossi, Chapter 4. http://gsbwww.uchicago.edu/fac/peter.rossi/research/bsm.html

See Also

rmvpGibbs

Examples

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

set.seed(66)
p=3
n=500
beta=c(-1,1,1,2)
Sigma=matrix(c(1,.5,.5,1),ncol=2)
k=length(beta)
x1=runif(n*(p-1),min=-1,max=1); x2=runif(n*(p-1),min=-1,max=1)
I2=diag(rep(1,p-1)); xadd=rbind(I2)
for(i in 2:n) { xadd=rbind(xadd,I2)}
X=cbind(xadd,x1,x2)
simout=simmnp(X,p,500,beta,Sigma)

Data=list(p=p,y=simout$y,X=simout$X)
Mcmc=list(R=R,keep=1)

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

cat("Betadraws ",fill=TRUE)
mat=apply(out$betadraw/sqrt(out$sigmadraw[,1]),2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat)
cat("Sigmadraws ",fill=TRUE)
mat=apply(out$sigmadraw/out$sigmadraw[,1],2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(as.vector(Sigma),mat); rownames(mat)[1]="sigma"; print(mat)