approx_p <- function(n0, a, b) {
    x <- rnorm(n0)
    sum(x >= a & x <= b) / n0
}
# berechne wahres p durch Integration
p <- integrate(dnorm, lower=-1, upper=3)$value
p
# berechne Differenz zwischen Näherung und wahrem Wert für verschiedene n0
abs(p - approx_p(1e1, -1, 3))
abs(p - approx_p(1e3, -1, 3))
abs(p - approx_p(1e5, -1, 3))