[proofplan]
Set $Y_i := f(X_i)$ and reduce the statement to the elementary variance computation for the average of independent, identically distributed square-integrable real random variables. The square-integrability assumption $\int_E f(x)^2\,d\pi(x)<\infty$ guarantees that the mean and variance of each $Y_i$ are finite. Linearity of expectation gives unbiasedness, independence removes all cross terms in the variance of the average, and the definition of root mean square error then gives the stated $N^{-1/2}$ rate.
[/proofplan]
[step:Check that the sampled function values have finite mean and variance]
For each $i\in\mathbb N$, define the real-valued [random variable](/page/Random%20Variable) $Y_i:\Omega\to\mathbb R$ by
\begin{align*}
Y_i(\omega) := f(X_i(\omega)).
\end{align*}
Since $X_i$ is $\mathcal F/\mathcal E$-measurable and $f$ is $\mathcal E/\mathcal B(\mathbb R)$-measurable, $Y_i$ is a real-valued random variable. Since $X_i$ has law $\pi$,
\begin{align*}
\mathbb E[Y_i^2] = \int_E f(x)^2\,d\pi(x) < \infty.
\end{align*}
Since $\pi(E)=1$, the [Cauchy-Schwarz inequality](/theorems/432) applied to the real functions $|f|$ and $1$ on the probability space $(E,\mathcal E,\pi)$ gives
\begin{align*}
\int_E |f(x)|\,d\pi(x) \leq \left(\int_E f(x)^2\,d\pi(x)\right)^{1/2}\left(\int_E 1^2\,d\pi(x)\right)^{1/2} = \left(\int_E f(x)^2\,d\pi(x)\right)^{1/2}.
\end{align*}
Thus $Y_i$ is integrable and square-integrable with respect to $\mathbb P$, and
\begin{align*}
\mathbb E[Y_i] = \int_E f(x)\,d\pi(x) = I.
\end{align*}
Consequently $\operatorname{Var}(Y_i)$ is finite for every $i\in\mathbb N$.
[/step]
[step:Use linearity of expectation to prove unbiasedness]
Since each $Y_i$ is integrable and $\widehat I_N = N^{-1}\sum_{i=1}^N Y_i$, linearity of expectation gives
\begin{align*}
\mathbb E[\widehat I_N] = \frac{1}{N}\sum_{i=1}^N \mathbb E[Y_i].
\end{align*}
Using $\mathbb E[Y_i]=I$ for every $i\in\{1,\dots,N\}$, we obtain
\begin{align*}
\mathbb E[\widehat I_N] = \frac{1}{N}\sum_{i=1}^N I = I.
\end{align*}
Therefore $\widehat I_N$ is unbiased for $I$.
[/step]
[step:Use independence to compute the variance of the sample average]
Define the centered random variables $Z_i:\Omega\to\mathbb R$ by
\begin{align*}
Z_i := Y_i-I.
\end{align*}
Then $\mathbb E[Z_i]=0$ and $\operatorname{Var}(Y_i)=\mathbb E[Z_i^2]$. Since the random variables $X_1,\dots,X_N$ are independent, the random variables $Y_1,\dots,Y_N$ are independent, hence $Z_1,\dots,Z_N$ are independent. For $i\neq j$, independence and integrability give
\begin{align*}
\mathbb E[Z_iZ_j] = \mathbb E[Z_i]\mathbb E[Z_j] = 0.
\end{align*}
Therefore
\begin{align*}
\operatorname{Var}(\widehat I_N) = \mathbb E\left[\left(\widehat I_N-I\right)^2\right].
\end{align*}
Since $\widehat I_N-I = N^{-1}\sum_{i=1}^N Z_i$, expanding the square gives
\begin{align*}
\operatorname{Var}(\widehat I_N) = \frac{1}{N^2}\mathbb E\left[\sum_{i=1}^N\sum_{j=1}^N Z_iZ_j\right].
\end{align*}
Linearity of expectation is valid because each $Z_iZ_j$ is integrable: the Cauchy-Schwarz inequality for real random variables on $(\Omega,\mathcal F,\mathbb P)$ gives $\mathbb E[|Z_iZ_j|]\leq (\mathbb E[Z_i^2])^{1/2}(\mathbb E[Z_j^2])^{1/2}<\infty$. Hence
\begin{align*}
\operatorname{Var}(\widehat I_N) = \frac{1}{N^2}\sum_{i=1}^N\sum_{j=1}^N \mathbb E[Z_iZ_j].
\end{align*}
The off-diagonal terms vanish, and the diagonal terms equal $\operatorname{Var}(Y_i)$. Since the $Y_i$ are identically distributed,
\begin{align*}
\operatorname{Var}(\widehat I_N) = \frac{1}{N^2}\sum_{i=1}^N \operatorname{Var}(Y_i) = \frac{\operatorname{Var}(Y_1)}{N}.
\end{align*}
[guided]
The variance computation is where both independence and identical distribution are used. Define $Z_i:\Omega\to\mathbb R$ by
\begin{align*}
Z_i := Y_i-I.
\end{align*}
These are the errors of the individual samples around the target mean. Since $\mathbb E[Y_i]=I$, each $Z_i$ is centered:
\begin{align*}
\mathbb E[Z_i]=0.
\end{align*}
Also,
\begin{align*}
\operatorname{Var}(Y_i)=\mathbb E[(Y_i-\mathbb E[Y_i])^2]=\mathbb E[Z_i^2].
\end{align*}
The estimator error is the average of these centered errors:
\begin{align*}
\widehat I_N-I = \frac{1}{N}\sum_{i=1}^N Z_i.
\end{align*}
Because $\widehat I_N$ is unbiased, its variance is the second moment of this centered error:
\begin{align*}
\operatorname{Var}(\widehat I_N)=\mathbb E[(\widehat I_N-I)^2].
\end{align*}
Substituting the expression above gives
\begin{align*}
\operatorname{Var}(\widehat I_N)=\frac{1}{N^2}\mathbb E\left[\left(\sum_{i=1}^N Z_i\right)^2\right].
\end{align*}
Expanding the square produces diagonal terms $Z_i^2$ and cross terms $Z_iZ_j$:
\begin{align*}
\operatorname{Var}(\widehat I_N)=\frac{1}{N^2}\mathbb E\left[\sum_{i=1}^N\sum_{j=1}^N Z_iZ_j\right].
\end{align*}
The expectation may be passed through the finite sum by linearity. Each product $Z_iZ_j$ is integrable because the Cauchy-Schwarz inequality for real random variables on $(\Omega,\mathcal F,\mathbb P)$ gives $\mathbb E[|Z_iZ_j|]\leq (\mathbb E[Z_i^2])^{1/2}(\mathbb E[Z_j^2])^{1/2}<\infty$. Hence
\begin{align*}
\operatorname{Var}(\widehat I_N)=\frac{1}{N^2}\sum_{i=1}^N\sum_{j=1}^N\mathbb E[Z_iZ_j].
\end{align*}
Now independence removes the cross terms. Since $X_1,\dots,X_N$ are independent, the [measurable functions](/page/Measurable%20Functions) $Y_1,\dots,Y_N$ are independent, and subtracting the constant $I$ preserves independence, so $Z_1,\dots,Z_N$ are independent. Thus, for $i\neq j$,
\begin{align*}
\mathbb E[Z_iZ_j]=\mathbb E[Z_i]\mathbb E[Z_j]=0.
\end{align*}
Only the diagonal terms remain:
\begin{align*}
\operatorname{Var}(\widehat I_N)=\frac{1}{N^2}\sum_{i=1}^N\mathbb E[Z_i^2].
\end{align*}
Because the $X_i$ have the same law $\pi$, the random variables $Y_i=f(X_i)$ have the same distribution, so all variances $\mathbb E[Z_i^2]$ equal $\operatorname{Var}(Y_1)$. Therefore
\begin{align*}
\operatorname{Var}(\widehat I_N)=\frac{1}{N^2}\sum_{i=1}^N\operatorname{Var}(Y_1)=\frac{\operatorname{Var}(Y_1)}{N}.
\end{align*}
[/guided]
[/step]
[step:Convert the variance formula into the root mean square error formula]
By definition,
\begin{align*}
\operatorname{RMSE}(\widehat I_N)=\left(\mathbb E[(\widehat I_N-I)^2]\right)^{1/2}.
\end{align*}
Since $\mathbb E[\widehat I_N]=I$, the mean square error equals the variance:
\begin{align*}
\mathbb E[(\widehat I_N-I)^2]=\operatorname{Var}(\widehat I_N).
\end{align*}
Using the variance formula from the previous step and $Y_1=f(X_1)$, we obtain
\begin{align*}
\operatorname{RMSE}(\widehat I_N)=\left(\frac{\operatorname{Var}(Y_1)}{N}\right)^{1/2}=\frac{\sqrt{\operatorname{Var}(f(X_1))}}{\sqrt{N}}.
\end{align*}
This is the claimed crude Monte Carlo root mean square error rate.
[/step]
