# file "Rrej" -- a rejection sampler

pi = function(x) { dnorm(x) }  # target distribution

h = function(x) { return( x^4 ) }  # functional of interest

f = function(x) { 0.5 * exp(-abs(x)) } # double-exponential density

K = 8  # bounding constant for rejection sampler, so pi <= K f

# fsample: R routine to sample from f (using only "runif")
fsample = function() { 
   if (runif(1) < 0.5)  # probability 1/2
       return( -log(runif(1)) )  # positive Exponential r.v.
   else
       return( log(runif(1)) )  # negative Exponential r.v.
}

M = 10000  # number of attempts
hlist = NULL  # for keeping track of h function values

for (i in 1:M) {
    X = fsample()  # sample from f
    U = runif(1)  # for accept/reject
    alpha = pi(X) / (K * f(X))  # for accept/reject
    if (U < alpha)
	hlist = c(hlist, h(X))  # keep h(X) value if X accepted
}

cat("Out of", M, "attempts, obtained", length(hlist), "samples\n")

cat("mean of h is about", mean(hlist), "\n")
# plot(hlist)

se =  sd(hlist) / sqrt(length(hlist))
cat("standard error is about", se, "\n")