## A minimalist proof of the linearity of E

Let's dig into the inexplicit steps of this calculation of $\mathsf E[X+Y] = \mathsf E[Z]$!

\begin{align}
&\mathsf E[Z]\\
&=\int_\Bbb R z~p(z) ~ \mathsf d z && \text{by definition of the expected value}\\
&=\int_\Bbb R z ~ \Big( \int_\Bbb R p(z \cap y) ~ dy \Big) \mathsf d z && \text{by the law of total probability}\\
&=\int_\Bbb R z ~ \Big( \int_\Bbb R p(z | y) ~ p(y) ~ dy \Big) \mathsf d z && \text{by definition of the conditional probability}\\
&=\int_\Bbb R z ~ \Big( \int_\Bbb R p(z - y | y) ~ p(y) ~ dy \Big) \mathsf d z && \text{since p(Z=x+y|Y=y) = p(X=x|Y=y)}\\
&=\int_\Bbb R \Big( \int_\Bbb R z ~ p(z - y | y) ~ p(y) ~ dy \Big) \mathsf d z && \text{by linearity of the integral}\\
&=\int_\Bbb R \Big( \int_\Bbb R (x + y) ~ p(x | y) ~ p(y) ~ dy \Big) \mathsf d x && \text{by substituting $x = z - y$ (note that $dx = dz$)}
\end{align}

By linearity of the integral, this last line is equal to
$$\int_\Bbb R x \Big( \int_\Bbb R p(x | y) p(y) dy \Big) \mathsf d x
+ \int_\Bbb R \Big( \int_\Bbb R y ~ p(x | y) p(y) dy \Big) \mathsf d x$$

Let's deal with the first term:
\begin{align}
\int_\Bbb R x \Big( \int_\Bbb R p(x | y) p(y) dy \Big) \mathsf d x
&=\int_\Bbb R x \Big( \int_\Bbb R p(x \cap y) dy \Big) \mathsf d x && \text{conditional prob}\\
&=\int_\Bbb R x p(x) \mathsf d x && \text{total prob}\\
&=\mathsf E[X]
\end{align}

Let's deal with the second term. We'll make the reasonable assumption that X and Y are measurable so as to use [Fubini's theorem](https://en.wikipedia.org/wiki/Fubini%27s_theorem). This theorem allows us to change the order of integrals:
\begin{align}
\int_\Bbb R \Big( \int_\Bbb R y ~ p(x | y) ~ p(y) ~ dy \Big) \mathsf d x
&=\int_\Bbb R y \Big( \int_\Bbb R p(y) ~ p(x | y) ~ dx \Big) \mathsf d y &&\text{with Fubini's theorem}\\
&=\int_\Bbb R y \Big( \int_\Bbb R p(x \cap y) ~ dx \Big) \mathsf d y &&\text{conditional prob}\\
&=\int_\Bbb R y ~ p(y) ~ \mathsf d y &&\text{total prob}\\
&=\mathsf E[Y]
\end{align}

## About $dz$

> Because in my point of view I may write
\begin{align}\mathsf E(X+Y)&=\int_\Bbb R (x+y)~f_{X+Y}(x+y)\mathsf d(x+y)\\&=  \int_\Bbb R (x+y)~f_{X+Y}(x+y)\mathsf (dx+dy) \\&=  \int_\Bbb R ~f_{X+Y}(x+y) (xdx+ydx+xdy+ydy)\\&\vdots \end{align}
> 
> Can someone explain this? Especially second step suddenly introduce one more variable becoming 2 dimensional integral and the third step replace  𝑓𝑋+𝑌,𝑌(𝑧,𝑦)
 to  𝑓𝑋,𝑌(𝑧−𝑦,𝑦)
. I am quite confused of the change of dz -> dzdy -> dxdy, which seems like change freely...

Indeed, as you correctly stated, one cannot change the order of these notations freely, nor can they separate $d(x+y)$ into $dx + dy$.

In my opinion, you can think of $dz$ as the part of these notations that indicates what is the name of the variable that "scans" your support. Since you are integrating once over a given support or space (in this case $\mathbb{R}$), $d(x+y)$ is a confusing notation which gives you the impression that you could split it in two ($dx + dy$), when in fact it is not the case.

When substituting a variable for another, it's in fact a _function_ of the variable that you are substituting. Check the wikipedia article about substitution for what I think is a clear formulation: https://en.wikipedia.org/wiki/Integration_by_substitution#Statement_for_definite_integrals

In our case, we are substituting $x$ to $\phi(z) = z - y$, and $dx$ is equal to $\frac{d\phi}{dz}(z) ~ dz = 1 \cdot dz$.