Proving that $\int_{\mathbb{R}} f \ d\mu = \frac{1}{N}\sum_{i=1}^N f(\lambda_i)$ I want to know if my proof is correct and if there is some easier way to prove this, in this case I would like to see the proof.
Consider some fixed real numbers $\lambda_1\leq\ldots\leq\lambda_N$ and the measurable space $(\mathbb{R},\mathcal{B}_{\mathbb{R}},\mu)$, such that the measure $\mu:\mathcal{B}_{\mathbb{R}}\to[0,1]$ is given by $\mu(E) = \frac{1}{N}\sum_{i=1}^N\delta_{\lambda_i}(E)$. For $i=1\ldots N$, $\delta_{\lambda_i}$ is the Dirac measure.
PS: I will not deal with the set of the $\lambda_j$'s, but with the multiset of them. If some $\lambda_j$ appears more than once, this has to be considered. All the manipulations below were made with this consideration.
Let $f:\mathbb{R}\to\mathbb{R}$ be a $(\mathcal{B}_{\mathbb{R}},\mathcal{B}_{\mathbb{R}})$-measurable function such that $f\in L^1(\mu)$. I would like to prove that $$\int_{\mathbb{R}} f \ d\mu = \frac{1}{N}\sum_{i=1}^N f(\lambda_i).$$
My Proof: Let's begin showing that the result is valid for $f^+$, the proof for $f^-$ is the same. We know that $$\int_{\mathbb{R}} f^+ \ d\mu = \sup\bigg\{\int_{\mathbb{R}}\Phi \ d\mu: \ 0\leq\Phi\leq f^+, \Phi \textrm{ is simple}\bigg\}.$$
Let $\Phi = \sum_{i=1}^na_i\textbf{I}_{E_i}$, with $a_i\in \mathbb{R}$ and $\textbf{I}_{E_i}$ the indicator function of $E_i\in\mathcal{B}_{\mathbb{R}}$, then we have $$\int_{\mathbb{R}}\Phi \ d\mu = \int_{\mathbb{R}}\sum_{i=1}^na_i\textbf{I}_{E_i} \ d\mu = \sum_{i=1}^na_i\mu(E_i).$$
We can suppose, without losing generality, that $\Phi$ is in standard representation. Note that $\mu(E_i) = \frac{1}{N}|E_i\cap\{\lambda_1,\ldots,\lambda_N\}| = \frac{1}{N}\cdot(\textrm{number of }\lambda_j\textrm{'s in }E_i)$, so $$\sum_{i=1}^na_i\mu(E_i) = \frac{1}{N}\sum_{i=1}^na_i|E_i\cap\{\lambda_1,\ldots,\lambda_N\}| =\frac{1}{N}\sum_{i=1}^na_i(\textrm{number of }\lambda_j\textrm{'s in }E_i).$$ 
We don't know if each $E_i$ will contain some $\lambda_j$, it's possible that some $E_i$ contains no $\lambda_j$, in that case $a_i\mu(E_i) = 0$. With this in mind, let $E_{i_1},\ldots,E_{i_m}$ be the sets containing at least one $\lambda_j$. Note that some $E_i$ may contain more than one $\lambda_j$. The sum becomes $$\sum_{k=1}^ma_{i_k}\mu(E_{i_k}) = \frac{1}{N}\sum_{k=1}^ma_{i_k}|E_{i_k}\cap\{\lambda_1,\ldots,\lambda_N\}| =\frac{1}{N}\sum_{k=1}^ma_{i_k}(\textrm{number of }\lambda_j\textrm{'s in }E_{i_k}).$$ 
For each $\lambda_j$, we know that $\Phi(\lambda_j) = \sum_{i=1}^na_i\textbf{I}_{E_i}(\lambda_j) = a_{i_0}$ for some $1\leq i_0\leq N$. If some $E_i$ contains no $\lambda_j$, we have $a_i\textbf{I}_{E_i}(\lambda_j) = 0$ for all $\lambda_j$. It tell us that some $a_i$ may not appear in $\sum_{i=1}^N\Phi(\lambda_i)$. For what we observed before, we know that $a_{i_1},\ldots,a_{i_m}$ are the terms appearing in this sum. Note that each one of these terms may appear more than once, in fact, the number of times $a_{i_k}$ appears is the number of $\lambda_j$'s in $E_{i_k}$. Then we have that $$\frac{1}{N}\sum_{i=1}^N\Phi(\lambda_i) = \frac{1}{N}\sum_{k=1}^ma_{i_k}|E_{i_k}\cap\{\lambda_1,\ldots,\lambda_N\}| = \frac{1}{N}\sum_{k=1}^ma_{i_k}\cdot(\textrm{number of }\lambda_j\textrm{'s in }E_i).$$
Therefore, $$\int_{\mathbb{R}}\Phi \ d\mu = \frac{1}{N}\sum_{k=1}^ma_{i_k}|E_{i_k}\cap\{\lambda_1,\ldots,\lambda_N\}| = \frac{1}{N}\sum_{i=1}^N\Phi(\lambda_i).$$
So the result is valid for simple functions. Now we can write $$\int_{\mathbb{R}} f^+ \ d\mu = \sup\bigg\{\frac{1}{N}\sum_{i=1}^N\Phi(\lambda_i): \ 0\leq\Phi\leq f^+, \Phi \textrm{ is simple}\bigg\}$$
and the $\sup$ of is attained at $\frac{1}{N}\sum_{i=1}^N f^+(\lambda_i)$. Using the same reasoning we show that $\int_{\mathbb{R}} f^- \ d\mu = \frac{1}{N}\sum_{i=1}^N f^-(\lambda_i)$.
Finally, we can write $$\int_{\mathbb{R}}f \ d\mu = \int_{\mathbb{R}}f^+ \ d\mu - \int_{\mathbb{R}}f^- \ d\mu = \frac{1}{N}\sum_{i=1}^N f^+(\lambda_i) - \frac{1}{N}\sum_{i=1}^N f^⁻(\lambda_i) = $$
$$= \frac{1}{N}\sum_{i=1}^N (f^+(\lambda_i) - f^-(\lambda_i)) = \frac{1}{N}\sum_{i=1}^N f(\lambda_i).\hspace{5cm}\square$$
I have seen and used this result in the context of random matrices, but never proved  it before. So I feel like this is something necessary to know. Sorry for the long text and thank you very much for the help!
 A: Your proof seems correct. The one critique I have is that I don't see the need for the part

