bayesm (version 2.2-1)

rordprobitGibbs: Gibbs Sampler for Ordered Probit

Description

rordprobitGibbs implements a Gibbs Sampler for the ordered probit model.

Usage

rordprobitGibbs(Data, Prior, Mcmc)

Arguments

Data

list(X, y, k)

Prior

list(betabar, A, dstarbar, Ad)

Mcmc

list(R, keep, s, change, draw)

Value

• 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

concept

• bayes
• MCMC
• probit
• Gibbs Sampling

Details

Model: $z = X\beta + e$. $e$ $\sim$ $N(0,I)$. y=1,..,k. cutoff=c( c [1] ,..c [k+1] ). y=k, if c [k] <= z="" <="" c="" [k+1]="" .="" p="">

Prior: $\beta$ $\sim$ $N(betabar,A^{-1})$. $dstar$ $\sim$ $N(dstarbar,Ad^{-1})$.

List arguments contain [object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object],[object Object]

References

Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch http://faculty.chicagogsb.edu/peter.rossi/research/bsm.html

See Also

rbprobitGibbs

Examples

##
## 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)
}