# 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 = 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
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")
# plot(xlist)
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) )
# Equivalently:
# se = (1/sqrt(num)) * sqrt( ( numless*(1-numless/num)^2 +
# (num-numless)*(numless/num)^2 ) /(num-1) )
cat("P[X < -1] is about", numless/num, "\n")
cat("standard error of P[X < -1] is about", se, "\n")