#source("c:/Work07/ShortCourse07/Programs/CauchyPosterior.txt",print.eval=TRUE)#
#Gibbs sampler for Cauchy Example#
#--------------Initial values--------------------------------
p<-function(t)exp(-(t^2)/(2*s2))/((1+(t-x[1])^2)*(1+(t-x[2])^2)*(1+(t-x[3])^2))/0.00406726
x<-c(-6,0,5);nx<-length(x)
nsim<-5000 #run 1000 then 5000
s2<-50		#prior variance
eta<-array(1,c(nx,1))
theta<-array(mean(x),c(nsim,1))
RB<-array(mean(x),c(nsim,1))
#----------------Gibbs Loop--------------------------------
for(i in 2:nsim)
{
	for(j in 1:nx)eta[j]<-rexp(1,rate=.5*(1+(x[j]-theta[i-1])^2))
	RB[i]<-	sum(x*eta)/(sum(eta)+(1/s2))
	theta[i]<-rnorm(1,mean=RB[i],sd=sqrt(1/(sum(eta)+(1/s2))))	
}
#----------------Importance Sampler---------------------------
thetaIS<-rcauchy(nsim, location = mean(x))
thetaIS<-thetaIS*p(thetaIS)/dcauchy(thetaIS, location = mean(x))
#-----------------------Plot----------------------------------
yrange<-c(-2,2)
par(mfrow=c(1,2))
hist(theta,freq=F,xlim=c(-10,10),ylim=c(0,.25),breaks=50,col="grey",ylab="n",main="")
par(new=T)
plot(function(t)p(t),-10,10,ylim=c(0,.25),lwd=2,xlab="n")
plot(cumsum(theta)/c(1:nsim),type="l",ylim=yrange,ylab="")
par(new=T)
plot(cumsum(RB)/c(1:nsim),type="l",ylim=yrange,ylab="",col="red")
par(new=T)
plot(cumsum(thetaIS)/c(1:nsim),type="l",ylim=yrange,ylab="",col="blue")