Assess density of a function in non-closed form I am working with Rice's Mathematical Statistics and Data Analysis and on page 179 it introduces Monte Carlo integration. This question is a smallish revision based on comments from @Iterator (in R chatroom on stackoverflow.com).
What I'm trying to show the law of large numbers
$I(f) = E[f(X)]$
where $I(f) = \frac{1}{\sqrt{2 \pi}} \int_{0}^1 e^{-\frac{x^2}{2}}$ 
and 
$f(X) = \frac{1}{n} \frac{1}{\sqrt{2 \pi}} \sum_{i=1}^{n} e^{- \frac{x^2}{2}}$.
 Note that this is a normal distribution.
I am thus trying to estimate integral of $I(f)$ on the interval $[0, 1]$, which yields
integrate(dnorm, lower = 0, upper = 1)
0.3413447 with absolute error < 3.8e-15

and compare that to $E[f(X)]$, which gives similar results
mean(sapply(1:1000, function(z, n = 1000) {
        randn <- runif(n = n, min = 0, max = 1)
        ihat <- (1/n) * (1/sqrt(2 * pi)) * sum(exp(-(randn^2)/2))
        return(ihat)
    }))
[1] 0.3413853

So I wanted to give the above function for calculating $E[f(X)]$ a field test. Density of the normal distribution on $[0, 1.96]$ is about $0.475$ (95% CI)
integrate(dnorm, lower = 0, upper = 1.96)
0.4750021 with absolute error < 5.3e-15

times $2$ gives roughly $0.95$. Check!
When I try to apply similar logic for $f(X)$, my function squats and returns something else
mean(sapply(1:1000, function(z, n = 1000) {
        randn <- runif(n = n, min = 0, max = 1.96)
        ihat <- (1/n) * (1/sqrt(2 * pi)) * sum(exp(-(randn^2)/2))
        return(ihat)
    }))
[1] 0.24224

I'm not exactly sure what's going on here. Any tips on what am I doing wrong, failing to get ~$0.475$ for the last expression?
 A: It’s a bit late (Henry’s answer has been accepted!), but I’ll aggregate my comments in one answer, with one additional comment on the use of standard error.


*

*The simplest thing to integrate is a constant function $f(x) = 1$, let’s say on $[0, 1.96]$. You will sample $n=1000$ values $x_1,\dots,x_n$ uniformly in this interval, and compute the mean of the $f(x_i)$, which is 1. But you want to compute $\int_0^{1.96} 1 \mathrm dx = 1.96$... 

*In fact you’ve been computing [an approximation of]
$$ E\left( f(X) \right) = \int_{-\infty}^{+\infty} f(x) \phi(x) \mathrm dx,$$
where $\phi(x)$ is the density of the sampling distribution; here $\phi(x)$ is 
$$\phi(x) = \left\{ \begin{array}{ll}0 &\text{if } x\notin [0,1.96] \\
{1\over 1.96} &\text{if } x\in [0,1.96], \end{array}\right.$$
hence you’ve been computing ${1\over 1.96} \int_0^{1.96} f(x) \mathrm dx$. This shows you how to rectify your computation, just multiply by 1.96.

*Consider the random variables $Y_i = f(X_i)$. They are independent, identically distributed; by the central limit theorem, their mean (your estimation of $E(Y)$ which is, up to multiplication by a constant, the quantity of interest) are normally distributed. You can thus get a confidence interval on $E(Y)$ using the sample standard error 
$s$: the standard deviation of the mean is $\simeq s/\sqrt n$.
Illustration with a piece of R code (for 95% CI).
> x <- runif(10000, 0, 1.96)
> y <- dnorm(x)
> 1.96*mean(y)
[1] 0.4744983
> 1.96*(mean(y) + c(-1,1)*1.96*sd(y)/sqrt(length(y)) )
[1] 0.4701337 0.4788629

A: Your code returns an estimate of the average density over the interval.  The average density is not the integral, unless the interval is of length 1 as in the first example.
To emphasis Elvis's point, note 
> 0.24224*1.96
[1] 0.4747904

