Gibbs Sampler for Ordered Probit
rordprobitGibbs implements a Gibbs Sampler for the ordered probit model with a RW Metropolis step for the cut-offs.
rordprobitGibbs(Data, Prior, Mcmc)
Data |
list(y, X, k) |
Prior |
list(betabar, A, dstarbar, Ad) |
Mcmc |
list(R, keep, nprint, s) |
z = Xβ + e with e ~ N(0, I)
y = k if c[k] ≤ z < c[k+1] with k = 1,…,K
cutoffs = {c[1], …, c[k+1]}
β ~ N(betabar, A^{-1})
dstar ~ N(dstarbar, Ad^{-1})
Be careful in assessing prior parameter Ad: 0.1 is too small for many applications.
Data = list(y, X, k)
y: |
n x 1 vector of observations, (1, …, k) |
X: |
n x p Design Matrix |
k: |
the largest possible value of y |
Prior = list(betabar, A, dstarbar, Ad) [optional]
betabar: |
p x 1 prior mean (def: 0) |
A: |
p x p prior precision matrix (def: 0.01*I) |
dstarbar: |
ndstar x 1 prior mean, where ndstar=k-2 (def: 0) |
Ad: |
ndstar x ndstar prior precision matrix (def: I) |
Mcmc = list(R, keep, nprint, s) [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) |
s: |
scaling parameter for RW Metropolis (def: 2.93/sqrt(p))
|
A list containing:
betadraw |
R/keep x p matrix of betadraws |
cutdraw |
R/keep x (k-1) matrix of cutdraws |
dstardraw |
R/keep x (k-2) matrix of dstardraws |
accept |
acceptance rate of Metropolis draws for cut-offs |
set c[1] = -100 and c[k+1] = 100. c[2] is set to 0 for identification.
The relationship between cut-offs and dstar is:
c[3] = exp(dstar[1]),
c[4] = c[3] + exp(dstar[2]), ...,
c[k] = c[k-1] + exp(dstar[k-2])
Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.
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)
## simulate data for ordered probit model
simordprobit=function(X, betas, cutoff){
z = X%*%betas + rnorm(nobs)
y = cut(z, br = cutoff, right=TRUE, include.lowest = TRUE, labels = FALSE)
return(list(y = y, X = X, k=(length(cutoff)-1), betas= betas, cutoff=cutoff ))
}
nobs = 300
X = cbind(rep(1,nobs),runif(nobs, min=0, max=5),runif(nobs,min=0, max=5))
k = 5
betas = c(0.5, 1, -0.5)
cutoff = c(-100, 0, 1.0, 1.8, 3.2, 100)
simout = simordprobit(X, betas, cutoff)
Data=list(X=simout$X, y=simout$y, k=k)
## set Mcmc for ordered probit model
Mcmc = list(R=R)
out = rordprobitGibbs(Data=Data, Mcmc=Mcmc)
cat(" ", fill=TRUE)
cat("acceptance rate= ", accept=out$accept, fill=TRUE)
## outputs of betadraw and cut-off draws
cat(" Summary of betadraws", fill=TRUE)
summary(out$betadraw, tvalues=betas)
cat(" Summary of cut-off draws", fill=TRUE)
summary(out$cutdraw, tvalues=cutoff[2:k])
## plotting examples
if(0){plot(out$cutdraw)}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.