# 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")