Gibbs Sampler (Albert and Chib) for Binary Probit
rbprobitGibbs implements the Albert and Chib Gibbs Sampler for the binary probit model.
rbprobitGibbs(Data, Prior, Mcmc)
Data |
list(y, X) |
Prior |
list(betabar, A) |
Mcmc |
list(R, keep, nprint) |
z = Xβ + e with e ~ N(0, I)
y = 1 if z > 0
β ~ N(betabar, A^{-1})
Data = list(y, X)
y: |
n x 1 vector of 0/1 outcomes |
X: |
n x k design matrix |
Prior = list(betabar, A) [optional]
betabar: |
k x 1 prior mean (def: 0) |
A: |
k x k prior precision matrix (def: 0.01*I) |
Mcmc = list(R, keep, nprint) [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)
|
A list containing:
betadraw |
R/keep x k matrix of betadraws |
Peter Rossi, Anderson School, UCLA, perossichi@gmail.com.
For further discussion, see Chapter 3, 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)
## function to simulate from binary probit including x variable
simbprobit = function(X, beta) {
y = ifelse((X%*%beta + rnorm(nrow(X)))<0, 0, 1)
list(X=X, y=y, beta=beta)
}
nobs = 200
X = cbind(rep(1,nobs), runif(nobs), runif(nobs))
beta = c(0,1,-1)
nvar = ncol(X)
simout = simbprobit(X, beta)
Data1 = list(X=simout$X, y=simout$y)
Mcmc1 = list(R=R, keep=1)
out = rbprobitGibbs(Data=Data1, Mcmc=Mcmc1)
summary(out$betadraw, tvalues=beta)
## plotting example
if(0){plot(out$betadraw, tvalues=beta)}Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.