# 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
xlist = hlist = fulllist = NULL  # for keeping track of sampled 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
    fulllist = c(fulllist, X)  # list of *all* the values
    if (U < alpha) {
	xlist = c(xlist, X)  # keep X value if accepted
	hlist = c(hlist, h(X))  # keep h(X) value if accepted
    }
}

cat("Out of", M, "attempts, obtained", length(xlist), "samples\n")
cat("mean of X is about", mean(xlist), "\n")

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

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

num = length(xlist);
numless = length(xlist[xlist < -1]);
se = (1/sqrt(num)) * sqrt( numless*(num-numless)/num / (num-1) )
cat("P[X < -1] is about", numless/num, "\n")
cat("standard error of P[X < -1] is about", se, "\n")

plot( fulllist, rep(0,length(fulllist)), pch="|" )
points(xlist, rep(0,length(xlist)), col="red", pch="|")