#Markov chain Monte Carlo simulation for Normal distribution

#Random variable generator
tmean <- 50
tsig2 <- 5
y <- rnorm(100,tmean,tsig2)
n <- length(y)
ybar <- mean(y)
sigma2 <- var(y)

#Log posterior distribution
logposterior <- function(x1, x2){
-1/2*(n*x2+((n-1)*sigma2+n*(ybar-x1)^2)/exp(x2))
}

#Simulation conditions
x1 <- 1
x2 <- log(1)
lp <- logposterior(x1,x2)
N=100000
x1trace <- x2trace <- numeric(N)

#MCMC simulation
for(i in 1:N){
　y1 <- rnorm(1, x1, 0.5)
　y2 <- rnorm(1, x2, 0.5)
　lq <- logposterior(y1,y2)
　if(lp-lq < rexp(1)){
　　x1 <- y1
　　x2 <- y2
　　lp <- lq
　}
　x1trace[i] <- x1
　x2trace[i] <- x2
}

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

#X-Y plot
plot(x1trace, x2trace, type="l", col="gray45")

#Histgram
hist((x1trace-ybar)/sqrt(sigma2/n),freq=FALSE, main="",col="cyan", xlim = c(-5,5),breaks=seq(-1000,1000,0.2))
curve(dt(x,n-1),add=TRUE,lwd=2,col="red") #t-distribution with d.f.=n-1

#Bayesian estimator
cat("Original mean ")
print(tmean)
cat("Original variance ")
print(tsig2*tsig2)
cat("Bayesian mean ")
print(median(x1trace))
cat("Bayesian variance ")
print(median(sigma2))

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