# rmvpGibbs: Gibbs Sampler for Multivariate Probit In bayesm: Bayesian Inference for Marketing/Micro-Econometrics

## Description

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

## Usage

 `1` ```rmvpGibbs(Data, Prior, Mcmc) ```

## Arguments

 `Data ` list(y, X, p) `Prior` list(betabar, A, nu, V) `Mcmc ` list(R, keep, nprint, beta0 ,sigma0)

## Details

#### Model and Priors

w_i = X_iβ + e with e ~ N(0,Σ). Note: w_i is p x 1.
y_{ij} = 1 if w_{ij} > 0, else y_i = 0. j = 1, …, p

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

β ~ N(betabar, A^{-1})
Σ ~ IW(nu, V)

To make X matrix use `createX`

#### Argument Details

`Data = list(y, X, p)`

 `X: ` n*p x k Design Matrix `y: ` n*p x 1 vector of 0/1 outcomes `p: ` dimension of multivariate probit

`Prior = list(betabar, A, nu, V)` [optional]

 `betabar: ` k x 1 prior mean (def: 0) `A: ` k x k prior precision matrix (def: 0.01*I) `nu: ` d.f. parameter for Inverted Wishart prior (def: (p-1)+3) `V: ` PDS location parameter for Inverted Wishart prior (def: nu*I)

`Mcmc = list(R, keep, nprint, beta0 ,sigma0)` [only `R` required]

 `R: ` number of MCMC draws `keep: ` MCMC thinning parameter -- keep every `keep`th draw (def: 1) `nprint: ` print the estimated time remaining for every `nprint`'th draw (def: 100, set to 0 for no print) `beta0: ` initial value for beta `sigma0: ` initial value for sigma

## Value

A list containing:

 `betadraw ` R/keep x k matrix of betadraws `sigmadraw ` R/keep x p*p matrix of sigma draws – each row is the vector form of sigma

## Author(s)

Peter Rossi, Anderson School, UCLA, [email protected].

## References

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

`rmnpGibbs`

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45``` ```if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} set.seed(66) 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)) } p = 3 n = 500 beta = c(-2,0,2) Sigma = matrix(c(1, 0.5, 0.5, 0.5, 1, 0.5, 0.5, 0.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) 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) ## plotting examples if(0){plot(betatilde, tvalues=beta/sqrt(diag(Sigma)))} ```

### Example output

```Table of y values
y
0   1
752 748

Starting Gibbs Sampler for MVP
500  obs of  3  binary indicators;  3  indep vars (including intercepts)

Prior Parms:
betabar
[1] 0 0 0
A
[,1] [,2] [,3]
[1,] 0.01 0.00 0.00
[2,] 0.00 0.01 0.00
[3,] 0.00 0.00 0.01
nu
[1] 6
V
[,1] [,2] [,3]
[1,]    6    0    0
[2,]    0    6    0
[3,]    0    0    6

MCMC Parms:
10  reps; keeping every  1 th draw  nprint=  100
initial beta=  0 0 0
initial sigma=
[,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

MCMC Iteration (est time to end - min)
Total Time Elapsed: 0.00