In the case of a standard normal distribution, I just read that a good estimator for E[f(x)] is $\frac{1}{M}\sum_{i=1}^M f(X_i)$ (where each $X_i$ is standard normally distributed and independent). However, I don't know why. I cannot find a proof of why that approximates the integral: $$\int_{-\infty}^\infty f(x)g(x) \, dx$$
where $g(x)$ is the standard normal distribution.
Given $X$ - a continuous random variable (eg. $\mathcal{N}(0,1)$), we have $\mathbb{E}\left[f(X)\right] = \int_{-\infty}^{+\infty}f(x)g(x)dx$, where $g(x)$ is the probability $\Bbb{density}$ function (pdf) of a $\mathcal{N}(0,1)$ variable and $f$ is a measurable function of $X$.
Since you mention the cumulative density function (cdf), I will refer to that as well: given a random variable $Z$ such that $\mathbb{P}(Z\ge 0) = 1$ (meaning it is non-negative in probability) we have $\mathbb{E}[Z] = \int_{0}^{+\infty}(1-F(z))dz$, where $F(z)$ is the cdf of a $Z$-distributed random variable.
Unfortunately, a normal-distributed random variable doesn't satisfy this condition, so we cannot apply its cdf when dealing with the expected value of a normal variable.
Naive Monte Carlo is probably the most basic of probabilistic approximations of an expected value (since an expected value can be used to model many different things - scalars, vectors, matrices, functions, other random variables, name it - it is very universal in applications). The following observation is the key insight behind using Monte Carlo for approximations (I assume the notation from the OP's post):
$$ \mathbb{E}\left[\frac{1}{M}\sum_{i=1}^{M}f(X_i)\right] = \frac{1}{M}\sum_{i=1}^{M}\mathbb{E}[f(X_i)] = \frac{1}{M}\sum_{i=1}^{M}\mathbb{E}[f(X_1)] = \mathbb{E}[f(X_1)] $$
Since $\mathbb{E}[f(X_1)] = \int_{-\infty}^{+\infty}f(x)g(x)dx$, we have an equality in terms of expectation. The first equality follows from the linearity of the expectation operator, the second from the iid property of the random variables $X_1, X_2,...,X_M$ and the last from simple arithmetics.
Thus the estimator is not biased. This suggests, that if one simulates in an iid manner from certain distribution, evaluates the results, aggregates them and averages over them, then in expectation one gets the same as if one would directly compute the value.
Since simulations are subject to eg. computer related errors (we don't have a true random number generator, numbers on a computer aren't as "exact" as their mathematical ideal counterparts etc.), we have some sources of bias, that cause our computation to be an approximation rather than a direct result. For more information, please read eg. the wikipedia article on the Monte Carlo method - ploughing through it all and the related links and sources should provide you with all the knowledge you need.