#######  Code to fit the model    ########:
#y[i] =  Bern[Phi(x[i,]%*%beta)]
#x[i,] = covariates for sub i


#Draw samples from a truncated normal:
rtnorm<-function(n,mu,sigma,lower,upper){ 
   lp<-pnorm(lower,mu,sigma) 
   up<-pnorm(upper,mu,sigma)  
   qnorm(runif(n,lp,up),mu,sigma) 
}


probit<-function(y,x,sd.beta=100,
            iters=10000,burn=1000,update=10){

  n<-length(y)
  p<-ncol(x)
  low<-ifelse(y==1,0,-Inf)
  high<-ifelse(y==1,Inf,0)

  #Initial values
  z<-y-.5
  beta<-rep(0,p)

  #store samples here
  keep.beta<-matrix(0,iters,p)
  keep.z<-matrix(0,iters,10)

  #Do some matrix manipulations offline
  txx<-t(x)%*%x
  cov.beta<-solve(txx+diag(p)/sd.beta^2)
  P1<-cov.beta%*%t(x)
  P2<-t(chol(cov.beta))

  #Let's go!
  for(i in 1:iters){

    #update the latent probit variables, z:
    mn<-x%*%beta
    z<-rtnorm(n,mn,1,low,high)

    #update beta:
    beta<-P1%*%z+P2%*%rnorm(p)

    keep.beta[i,]<-beta
    keep.z[i,]<-z[1:10]
   
  }


list(beta=keep.beta,
     z=keep.z)}


#Simulated binary data

set.seed(0820)

n<-500
p<-10
true.beta<-seq(0,1,length=p)
x<-matrix(rnorm(n*p),n,p)
x[,1]<-1
p<-pnorm(x%*%true.beta)
y<-rbinom(n,1,p)


#Fit the model
fit<-probit(y,x)

#Plot results
boxplot(fit$beta,main="Posterior of beta")
lines(true.beta)
boxplot(fit$z,main="Posterior of z for 10 observations")