We don't know if each $E_i$ will contain some $\lambda_j$, it's possible that some $E_i$ contains no $\lambda_j$, in that case $a_i\mu(E_i) = 0$. With this in mind...

where you removed the terms equal to zero in the preceding sum. Those zero terms are harmless. All that getting rid of them did essentially was to add an extra paragraph to your proof and complicate the notation.

Here's how I would prove it. Let $A = \{\Lambda_1, \dotsc, \Lambda_M\}$ be the underlying set of your multiset $\{\lambda_1, \dotsc, \lambda_N\}$. Observe that the simple function
$$
\sum_{j=1}^M f(\Lambda_j) \mathbf{I}_{\{\Lambda_j\}}(x)
= \begin{cases}f(x) & x \in A \\ 0 & x \notin A\end{cases}
$$
equals $f$ almost everywhere. Therefore
$$
\int_{\mathbb R} f \, d\mu
= \int_{\mathbb R} \sum_{j=1}^M f(\Lambda_j) \mathbf{I}_{\{\Lambda_j\}} \, d\mu
= \sum_{j=1}^M f(\Lambda_j) \mu(\{\Lambda_j\})
= \frac{1}{N}\sum_{i=1}^N f(\lambda_i)
$$
as desired.
Alternatively, prove the more general statement that if $a,b \geq 0$, and $\mu,\nu$ are measures on the same $\sigma$-algebra, then $L^1(a\mu + b\nu) = L^1(\mu) \cap L^1(\nu)$ and we have
$$
\int f \, d(a\mu + b\nu) = a \int f \, d\mu + b \int f \, d\nu
$$
for all $f \in L^1(a\mu + b\nu)$.
