aa<-c(0,.2,.4,.6,.8,1)
ttt<-scan("http://madison.byu.edu/bayes/ballbearing2.dat")
n<-length(ttt)
ttt<-ttt/1000
nburn<-1000
nsamp<-10000
niter<-nburn+nsamp
lambda<-rep(0,niter)
beta<-rep(0,niter)
lambda[1]<-2
beta[1]<-.01
cs.l<-1.0
cs.b<-0.02
a.l<-2
b.l<-1
a.b<-1
b.b<-1
ar<-rep(0,2)
for(i in 2:niter){
  y<-rnorm(1,lambda[i-1],cs.l)
  lambda[i]<-lambda[i-1]
  if(y>0){
    llo<-sum(dweibull(ttt,lambda[i-1],beta[i-1],log=T))+
      dgamma(lambda[i-1],a.l,scale=b.l,log=T)
    lln<-sum(dweibull(ttt,y,beta[i-1],log=T))+
      dgamma(y,a.l,scale=b.l,log=T)
    if((lln-llo)>log(runif(1,0,1))){ lambda[i]<-y;ar[1]<-ar[1]+1}
  }
  y<-rnorm(1,beta[i-1],cs.b)
  beta[i]<-beta[i-1]
  if(y>0){ 
    llo<-sum(dweibull(ttt,lambda[i],beta[i-1],log=T))+
      dgamma(beta[i-1],a.b,scale=b.b,log=T)
    lln<-sum(dweibull(ttt,lambda[i],y,log=T))+
      dgamma(y,a.b,scale=b.b,log=T)
    if((lln-llo)>log(runif(1,0,1))) {beta[i]<-y;ar[2]<-ar[2]+1}
  }
}
ar<-ar/niter
good<-(nburn+1):niter
ngood<-length(good)
bincnts<-rep(0,length(aa)-1)
pk<-diff(aa)

BX2<-rep(0,length(ngood))
for(i in 1:ngood){
  quants<-qweibull(aa,lambda[i+nburn],beta[i+nburn])
  
  for(j in 1:(length(aa)-1)){
    bincnts[j]<-sum(ttt>quants[j]&ttt<quants[j+1])
  }
  BX2[i]<-sum(((bincnts-n*pk)/sqrt(n*pk))^2)
}