#1D Markov chain Monte Carlo simulation for Cauchy distribution

#Simulation conditions
N <- 100000
a <- numeric(N)
x <- -50 #Initial value of x
p <- 1/(1+x^2) #Kernel of standard Cauchy distribution
accept <- 0
b <- rep(NA,N)

#MCMC simulation
for(i in 1:N){
　y <- rnorm(1,x,1)
　q <- 1/(1+y^2)
　if(runif(1) < q/p){
　　x <- y
　　p <- q
　　accept <- accept+1
　}else b[i] <- y
　a[i] <- x
}

#Plotting column
par(mfrow=c(1,2))

#X-Y plot
plot(0:N, c(0,a), pch=16, type="o", ylim=range(c(a,b),na.rm=TRUE))
points(1:N,b,pch=1)
segments(0:(N-1),a,1:N,b)

#Histgram
hist(a,freq=FALSE,breaks=seq(-1000,1000,0.4),xlim=c(-10,10),main="",col="cyan")
curve(dcauchy(x),add=TRUE,lwd=2,col="red")

#Plotting column
par(mfrow=c(1,1))

#Accept ratio
ar<-accept/N
cat("Accept ratio")
print(ar)