sir-aljo.r

rz_helper = function() {
    y1 = rexp(1, 1) # step 1
    y2 = rexp(1, 1) # step 2

    # step 3
    while (y2 <= ((y1 - 1)^2)/2) {
        y1 = rexp(1, 1)
    }

    # step 4
    u = runif(1)

    if (u <= 0.5) {
        return(y1)
    } else {
        return(-y1)
    }
}

rz = function(n) {
    out = c()
    for (i in 1:n) {
        out[i] = rz_helper()
    }
    return(out)
}

# Generate 10000 standard normal random samples
s = rz(10000)

# Take the histogram
hist(s)

# Check the mean if close to 0
mean(s)

# Check the variance if close to 1
var(s)