bayesm (version 1.1-2)

rmvpGibbs: Gibbs Sampler for Multivariate Probit

Description

rmvpGibbs implements the Edwards/Allenby Gibbs Sampler for the multivariate probit model.

Usage

rmvpGibbs(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*p array of sigma draws -- each row is in vector form

concept

• bayes
• multivariate probit
• MCMC
• Gibbs Sampling

Details

model: $w_i = X_i beta + e$. $e$ $\sim$ N(0,Sigma). note: $w_i$ is p x 1. $y_{ij} = 1$, if $w_{ij} > 0$, else $y_i=0$. j=1,...,p. priors: $beta$ $\sim$ $N(betabar,A^{-1})$ $Sigma$ $\sim$ IW(nu,V) to make up X matrix use createX List arguments contain

• p

{dimension of multivariate probit} X{n*p x k Design Matrix} y{n*p x 1 vector of 0,1 outcomes} 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

rmnpGibbs

Examples

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

set.seed(66)
p=3
n=500
beta=c(-2,0,2)
Sigma=matrix(c(1,.5,.5,.5,1,.5,.5,.5,1),ncol=3)
k=length(beta)
I2=diag(rep(1,p)); xadd=rbind(I2)
for(i in 2:n) { xadd=rbind(xadd,I2)}; X=xadd
simout=simmvp(X,p,500,beta,Sigma)

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

ind=seq(from=0,by=p,length=k)
inda=1:3
ind=ind+inda
cat("Betadraws ",fill=TRUE)
mat=apply(out$betadraw/sqrt(out$sigmadraw[,ind]),2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(beta,mat); rownames(mat)[1]="beta"; print(mat)
rdraw=matrix(double((R)*p*p),ncol=p*p)
rdraw=t(apply(out$sigmadraw,1,nmat))
cat("Draws of Correlation Matrix ",fill=TRUE)
mat=apply(rdraw,2,quantile,probs=c(.01,.05,.5,.95,.99))
mat=rbind(as.vector(Sigma),mat); rownames(mat)[1]="Sigma"; print(mat)