n = 100000 Z = rnorm(n) X = Z^4 * cos(Z^3) m = mean(X) se = sd(X) / sqrt(n) cat("MC: ", m, " +- ", se, " (n=", n, ")", "\n", sep='') f = function(z) { z^4 * cos(z^3) * exp(-z^2/2)/sqrt(2*pi) } ival = integrate(f,-Inf,Inf) cat("Int: ") print(ival) dif = m - ival$value cat("Dif: ", dif, "\n")