# 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

# doubleexp: R routine to sample from double-exponential
doubleexp = 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
xlist = hlist = fulllist = rep(NA, M) # for keeping track of sampled values
numsamples = 0

for (i in 1:M) {
    X = doubleexp()  # sample from f
    U = runif(1)  # for accept/reject
    alpha = pi(X) / (K * f(X))  # for accept/reject
    fulllist[i] = X  # list of *all* the values
    if (U < alpha) {
	xlist[i] = X  # keep X value if accepted
	hlist[i] = h(X)  # keep h(X) value if accepted
	numsamples = numsamples + 1
    }
}

cat("Out of", M, "attempts, obtained", numsamples, "samples\n")
cat("mean of X is about", mean(xlist, na.rm=TRUE), "\n")

se =  sd(xlist, na.rm=TRUE) / sqrt(numsamples)
cat("standard error of X is about", se, "\n")

cat("mean of h(X) is about", mean(hlist, na.rm=TRUE), "\n")
se =  sd(hlist, na.rm=TRUE) / sqrt(numsamples)
cat("standard error of h(X) is about", se, "\n")

plot( fulllist, rep(0,M), pch="|", ylim=c(-1,2) )
points(xlist, rep(0.2,M), col="red", pch="|")