N <- 1e6 # 1000000
X <- rnorm(N)

# Z <- mean(ifelse(X >= 4, 1, 0))
Z.MC <- mean(X >= 4)
RMSE.MC <- sqrt(var(X >= 4)/N)

a <- 4
Y <- rnorm(N, a)
Z.IS <- mean((Y >= 4) * dnorm(Y) / dnorm(Y, a))
RMSE.IS <- sqrt(var((Y >= 4) * dnorm(Y) / dnorm(Y, a)) / N)

Y <- rexp(N) + 4
Z.IS <- mean((Y >= 4) * dnorm(Y) / dexp(Y - 4))
RMSE.IS <- sqrt(var((Y >= 4) * dnorm(Y) / dexp(Y - 4)) / N)