In [1]:

    

rnorm_aa <- function(n, mean=0, sd=1) {
    # Generate normal distributed random numbers
    # via method of acceptance and discardment
    
    C <- pi / sqrt(2 * pi) * 2 * exp(-1/2)
    
    x_norm <- c()
    
    while (length(x_norm) < n) {
        u1 <- runif(1)
        u2 <- runif(1)
        cauchy <- tan((u1 - 0.5) * pi)
        if (u2 <= dnorm(cauchy)/(C * dcauchy(cauchy))) {
            x_norm = c(x_norm, cauchy)
        }
    }
    return(mean + x_norm * sd)
}


    

In [2]:

    

y <- rnorm_aa(500, 10, 2)
low <- floor(min(y))
high <- ceiling(max(y))


    

In [3]:

    

hist(y, breaks=seq(low, high, 1))
x <- seq(low, high, 0.1)
lines(x, 500*dnorm(x, 10, 2))


    

In [ ]: