rmvpGibbs: Gibbs Sampler for Multivariate Probit

In bayesm: Bayesian Inference for Marketing/Micro-Econometrics

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(y, X, p)
Prior	list(betabar, A, nu, V)
Mcmc	list(R, keep, nprint, beta0 ,sigma0)

Details

Model and Priors

w_i = X_i\beta + e with 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, \ldots, 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.

\beta \sim N(betabar, A^{-1})
\Sigma \sim 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 keepth 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, [e-mail].

References

For further discussion, see Chapter 4, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.

See Also

rmnpGibbs

Examples

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

bayesm documentation built on Sept. 24, 2023, 1:07 a.m.