# simulations suivant une loi exponentielle de
# parametre 2

u <-  runif(1000000)
x <- -log(u)/2
hist(x,nclass=100,prob=TRUE)
curve(dexp(x,2),0,5,col="red",lwd=2,add=TRUE)

hist(pexp(x,2))
x <-qexp(u,2)

# simulations suivant une loi normale centree
# reduite

N <- 10^6
x <- rep(0,N)
for (i in 1:N)
{
u1 <- runif(1)
u2 <- runif(1)
x[i] <- sqrt(-2*log(u1))*cos(u2*2*pi)
}
hist(x,nclass=100,prob=TRUE)
curve(dnorm,-3,3,col="red",lwd=2,add=TRUE)

# simulations suivant une loi gausienne 
# multivariee

Sigma <- matrix(c(1,0.9,0.9,1),2,2) 
Sigma

mu <- c(3,5)
Sigma12 <- t(chol(Sigma))
X <- matrix(x[1:2000],1000,2)

Y <- t(mu+Sigma12%*%t(X))
plot(Y)

# simulations suivant une loi gamma(2,1)
# algorithme acceptation-rejet
# loi instrumentale optimale parmi les lois
# exponentielles : Exp(1/2)

test <- TRUE ; k <- 4*exp(-1) ; N <- 100000
x <- rep(0,N) ; compteur <- 0
for (i in 1:N)
{
  while (test)
  {
    y <- rexp(1,1/2)
    paccept <- dgamma(y,2,1)/(dexp(y,1/2)*k)
    if (runif(1) <= paccept) test <- FALSE
    compteur <- compteur+1
  }
  x[i] <- y
  test <- TRUE
}
hist(x,prob=TRUE,nclass=50)
curve(dgamma(x,2,1),col="red",lwd=4,add=TRUE)
curve(dexp(x,1/2),col="green",lwd=4,add=TRUE)
N/compteur

# loi instrumentale Exp(1/10)

test <- TRUE ; k <- 90 ; N <- 1000
x <- rep(0,N) ; compteur <- 0
for (i in 1:N)
{
  while (test)
  {
    y <- rexp(1,1/10)
    paccept <- dgamma(y,2,1)/(dexp(y,1/10)*k)
    if (runif(1) <= paccept) test <- FALSE
    compteur <- compteur+1
  }
  x[i] <- y
  test <- TRUE
}
hist(x,prob=TRUE,nclass=50)
curve(dgamma(x,2,1),0.1,8,col="red",lwd=4,add=TRUE)
curve(dexp(x,1/10),0.1,8,col="green",lwd=4,add=TRUE)
N/compteur

# cible (à une constante près ) 
# x*exp(-x)(1+cos^2(x))
# loi instrumentale Exp(1/2)

test <- TRUE ; ktilde <- 8 ; N <- 100000
x <- rep(0,N) ; compteur <- 0
for (i in 1:N)
{
  while (test)
  {
    y <- rexp(1,1/2)
    paccept <- dgamma(y,2,1)*(1+cos(y)^2)/
      (dexp(y,1/2)*ktilde)
    if (runif(1) <= paccept) test <- FALSE
    compteur <- compteur+1
  }
  x[i] <- y
  test <- TRUE
}

ptilde <- function(x) dgamma(x,2,1)*(1+cos(x)^2)
c <- integrate(ptilde,0,20)$val

hist(x,prob=TRUE,nclass=50)
curve(ptilde(x)/c,col="red",lwd=4,add=TRUE)
curve(dexp(x,1/2),col="green",lwd=4,add=TRUE)
N/compteur*Mtilde