set.seed(323)
n<-15
Y<-rweibull(n,shape=5,scale=10)

(1/10)^5

gpts<-seq(0,1,by=0.01)
al.vec<-gpts*6+2
be.vec<-gpts*10+5

like.mat<-matrix(0,nrow=length(be.vec),ncol=length(al.vec))

for(i in 1:length(be.vec)){
	for(j in 1:length(al.vec)){
		a<-al.vec[j]
		b<-be.vec[i]
		like.mat[j,i]<-sum(dweibull(Y,shape=a,scale=b,log=T))-0.01*a-0.1*(1/b)^a+log(a)-(a+1)*log(b)
	}
}


#col=gray(rev((13:32)/32)),add=T

postscript("Gibbs-Weibull1.ps")
par(pty="s")
plot(al.vec,be.vec,type="n",xlab=expression(gamma),ylab=expression(lambda))
image(al.vec,be.vec,like.mat,zlim=range(-100,-20),col=gray(rev((8:32)/32)),add=T)
dev.off()

ahat<-al.vec[which.max(apply(like.mat,1,max))]
bhat<-be.vec[which.max(apply(like.mat,2,max))]
abline(v=ahat)
abline(h=bhat)

########################################################################

set.seed(23)
nreps<-2000
old.al<-5.5
old.be<-8
sig.al<-sig.be<-1


old.like<-sum(dweibull(Y,shape=old.al,scale=old.be,log=T))
old.prior<--0.01*old.al-0.01*(1/old.be)^old.al+log(old.al)-(old.al+1)*log(old.be)
old.post<-old.like+old.prior

post.samp1<-matrix(0,nrow=nreps,2)
for(irep in 1:nreps){

	new.al<-abs(rnorm(1,old.al,sig.al))
	new.be<-old.be
	new.like<-sum(dweibull(Y,shape=new.al,scale=new.be,log=T))
	new.prior<--0.01*new.al-0.01*(1/new.be)^new.al+log(new.al)-(new.al+1)*log(new.be)
	new.post<-new.like+new.prior
	if(log(runif(1)) < new.post-old.post){
		old.al<-new.al
		old.be<-new.be
		old.like<-new.like
		old.prior<-new.prior
		old.post<-new.post
	}
	post.samp1[irep,1]<-old.al

	new.al<-old.al
	new.be<-abs(rnorm(1,old.be,sig.be))
	new.like<-sum(dweibull(Y,shape=new.al,scale=new.be,log=T))
	new.prior<--0.01*new.al-0.01*(1/new.be)^new.al+log(new.al)-(new.al+1)*log(new.be)
	new.post<-new.like+new.prior
	if(log(runif(1)) < new.post-old.post){
		old.al<-new.al
		old.be<-new.be
		old.like<-new.like
		old.prior<-new.prior
		old.post<-new.post
	}
	post.samp1[irep,2]<-old.be
}

postscript("Gibbs-Weibull1-Post.ps")
par(pty="s")
plot(al.vec,be.vec,type="n",xlab=expression(gamma),ylab=expression(lambda))
image(al.vec,be.vec,like.mat,zlim=range(-100,-20),col=gray(rev((8:32)/32)),add=T)
points(post.samp1,pch=19,cex=0.75)
dev.off()

postscript("Gibbs-Weibull1-Trace.ps")
par(mfrow=c(1,2))
plot(1:nreps,post.samp1[,1],type="l",xlab="Iteration",ylab=expression(gamma[0]))
plot(1:nreps,post.samp1[,2],type="l",xlab="Iteration",ylab=expression(lambda[0]))
dev.off()

postscript("Gibbs-Weibull1-acf.ps")
par(mfrow=c(1,2))
acf(post.samp1[,1],main=expression(gamma[0]))
acf(post.samp1[,2],main=expression(lambda[0]))
dev.off()

qqplot(post.samp0[,1],post.samp1[,1])


yvec<-seq(0,15,by=0.01)
Smat<-matrix(0,nrow=nreps,ncol=length(yvec))
for(i in 1:nreps){
	Smat[i,]<-exp(-(yvec/post.samp1[i,2])^post.samp1[i,1])
}

