prior_mean <- 0
prior_sd <- 1

#####################
data <- c(-1,0,1)
#####################

mean_of_data <- mean(data)
assumed_sd_of_data <- sqrt(var(data))
se_data <- assumed_sd_of_data/sqrt(3)

step <- 2 


iter <-  10000
X<- Y <- x_prior <- y_prior<- x_likelihood <- y_likelihood <- R <-rep(0,iter)
X[1] <- round(runif(1,1,5),0)

for (i in 1:iter) {
  if (i > 1){
  X[i] <- ifelse(runif(1) < R[i-1], Y[i-1], X[i-1]) }
  #v <- (runif(1)-0.5)*step
  v <- runif(1,-1,1)
  Y[i] <- X[i]+v
  x_prior[i] <- dnorm(X[i],prior_mean,prior_sd)
  x_likelihood[i] <- dnorm(data[1],X[i],assumed_sd_of_data)*
                     dnorm(data[2],X[i],assumed_sd_of_data)*
                     dnorm(data[3],X[i],assumed_sd_of_data)
  y_prior[i] <- dnorm(Y[i],prior_mean,prior_sd)
  y_likelihood[i] <- dnorm(data[1],Y[i],assumed_sd_of_data)*
                     dnorm(data[2],Y[i],assumed_sd_of_data)*
                     dnorm(data[3],Y[i],assumed_sd_of_data)
  R[i] <- exp(log(y_prior[i]) + log(y_likelihood[i])-
              log(x_prior[i]) - log(x_likelihood[i]))
  #print(c(X[i],Y[i],x_prior[i],x_likelihood[i],y_prior[i],y_likelihood[i],R[i]))
  
}

par(mfrow=c(1,2))
plot(1:iter,X,ylim=c(-4,4),type="l")
abline(h=0,lty=2)
mean(X[1500:3000])
sd(X[1500:3000])

plot(density(X),xlim=c(-3,5))
lines(density(rnorm(20000,prior_mean,prior_sd)),col="red")
abline(v=mean(X),lty=2)