rmvpGibbs {bayesm} | R Documentation

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

### Details

model:
w_i = X_i beta + e. e ~ 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 ~ N(betabar,A^{-1})
Sigma ~ 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)

### 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

### Note

beta and Sigma are not identifed. Correlation matrix and the betas divided by the appropriate standard deviation are. See Allenby et al for details or example below.

### Author(s)

Peter Rossi, Graduate School of Business, University of Chicago, Peter.Rossi@ChicagoGsb.edu.

### References

For further discussion, see Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch, 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)

[Package bayesm version 2.0-8 Index]