lambda <- 1

c <- exp(lambda^2/2) / lambda
g <- function(x) dexp(x, lambda)
f <- function(x) exp(-x^2 / 2)

x <- seq(0, 5, by = 0.01)
plot(x, c*g(x), type="l", log="y")
lines(x, f(x), col = "red")

n <- 10000

## half-normal

out <- NULL
num_proposal <- 0
while (length(out) < n) {
  X <- rexp(1, lambda)
  num_proposal <- num_proposal + 1
  U <- runif(1)
  if (c * g(X) * U <= f(X)) {
    out <- c(out, X) # accept X
  }
}

hist(out, probability = TRUE)
lines(x, 2/sqrt(2*pi)*f(x), col = "blue")

## standard normal

out <- NULL
while (length(out) < n) {
  X <- rexp(1, lambda)
  U <- runif(1)
  if (c * g(X) * U <= f(X)) {
    # accept X
    if (runif(1) <= 1/2) {
      out <- c(out, X)
    } else {
      out <- c(out, -X)
    }
  }
}

x <- seq(-4, 4, by = 0.01)
hist(out, probability = TRUE)
lines(x, dnorm(x), col = "blue")

## minimum of c(lambda)

ll <- seq(0.5, 2, by = 0.01)
plot(ll, 1/ll*exp(ll^2/2), type="l")

## new distribution: f(x) = cos(x)^2, x in [0, 2*pi]
f <- function(x) cos(x)^2
g <- function(x) dunif(x, min=0, max=2*pi)
c <- 2*pi

x <- seq(0, 2*pi, length.out=100)
plot(x, c * g(x), type="l", ylim=c(0,1))
lines(x, f(x))

n <- 1e5

out <- NULL
num_proposal <- 0
while (length(out) < n) {
  X <- runif(1, 0, 2*pi)
  num_proposal <- num_proposal + 1
  U <- runif(1)
  if (c * g(X) * U <= f(X)) {
    out <- c(out, X) # accept X
  }
}

hist(out, breaks = 50, probability = TRUE)
x <- seq(0, 2*pi, length.out=100)
lines(x, 1/pi*f(x), col = "blue")