# Generate values from the half normal distribution (|N(0,1)|) by an 
# Acceptance-Rejection scheme with an Exponential(mu) majorizing distribution
#
# n = number of values to be generated
# mu = mean of exponential distribution

rhalfnorm <- function(n, mu=1) {
  y <- rep(0,n)
  c <- mu * exp(0.5 / mu^2) * sqrt(2/pi)

  for(i in 1:n) {
    accept <- F
    while (!accept) {
      x <- rexp(1,mu)
      r <- exp(-(x - 1/mu)^2 / 2)
      accept <- (runif(1) <= r)
    }
    y[i] <- x
  }
  y
}