Hastings Metropolis for N(0,1)

# cible N(0,1)
# Metropolis-Hastings

N <- 10^5
x <- rep(0,N)
x[1] <- -3
pi <- function(x) {exp(-x^2/2)}
sigma2 <- 2

for (t in 2:N)
{
  xtilde <- rnorm(1,x[t-1],sqrt(sigma2))
  # xtilde <- rt(1,3)
  rho <- pi(xtilde)/pi(x[t-1])
  # rho <- pi(xtilde)/pi(x[t-1])*dt(x[t-1],3)/dt(xtilde,3)
  if (runif(1)<=rho) x[t] <- xtilde else x[t] <- x[t-1]
}

plot(1:N,x,type="l")
plot(acf(x,lag=100))

sub <- seq(101,N,by=20)
length(sub)
hist(x[sub],prob=TRUE,nclass=40)
curve(dnorm,-4,4,col="red",add=TRUE)