rmvpGibbs {bayesm} | R Documentation |
rmvpGibbs
implements the Edwards/Allenby Gibbs Sampler for the multivariate probit model.
rmvpGibbs(Data, Prior, Mcmc)
Data |
list(p,y,X) |
Prior |
list(betabar,A,nu,V) (optional) |
Mcmc |
list(beta0,sigma0,R,keep) (R required) |
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
X
y
betabar
A
nu
V
beta0
sigma0
R
keep
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 |
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.
Peter Rossi, Graduate School of Business, University of Chicago, Peter.Rossi@ChicagoGsb.edu.
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
## 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)