# program to estimate E[Z^4 cos(Z)] where Z ~ N(0,1) by simple Monte Carlo # choose simulation size n = 100 # make a vector of Z and of X Z = rnorm(n) X = Z^4 * cos(Z) # compute and output the mean and standard error m = mean(X) se = sd(X) / sqrt(n) cat("MC: ", m, " +- ", se, " (n=", n, ")", "\n", sep='') # \n is "newline" # compute an approximate 95% confidence interval cat(" 95% C.I.: (", m-1.96*se, ",", m+1.96*se, ")\n", sep='') # for comparison purposes, perform numerical integration f = function(z) { z^4 * cos(z) * exp(-z^2/2)/sqrt(2*3.14159265) } ival = integrate(f,-Inf,Inf)$value cat("Int: ", ival, "\n") # compute and output the difference between the two methods dif = m - ival cat("Dif: ", dif, "\n")