Show that the sum of two random variables is a mixture Take any  $({\lambda},{\mu},F,G)$ such that
1) $\lambda\equiv (\lambda_1,..., \lambda_J)$, $\lambda_j\in (0,1)$ for each $j=1,...,J$ and $\sum_{j=1}^J \lambda_j=1$
2) $\mu\equiv (\mu_1,..., \mu_J)$ with $\mu_1<...<\mu_J$.
3) $F$ is a cumulative distribution function with steps of height $\lambda_j$ at each $\mu_j$. 
4) $G$ is a cumulative distribution function whose associated probability mass function of probability density function is symmetric around zero.
Consider the mutually independent random variables, $Y\stackrel{d}{\sim}   F$ and $Z\stackrel{d}{\sim} G$.
Question: Show that the cumulative distribution function of  $Y+Z$ is the mixture $
\sum_{j=1}^{J} \lambda_j G(x-\mu_j)$ at each $ x \in \mathbb{R}$. 
I have done some simulations and realised that the result seems in fact to hold. However, when I try to prove it formally, I'm completely stuck. Could you help? Even some informal intuition would be very useful.
 A: This ultimately is Fubini's Theorem, but to keep the analysis elementary let's stick to a finite mixture of arbitrary distributions.  By definition, this means $Y$ can be considered in terms of other variables $Y_j$ with cumulative distribution functions (CDFs) $F_j$ and that for any number $y,$
$$F(y) = \Pr(Y \le y) = \sum_{j} \lambda_j \Pr(Y_j \le y) = \sum_j \lambda_j F_j(y).$$
Now let $x$ be any number and compute the CDF of $Y+Z$ at $x$ in terms of the CDF $G$ of $Z$ (defined by $G(z) = \Pr(Z \le z)$) as
$$\eqalign{
{\Pr}_{Y,Z}(Y+Z \le x) &= {\Pr}_{Y,Z}( Z \lt x - Y ) \\
&= \mathbb{E}_Y (G(x - Y)) \\
&= \int G(x-y)\,\mathrm{d}F(y) \\
&= \int G(x-y)\, \mathrm{d}\left(\sum_j \lambda _j F_j(y)\right) \\
&= \sum_j \lambda_j \int G(x-y)\,\mathrm{d}F_j(y).\tag{*}
}$$
(The switching of the order of integration and summation merely expressed the linearity of integration but can be seen as an example of Fubini's Theorem.)
In the question itself, the component distributions are atoms at the $\mu_j$ and their CDFs $F_j,$ which jump from a value of $0$ to a value of $1$ at $\mu_j,$ have the property that for any piecewise differentiable function $H$ with ${\lim}_{x\to\infty}H(x)=0,$
$$\eqalign{
\int H(y)\,\mathrm{d}F_j(y) &= H(y)F_j(y)\mid_{-\infty}^\infty - \int H^{\prime}(y) F_j(y)\,\mathrm{d}y \\
&=(0 - 0) - \left(\int_{-\infty}^{\mu_j}H^\prime(y)(0)\mathrm{d}y + \int_{\mu_j}^\infty H^\prime(y)(1)\mathrm{d}y\right) \\
&= (0 - 0) - \left(0 + (0 - H(\mu_j))\right) \\
&= H(\mu_j).
}$$
The integrals are all in the sense of Riemann or Lebesgue and the initial equality is integration by parts.  The zeros arise from the assumed limiting value of $H$ at $\infty$ and from the fact that $F_j$ is identically zero for very negative arguments.
Consequently, applying this to the function $H: y \to G(x-y)$ (whose limit as $y\to\infty$ clearly is zero), the general result $(*)$ reduces to
$$\Pr(Y + Z \le x) = \sum_j \lambda_j H(\mu_j) = \sum_j \lambda_j G(x-\mu_j),$$
QED.
Note $G$ does not have to be symmetric, but it is important that it be piecewise differentiable.  This includes the CDFs of continuous random variables, the usual discrete random variables (with no accumulation points among their supports), and mixtures thereof.
