# Simulations following an exponential distribution with parameter 2

u <- runif(1000000)  # Generate 1,000,000 uniform random numbers between 0 and 1
x <- -log(u) / 2  # Transform using the inverse CDF of the exponential distribution
hist(x, nclass = 100, prob = TRUE)  # Plot histogram of the exponential samples
curve(dexp(x, 2), 0, 5, col = "red", lwd = 2, add = TRUE)  # Overlay the theoretical density curve of Exp(2)

hist(pexp(x, 2))  # Histogram of the CDF values of the exponential samples
x <- qexp(u, 2)  # Quantile transformation using the inverse CDF of Exp(2)

# Simulations following a standard normal distribution

N <- 10^6  # Number of samples to generate
x <- rep(0, N)  # Initialize a vector to hold the normal samples
for (i in 1:N) {
  u1 <- runif(1)  # Generate two independent uniform random numbers
  u2 <- runif(1)
  x[i] <- sqrt(-2 * log(u1)) * cos(u2 * 2 * pi)  # Apply the Box-Muller transform to generate normal samples
}
hist(x, nclass = 100, prob = TRUE)  # Plot histogram of the normal samples
curve(dnorm, -3, 3, col = "red", lwd = 2, add = TRUE)  # Overlay the theoretical normal density

# Simulations following a multivariate Gaussian distribution

Sigma <- matrix(c(1, 0.9, 0.9, 1), 2, 2)  # Covariance matrix with high correlation (0.9) between variables
Sigma  # Display the covariance matrix

mu <- c(3, 5)  # Mean vector of the distribution
Sigma12 <- t(chol(Sigma))  # Cholesky decomposition to generate correlated normal samples
X <- matrix(x[1:2000], 1000, 2)  # Reshape normal samples into a matrix

Y <- t(mu + Sigma12 %*% t(X))  # Apply the transformation to obtain multivariate Gaussian samples
plot(Y)  # Plot the samples

# Simulations following a Gamma(2,1) distribution using the acceptance-rejection algorithm with the optimal instrumental distribution among exponential distributions: 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)  # Sample from the instrumental Exp(1/2) distribution
    paccept <- dgamma(y, 2, 1) / (dexp(y, 1/2) * k)  # Calculate the acceptance probability
    if (runif(1) <= paccept) test <- FALSE  # Accept or reject the sample
    compteur <- compteur + 1
  }
  x[i] <- y
  test <- TRUE
}
hist(x, prob = TRUE, nclass = 50)  # Plot the histogram of the accepted Gamma samples
curve(dgamma(x, 2, 1), col = "red", lwd = 4, add = TRUE)  # Overlay the theoretical Gamma density
curve(dexp(x, 1/2), col = "green", lwd = 4, add = TRUE)  # Overlay the instrumental distribution
N / compteur  # Calculate the efficiency of the acceptance-rejection method

# Instrumental law 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)  # Sample from the instrumental Exp(1/10) distribution
    paccept <- dgamma(y, 2, 1) / (dexp(y, 1/10) * k)  # Acceptance probability
    if (runif(1) <= paccept) test <- FALSE  # Accept or reject the sample
    compteur <- compteur + 1
  }
  x[i] <- y
  test <- TRUE
}
hist(x, prob = TRUE, nclass = 50)  # Plot the histogram of accepted Gamma samples
curve(dgamma(x, 2, 1), 0.1, 8, col = "red", lwd = 4, add = TRUE)  # Overlay the theoretical Gamma density
curve(dexp(x, 1/10), 0.1, 8, col = "green", lwd = 4, add = TRUE)  # Overlay the instrumental distribution
N / compteur  # Calculate the efficiency of the acceptance-rejection method

# Target distribution (up to a constant factor): x * exp(-x) * (1 + cos^2(x)) ; Instrumental distribution 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)  # Sample from the instrumental Exp(1/2) distribution
    paccept <- dgamma(y, 2, 1) * (1 + cos(y)^2) /
      (dexp(y, 1/2) * ktilde)  # Acceptance probability for the target distribution
    if (runif(1) <= paccept) test <- FALSE  # Accept or reject the sample
    compteur <- compteur + 1
  }
  x[i] <- y
  test <- TRUE
}

ptilde <- function(x) dgamma(x, 2, 1) * (1 + cos(x)^2)  # Target density function
c <- integrate(ptilde, 0, 20)$val  # Compute the normalization constant by integrating the density

hist(x, prob = TRUE, nclass = 50)  # Plot the histogram of accepted samples
curve(ptilde(x) / c, col = "red", lwd = 4, add = TRUE)  # Overlay the target distribution
curve(dexp(x, 1/2), col = "green", lwd = 4, add = TRUE)  # Overlay the instrumental distribution
N / compteur * Mtilde  # Calculate the efficiency of the algorithm