bayesm (version 3.0-2)

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, nprint 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

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] .

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

List arguments contain

X

n x nvar Design Matrix

y

n x 1 vector of observations, (1,...,k)

k

the largest possible value of y

betabar

nvar x 1 prior mean (def: 0)

A

nvar x nvar prior precision matrix (def: .01I)

dstarbar

ndstar x 1 prior mean, ndstar=k-2 (def: 0)

Ad

ndstar x ndstar prior precision matrix (def:I)

s

scaling parm for RW Metropolis (def: 2.93/sqrt(nvar))

R

number of MCMC draws

keep

thinning parameter - keep every keepth draw (def: 1)

nprint

print the estimated time remaining for every nprint'th draw (def: 100)

References

Bayesian Statistics and Marketing by Rossi, Allenby and McCulloch http://www.perossi.org/home/bsm-1

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