bayesm (version 3.0-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,nprint) (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

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

• pnumber of choices or possible multinomial outcomes

• yn x 1 vector of multinomial outcomes

• Xn*(p-1) x k Design Matrix

• betabark x 1 prior mean (def: 0)

• Ak 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)

• nprint print the estimated time remaining for every nprint'th draw (def: 100)

References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch, Chapter 4. http://www.perossi.org/home/bsm-1

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

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

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


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