n <- 20
mu <- 1.5
sigma2 <- 2
x <- rnorm(n,mu,sqrt(sigma2))
xbar <- mean(x)

## EXEMPLE 1

# Gibbs sampler pour
# mu uniforme entre -3 et 3
# sigma2 Inverse Gamma (1,1)

Nsimu <- 10000
mcmc.sigma2 <- rep(0,Nsimu)
mcmc.mu <- rep(0,Nsimu)
mcmc.sigma2[1] <- 4
# mcmc.sigma2[1] <- 1/rgamma(1,1,1)

for (i in 2:Nsimu)
{
  bool <- TRUE
  while (bool)
  {
    y <- rnorm(1,xbar,sqrt(mcmc.sigma2[i-1]/n))
    if (y<=3 & y>=-3) { mcmc.mu[i] <- y ; bool <- FALSE}
  }
  mcmc.sigma2[i] <- 1/rgamma(1,n/2+1,sum((x-mcmc.mu[i])^2)/2+1)
}

plot(mcmc.mu,type="l")
acf(mcmc.mu)
mean(mcmc.mu)

plot(mcmc.sigma2,type="l")
acf(mcmc.sigma2)
mean(mcmc.sigma2)

## EXEMPLE 2

# Metropolis-Within-Gibbs sampler pour
# mu uniforme entre -3 et 3
# sigma2 unifirme entre 1 et 5

Nsimu <- 10000
mcmc.sigma2 <- rep(0,Nsimu)
mcmc.mu <- rep(0,Nsimu)
mcmc.sigma2[1] <- 4

target <- function(sigma2,mu) exp(-1/(2*sigma2)*sum((x-mu)^2))*sigma2^(-n/2)

for (i in 2:Nsimu)
{
  bool <- TRUE
  while (bool)
  {
    y <- rnorm(1,xbar,sqrt(mcmc.sigma2[i-1]/n))
    if (y<=3 & y>=-3) { mcmc.mu[i] <- y ; bool <- FALSE}
  }
  sigma2tilde <- runif(1,1,5)
  rho <- target(sigma2tilde,mcmc.mu[i])/target(mcmc.sigma2[i-1],mcmc.mu[i])
  if (runif(1)<=rho) mcmc.sigma2[i] <- sigma2tilde else mcmc.sigma2[i] <- mcmc.sigma2[i-1]
}

plot(mcmc.mu,type="l")
acf(mcmc.mu)
mean(mcmc.mu)

plot(mcmc.sigma2,type="l")
acf(mcmc.sigma2)
mean(mcmc.sigma2)