ndraw <- 400

z <- rnorm(ndraw,mean=5,sd=3)

z2 <- subset(z,z>=0 & z <= 2)

print(sum(z2)/ndraw)

fnorm <- function(z) 1/sqrt(2*pi*9)*exp(-(z-5)^2/(2*9))

integrate(fnorm,lower=0,upper=2)