#---------------------------------------
# read the senility symptoms/wais data
#---------------------------------------
source('chap02/ex3/chap02_ex3_wais_data.r')
wais <- ex2_3.wais
# -------------------------------------------- 
y<-wais$senility; x<-wais$wais; n<-length(y)
positive<- y==1 
Iterations<-55000
mu.beta<-c(0,0); s.beta<-c(100,100) 
beta <- matrix(nrow=Iterations, ncol=2) 
acc.prob <- 0 
current.beta<-c(0,0); u<-numeric(n)  
for (t in 1:Iterations){  
  eta<-current.beta[1]+current.beta[2]*x
	U<-exp(y*eta)/(1+exp(eta))
	u<-runif( n, rep(0,n), U) 
	logitu<-log( u/(1-u) )
	logitu1<-  logitu[positive]
	logitu2<- -logitu[!positive]
	
	l0<- max( logitu1 - current.beta[2]*x[positive] )
	u0<- min( logitu2 - current.beta[2]*x[!positive] )
	unif.random<-runif(1,0,1)
	fa<- pnorm(l0, mu.beta[1], s.beta[1])
	fb<- pnorm(u0, mu.beta[1], s.beta[1])
	current.beta[1] <- qnorm( fa + unif.random*(fb-fa), mu.beta[1], s.beta[1])
	
	l1<- max( (logitu1 - current.beta[1])/x[positive] )
	u1<- min( (logitu2 - current.beta[1])/x[!positive] )
	unif.random<-runif(1,0,1)
	fa<- pnorm(l1, mu.beta[2], s.beta[2])
	fb<- pnorm(u1, mu.beta[2], s.beta[2])
	current.beta[2] <- qnorm( fa + unif.random*(fb-fa), mu.beta[2], s.beta[2])
	beta[t,]<-current.beta  
}

mcmc.output<-beta