# program to estimate E[Z^4 cos(Z^3)] where Z ~ N(0,1) by simple Monte Carlo

# choose simulation size
n = 100000

# make a vector of Z and of X
Z = rnorm(n)
X = Z^4 * cos(Z^3)

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

# for comparison purposes, perform numerical integration
f = function(z) { z^4 * cos(z^3) * exp(-z^2/2)/sqrt(2*pi) }
ival = integrate(f,-Inf,Inf)
cat("Int: ")
print(ival)

# compute and output the difference between the two methods
dif = m - ival$value
cat("Dif: ", dif, "\n")