# cible PI N(0,1)
# Q N(.,sigma2)
# 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))
  rho <- pi(xtilde)/pi(x[t-1])
  if (runif(1)<=rho) x[t] <- xtilde else x[t] <- x[t-1]
}

sub <- 100
plot(1:sub,x[1:sub],type="l")
acf(x,lag=30)

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

## TP gibbs sampling 

n <- 10
alpha <- 1
beta <- 2
gibbs <- function(niter)
{
  x <- matrix(0,niter+1,2)
  x[1,1] <- 5
  x[1,2] <- rbeta(1,x[1,1]+alpha,n-x[1,1]+beta)
  for (t in 2:(niter+1))
  {
    x[t,1] <- rbinom(1,n,x[t-1,2])
    x[t,2] <- rbeta(1,x[t,1]+alpha,n-x[t,1]+beta)
  }
  x
}

simu <- gibbs(10^4)

plot(simu[,1],type="l")
acf(simu[,1])

sub <- seq(1001,10^4,by=20)
length(sub)
barplot(table(as.factor(simu[,1])))

plot(simu[,2],type="l")
acf(simu[,2])
sub <- seq(1001,10^4,by=10)
hist(simu[sub,2],prob=TRUE)
curve(dbeta(x,alpha,beta),0,1,add=TRUE,col="red")

## TP Metropolis-Hastings

alpha <- 4
n <- 5
target <- function(x){x^(-n/2)*exp(-alpha/(2*x))}

hm1 <- function(niter)
{
  x <- rep(0,niter+1)
  x[1] <- rchisq(1,2)
  for (t in 2:(niter+1))
  {
    xtilde <- rchisq(1,2)
    rho <- (target(xtilde)*dchisq(x[t-1],2))/(target(x[t-1])*dchisq(xtilde,2))
    if (runif(1)<=rho) x[t]=xtilde else x[t]=x[t-1]
  }
  x
}

simu <- hm1(10^4)
plot(1:(10^4+1),simu,type="l")
acf(simu,lag=100)
sub <- seq(50,10^4,by=10)
hist(simu[sub],prob=TRUE,nclass=50)
curve(target(x)/integrate(target,0,40)$val,0,20,lwd=2,col="red",add=TRUE)
curve(dchisq(x,2),0,20,lwd=2,col="green",add=TRUE)

test <- acf(simu,lag=1000)
sum(test$acf)


hm2 <- function(niter)
{
  x <- rep(0,niter+1)
  x[1] <- 200
  for (t in 2:(niter+1))
  {
    xtilde <- rnorm(1,x[t-1],sqrt(5))
    if (xtilde>=0) rho <- target(xtilde)/target(x[t-1]) else rho <- 0
    if (runif(1)<=rho) x[t] <- xtilde else x[t] <- x[t-1]
  }
  x
}

simu <- hm2(10^4)
plot(1:10^4,simu[1:10^4],type="l")
acf(simu,lag=100)
sub <- seq(50,10^4,by=10)
hist(simu[sub],prob=TRUE,nclass=50)
curve(target(x)/integrate(target,0,40)$val,0,20,lwd=2,col="red",add=TRUE)

test <- acf(simu,lag=1000)
sum(test$acf)