# setting values 
set.seed(1)
s2<-1 
t2<-10 
mu<-5; 
n<-5

# defining the data
x<-c(9.37, 10.18, 9.16, 11.60, 10.33)
# mean of the normal posterior
mu.n<-( mean(x)*n/s2 + mu/t2 )/( n/s2+1/t2) 
# variance of the normal posterior 
t2.n<-1/(n/s2+1/t2)



####metropolis part####
##S = total num of simulations
theta<-0 ; delta<-1 ; S<-500 ; THETA<-NULL
count = 0;
for(s in 1:S)
{
  ## simulating our proposal
  theta.star<-rnorm(1,theta,sqrt(delta))
  ##taking the log of the ratio r
  log.r<-( sum(dnorm(x,theta.star,sqrt(s2),log=TRUE)) +
             dnorm(theta.star,mu,sqrt(t2),log=TRUE) ) -
    ( sum(dnorm(x,theta,sqrt(s2),log=TRUE)) +
        dnorm(theta,mu,sqrt(t2),log=TRUE) )
  if(log(runif(1))<log.r) { 
    theta<-theta.star 
    count = count +1;
  }
  THETA<-c(THETA,theta) ##updating THETA
  ## contains theta's that are accepted and rejected
}
count/S

# two plots: trace of theta and 
# comparing the empirical distribution
# of simulated values to the true posterior
pdf("metropolis_normal-500.pdf",
    family="Times",height=3.5,width=7)
par(mar=c(3,3,1,1),mgp=c(1.75,.75,0))
par(mfrow=c(1,2))
# creating a sequence
skeep<-seq(10,S,by=10)
# making a trace place
plot(skeep,THETA[skeep],type="l",
     xlab="iteration",ylab=expression(theta))
# making a histogram
hist(THETA[-(1:50)],prob=TRUE,main="",
     xlab=expression(theta),ylab="density")
th<-seq(min(THETA),max(THETA),length=100)
lines(th,dnorm(th,mu.n,sqrt(t2.n)) )
dev.off()


S <- 10000
for(s in 1:S)
{
  ## simulating our proposal
  theta.star<-rnorm(1,theta,sqrt(delta))
  ##taking the log of the ratio r
  log.r<-( sum(dnorm(x,theta.star,sqrt(s2),log=TRUE)) +
             dnorm(theta.star,mu,sqrt(t2),log=TRUE) ) -
    ( sum(dnorm(x,theta,sqrt(s2),log=TRUE)) +
        dnorm(theta,mu,sqrt(t2),log=TRUE) )
  if(log(runif(1))<log.r) { theta<-theta.star }
  ##updating THETA
  THETA<-c(THETA,theta)
}
# two plots: trace of theta and 
# comparing the empirical distribution
# of simulated values to the true posterior
pdf("metropolis_normal-10000.pdf",
    family="Times",height=3.5,width=7)
par(mar=c(3,3,1,1),mgp=c(1.75,.75,0))
par(mfrow=c(1,2))
# creating a sequence
skeep<-seq(10,S,by=10)
# making a trace place
plot(skeep,THETA[skeep],type="l",
     xlab="iteration",ylab=expression(theta))
# making a histogram
hist(THETA[-(1:50)],prob=TRUE,main="",
     xlab=expression(theta),ylab="density")
th<-seq(min(THETA),max(THETA),length=100)
lines(th,dnorm(th,mu.n,sqrt(t2.n)) )
dev.off()