mc_int.R

library(cowplot)
f <- function(x) {
  exp(-x^2) + sin(10*x)
}
n <- 3e3
bounds <- c(0, 1)
set.seed(123)
x <- runif(n, bounds[1], bounds[2])
y <- f(x)
p1 <- qplot(x, y, geom = "line")
val <- integrate(f, bounds[1], bounds[2])$value
mcint <- diff(bounds) * (cumsum(y) / (1:n))
mcse <- sqrt(cumsum((y - mcint)^2)) / (1:n)
df <- data.frame(n = 1:n, mcint, mcse)
p2 <- ggplot(df, aes(n, mcint)) +
  geom_line(size = 0.3) +
  geom_hline(yintercept = val, color = "red", size = 0.3) +
  geom_ribbon(aes(
    ymin = mcint - 2 * mcse,
    ymax = mcint + 2 * mcse
  ),
  fill = "grey70",
  alpha = 1/2) + ylim(0.5, 1.25) +
  xlab("iteration") + ylab("integral value")
plot_grid(p1, p2)