Résumé de section

    • # 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)