#data augmented Gibbs sampler probit
library(bayesm)
n<-1000
x1<-runif(n)
x2<-runif(n)
V<-rnorm(n)
theta<-c(-1,1,1)
X<-cbind(rep(1,n),x1,x2)
UD<-X%*%theta-V
D<-sign(UD)*.5+.5
sum(D)

#mle
beta.probit<-glm(D~0+X,family=quasi(var="mu(1-mu)", link="probit"), start=rep(0,3))
summary(beta.probit)
1-beta.probit$deviance/beta.probit$null

#Gibbs probit
R<-5000
Data1=list(X=X,y=D)
Mcmc1=list(R=R,keep=1)

out=rbprobitGibbs(Data=Data1,Mcmc=Mcmc1)
beta.sample<-out$betadraw[1001:R,]
par(mfrow=c(3,1),mar=c(1,2,1,1))
for(j in 1:3)plot(beta.sample[,j])	#trace plots
apply(beta.sample,2,mean)
apply(beta.sample,2,median)
apply(beta.sample,2,sd)
apply(beta.sample,2,quantile,c(.01,.025,.05,.1,.25,.75,.9,.95,.975,.99))
apply(beta.sample,2,max)
apply(beta.sample,2,min)