# auxiliary variable rejection sampling example g = function(y) { y^3 * sin(y^4) * cos(y^5) }; M = 10^6; xlist = runif(M); ylist = runif(M); pilist = xlist[ylist < g(xlist)]; len = length(pilist); cat("Obtained", len, "samples out of", M, "attempts.\n"); hlist = pilist^2; cat("E(h) estimate: ", mean(hlist), "+-", sd(hlist)/sqrt(len), "\n");