rordprobitGibbs {bayesm} | R Documentation |
rordprobitGibbs
implements a Gibbs Sampler for the ordered probit model.
rordprobitGibbs(Data, Prior, Mcmc)
Data |
list(X, y, k) |
Prior |
list(betabar, A, dstarbar, Ad) |
Mcmc |
list(R, keep, s, change, draw) |
Model: z = Xβ + e. e ~ N(0,I).
y=1,..,k. cutoff=c( c [1] ,..c [k+1] ).
y=k, if c [k] <= z < c [k+1] .
Prior: β ~ N(betabar,A^{-1}). dstar ~ N(dstarbar,Ad^{-1}).
List arguments contain
X
y
k
betabar
A
dstarbar
Ad
s
R
keep
betadraw |
R/keep x k matrix of betadraws |
cutdraw |
R/keep x (k-1) matrix of cutdraws |
dstardraw |
R/keep x (k-2) matrix of dstardraws |
accept |
a value of acceptance rate in RW Metropolis |
set c[1]=-100. 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(datsr[k-2])
Be careful in assessing prior parameter, Ad. .1 is too small for many applications.
Peter Rossi, Graduate School of Business, University of Chicago, Peter.Rossi@ChicagoGsb.edu.
Bayesian Statistics and Marketing
by Rossi, Allenby and McCulloch
http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html
## ## rordprobitGibbs example ## if(nchar(Sys.getenv("LONG_TEST")) != 0) {R=2000} else {R=10} ## 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 )) } set.seed(66) 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]) if(0){ ## plotting examples plot(out$cutdraw) }