Post.survivor<-apply(Smat,2,quantile,probs=c(0.025,0.5,0.975))
postscript("Gibbs-Weibull1-Surv.ps")
plot(yvec,Post.survivor[2,],type="n",xlab="y",ylab="Prob")
polygon(c(yvec,rev(yvec)),c(Post.survivor[1,],rev(Post.survivor[3,])),col="gray")
lines(yvec,Post.survivor[2,],col="red",lwd=2)
lines(yvec,Post.survivor[1,],col="red",lwd=2,lty=2)
lines(yvec,Post.survivor[3,],col="red",lwd=2,lty=2)
rug(Y)
lines(c(0,sort(Y)),1-c(0:n)/n,type="s")
dev.off()


########################################################################

gpts<-seq(0,1,by=0.01)
al.vec0<-gpts*4+4
be.vec0<-gpts*1e-5

like.mat0<-matrix(0,nrow=length(be.vec0),ncol=length(al.vec0))

for(i in 1:length(be.vec0)){
	for(j in 1:length(al.vec0)){
		a<-al.vec0[j]
		b<-be.vec0[i]
		like.mat0[j,i]<-sum(log(a*b*Y^(a-1))-b*Y^a)-0.01*(a+b)
	}
}

postscript("Gibbs-Weibull0.ps")
par(pty="s")
plot(al.vec0,be.vec0,type="n",xlab=expression(gamma[0]),ylab=expression(lambda[0]))
image(al.vec0,be.vec0,like.mat0,zlim=range(-100,-20),col=gray(rev((8:32)/32)),add=T)
dev.off()

ahat0<-al.vec0[which.max(apply(like.mat0,1,max))]
bhat0<-be.vec0[which.max(apply(like.mat0,2,max))]
abline(v=ahat0)
abline(h=bhat0)




########################################################################

set.seed(23)
nreps<-2000
old.al<-5.5
old.be<-1e-6
sig.al<-1
sig.be<-1e-6


old.like<-sum(log(old.al*old.be*Y^(old.al-1))-old.be*Y^old.al)
old.prior<--0.01*(old.al+old.be)
old.post0<-old.like+old.prior

post.samp0<-matrix(0,nrow=nreps,2)
for(irep in 1:nreps){

	new.al<-abs(rnorm(1,old.al,sig.al))
	new.be<-old.be
	new.like<-sum(log(new.al*new.be*Y^(new.al-1))-new.be*Y^new.al)
	new.prior<--0.01*(new.al+new.be)
	new.post<-new.like+new.prior
	if(log(runif(1)) < new.post-old.post){
		old.al<-new.al
		old.be<-new.be
		old.like<-new.like
		old.prior<-new.prior
		old.post<-new.post
	}
	post.samp0[irep,1]<-old.al

	new.al<-old.al
	new.be<-abs(rnorm(1,old.be,sig.be))
	new.like<-sum(log(new.al*new.be*Y^(new.al-1))-new.be*Y^new.al)
	new.prior<--0.01*(new.al+new.be)
	new.post<-new.like+new.prior
	if(log(runif(1)) < new.post-old.post){
		old.al<-new.al
		old.be<-new.be
		old.like<-new.like
		old.prior<-new.prior
		old.post<-new.post
	}
	post.samp0[irep,2]<-old.be
}


postscript("Gibbs-Weibull0-Post.ps")
par(pty="s")
plot(al.vec0,be.vec0,type="n",xlab=expression(gamma[0]),ylab=expression(lambda[0]))
image(al.vec0,be.vec0,like.mat0,zlim=range(-100,-20),col=gray(rev((8:32)/32)),add=T)
points(post.samp0,pch=19,cex=0.75)
dev.off()

postscript("Gibbs-Weibull0-Trace.ps")
par(mfrow=c(1,2))
plot(1:nreps,post.samp0[,1],type="l",xlab="Iteration",ylab=expression(gamma[0]))
plot(1:nreps,post.samp0[,2],type="l",xlab="Iteration",ylab=expression(lambda[0]))
dev.off()

postscript("Gibbs-Weibull0-acf.ps")
par(mfrow=c(1,2))
acf(post.samp0[,1],main=expression(gamma[0]))
acf(post.samp0[,2],main=expression(lambda[0]))
dev.off()


################################################################################

library(coda)
effectiveSize(post.samp0[,1])
effectiveSize(post.samp0[,2])
effectiveSize(post.samp1[,1])
effectiveSize(post.samp1[,2])


nreps/(1+2*sum(acf(post.samp1[,1],lag.max=100,plot=F)$acf))

rho.11<-acf(post.samp1[,1],plot=F)$acf[2]

1/(1+