rmnpGibbs {bayesm} | R Documentation |
rmnpGibbs
implements the McCulloch/Rossi Gibbs Sampler for the multinomial probit model.
rmnpGibbs(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β + e. e ~ N(0,Sigma). note: w_i, e are (p-1) x 1.
y_i = j, if w_{ij} > max(0,w_{i,-j}) j=1,...,p-1. w_{i,-j} means elements of w_i
other than the jth.
y_i = p, if all w_i < 0.
priors:
beta ~ N(betabar,A^{-1})
Sigma ~ IW(nu,V)
to make up X matrix use createX
with DIFF=TRUE
.
List arguments contain
p
y
X
betabar
A
nu
V
beta0
sigma0
R
keep
a list containing:
betadraw |
R/keep x k array of betadraws |
sigmadraw |
R/keep x (p-1)*(p-1) array of sigma draws – each row is in vector form |
beta is not identified. beta/sqrt(sigma_{11}) and Sigma/sigma_{11} are. See Allenby et al or example below for details.
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://faculty.chicagogsb.edu/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(-1,1,1,2) Sigma=matrix(c(1,.5,.5,1),ncol=2) k=length(beta) X1=matrix(runif(n*p,min=0,max=2),ncol=p); X2=matrix(runif(n*p,min=0,max=2),ncol=p) X=createX(p,na=2,nd=NULL,Xa=cbind(X1,X2),Xd=NULL,DIFF=TRUE,base=p) simmnp= function(X,p,n,beta,sigma) { indmax=function(x) {which(max(x)==x)} Xbeta=X%*%beta w=as.vector(crossprod(chol(sigma),matrix(rnorm((p-1)*n),ncol=n)))+ Xbeta w=matrix(w,ncol=(p-1),byrow=TRUE) maxw=apply(w,1,max) y=apply(w,1,indmax) y=ifelse(maxw < 0,p,y) return(list(y=y,X=X,beta=beta,sigma=sigma)) } simout=simmnp(X,p,500,beta,Sigma) Data1=list(p=p,y=simout$y,X=simout$X) Mcmc1=list(R=R,keep=1) out=rmnpGibbs(Data=Data1,Mcmc=Mcmc1) cat(" Summary of Betadraws ",fill=TRUE) betatilde=out$betadraw/sqrt(out$sigmadraw[,1]) attributes(betatilde)$class="bayesm.mat" summary(betatilde,tvalues=beta) cat(" Summary of Sigmadraws ",fill=TRUE) sigmadraw=out$sigmadraw/out$sigmadraw[,1] attributes(sigmadraw)$class="bayesm.var" summary(sigmadraw,tvalues=as.vector(Sigma[upper.tri(Sigma,diag=TRUE)])) if(0){ ## plotting examples plot(betatilde,tvalues=beta) }