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

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

• pdimension of multivariate probit

• Xn*p x k Design Matrix

• yn*p x 1 vector of 0,1 outcomes

• 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

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

simmvp= function(X,p,n,beta,sigma) {
  w=as.vector(crossprod(chol(sigma),matrix(rnorm(p*n),ncol=n)))+ X%*%beta
  y=ifelse(w<0,0,1)
  return(list(y=y,X=X,beta=beta,sigma=sigma))
}

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

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

ind=seq(from=0,by=p,length=k)
inda=1:3
ind=ind+inda
cat(" Betadraws ",fill=TRUE)
betatilde=out$betadraw/sqrt(out$sigmadraw[,ind])
attributes(betatilde)$class="bayesm.mat"
summary(betatilde,tvalues=beta/sqrt(diag(Sigma)))

rdraw=matrix(double((R)*p*p),ncol=p*p)
rdraw=t(apply(out$sigmadraw,1,nmat))
attributes(rdraw)$class="bayesm.var"
tvalue=nmat(as.vector(Sigma))
dim(tvalue)=c(p,p)
tvalue=as.vector(tvalue[upper.tri(tvalue,diag=TRUE)])
cat(" Draws of Correlation Matrix ",fill=TRUE)
summary(rdraw,tvalues=tvalue)

if(0){
plot(betatilde,tvalues=beta/sqrt(diag(Sigma)))
}