# Exemple x1,...,xn iid suivant N(theta,1)
# theta = 2
# n = 20

# loi a priori uniforme entre -3 et 3
# generation des donnees

n <- 20
x <- rnorm(20,2,1)

# cible loi a posteriori

target <- function(theta)
{
  num <- length(theta)
  res <- rep(0,num)
  for (i in 1:num)
    {
    if ((theta[i]<3) & (theta[i]>-3)) res[i] <- exp(-0.5*sum((x-theta[i])^2)) 
    else res[i] <- 0
  }
  res
}

curve(target,-4,4)

cst <- integrate(target,-3,3)$val
curve(target(x)/cst,-3,3)

# z <- seq(-4,4,length=100)
# dens <- rep(0,100)
# for (i in 1:100) dens[i] <- target(z[i])
# plot(z,dens,type="l")

T <- 1000
theta <- rep(0,T)

theta[1] <- 0 # x^(0)
for (t in 2:T)
{
thetatilde <- rnorm(1,theta[t-1],0.5)  # xtilde et Q
if ((thetatilde>3) | (thetatilde<(-3))) theta[t] <- theta[t-1]
else 
  {
  rho <- target(thetatilde)/target(theta[t-1])
  u <- runif(1)
  if (u<=rho) theta[t] <- thetatilde else theta[t] <- theta[t-1]
  }
}

plot(theta,type="l")
acf(theta)
sub <- seq(101,T,by=10)
hist(theta[sub],prob=TRUE)
curve(target(x)/cst,add=TRUE,col="red")

mean(theta[sub])
h <- function(x) target(x)*x
integrate(h,-3,3)$val/cst