Gibbs Sampler for Multinomial Probit
rmnpGibbs
implements the McCulloch/Rossi Gibbs Sampler for the multinomial probit model.
rmnpGibbs(Data, Prior, Mcmc)
Data |
list(y, X, p) |
Prior |
list(betabar, A, nu, V) |
Mcmc |
list(R, keep, nprint, beta0, sigma0) |
w_i = X_iβ + e with e ~ N(0, Σ).
Note: w_i and e are (p-1) x 1.
y_i = j if w_{ij} > max(0,w_{i,-j}) for j=1, …, p-1 and
where w_{i,-j} means elements of w_i other than the jth.
y_i = p, if all w_i < 0
β is not identified. However, β/sqrt(σ_{11}) and Σ/σ_{11} are identified. See reference or example below for details.
β ~ N(betabar,A^{-1})
Σ ~ IW(nu,V)
Data = list(y, X, p)
y: |
n x 1 vector of multinomial outcomes (1, ..., p) |
X: |
n*(p-1) x k design matrix. To make X matrix use createX with DIFF=TRUE |
p: |
number of alternatives |
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 (def: 0) |
sigma0: |
initial value for sigma (def: I) |
A list containing:
betadraw |
R/keep x k matrix of betadraws |
sigmadraw |
R/keep x (p-1)*(p-1) matrix of sigma draws – each row is the vector form of sigma |
Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.
For further discussion, see Chapter 4, Bayesian Statistics and Marketing by Rossi, Allenby, and McCulloch.
http://www.perossi.org/home/bsm-1
if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} set.seed(66) 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)) } p = 3 n = 500 beta = c(-1,1,1,2) Sigma = matrix(c(1, 0.5, 0.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) 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)])) ## plotting examples if(0){plot(betatilde,tvalues=beta)}
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.