Why does convolution work? So I know that if we want to find the probability distribution of a sum of independent random variables $X + Y$, we can compute it from the probability distributions of $X$ and $Y$, by saying
$$f_{X + Y}(a) = \int_{x = -\infty}^{\infty} f_{X, Y}(X = x, Y = a - x)~dx = \int_{x = -\infty}^{\infty} f_X(x) f_Y(a - x)~dx$$
Intuitively, this makes sense, because if we want to find the probability that two random variables sum to $a$, it's basically the sum of the probabilities of all the events that lead to those variables summing to $a$.  But how can I formally prove this statement?
 A: The more general solution considers $Z = X + Y$ where $X$ and $Y$ are not necessarily independent. A common solution strategy for problems where you are wondering where a PDF came from or how to justify it, is to find a cumulative probably instead, then differentiate to reduce the CDF to a PDF. 
It is quite easy to see that in that case $F_Z(z) = \mathrm{P}(Z \leq z) = \int \int_R f_{X,Y}(x,y)\,\mathrm{d}x\,\mathrm{d}y$ where $R$ is the region of the $x$-$y$ plane for which $x + y \leq z$.
This is the blue-hatched region in the diagram below. It's natural to integrate over this region by breaking it down into strips - I have done it with vertical strips but horizontal ones will do. Effectively I end up with a strip for each $x$ co-ordinate, ranging from $-\infty$ to $\infty$, and along each strip I want $y$ values not to rise above the line $x + y = z$, so $y \leq z - x$. 

Now we have obtained limits of integration in terms of $x$ and $y$, we can make a substitution $u=x$, $v=x+y$ as follows, with the aim of getting $z$ to appear as the upper limit of $v$. The maths is straightforward so long as you understand the use of the Jacobian to change the variables.
$$F_Z(z) = \int_{x = -\infty}^{x=\infty}\int_{y=-\infty}^{y=z-x}f_{X,Y}(x,y)\,\mathrm{d}x\,\mathrm{d}y = \int_{v = -\infty}^{v=z}\int_{u=-\infty}^{y=\infty}f_{X,Y}(u,v-u)\,\mathrm{d}u\,\mathrm{d}v $$
So long as certain conditions are met we can differentiate under the integral sign with respect to $z$ to obtain:
$$f_Z(z) = \int_{-\infty}^{\infty}f_{X,Y}(u, z-u)\,\mathrm{d}u$$
That works even if $X$ and $Y$ are not independent. But if they are, we can rewrite the joint density as the product of the two marginal ones:
$$f_Z(z) = \int_{-\infty}^{\infty}f_X(u)f_Y(z-u)\,\mathrm{d}u$$
The dummy variable $u$ can without harm be written as $x$ if desired.
My notation for the integrals exactly follows Section 6.4 of Geoffrey Grimmett and Dominic Walsh, Probability: An Introduction, Oxford University Press, New York, 2000.
A: The statement is true if and only if the right hand side acts like a density for $X+Y$; that is, 
$$F_{X+Y}(a)=\mathbb{P}(X+Y\le a) = \int_{-\infty}^a f_{X+Y}(z)\,\mathrm{d}z = \int_{-\infty}^a \left(\int f_X(x) f_Y(z-x)\,\mathrm{d}x\right)\mathrm{d}z$$
for all $a$.  Let's verify this by starting with the right hand side.
Apply Fubini's Theorem to change the order of integration and make the substitution $z = x+y$.  The determinant of its Jacobian is $1$, so no additional terms are introduced by this change of variables.  Note that because $z$ and $y$ are in one-to-one correspondence and $-\infty \lt z \le a$ if and only if $-\infty \lt y \lt a-x$, we may rewrite the integral as
$$=\int \left(\int_{-\infty}^{a-x}f_X(x)f_Y(y)\,\mathrm{d} y\right)\mathrm{d}x.$$
By definition this is the integral over $\mathbb{R}^2$ of
$$=\iint I(x+y\le a)f_X(x)f_Y(y)\,\mathrm{d}y\,\mathrm{d}x$$
where $I$ is the indicator function of a set.  Finally, since $X$ and $Y$ are independent, $f_{(X,Y)}(x,y) = f_X(x)f_Y(y)$ for all $(x,y)$, revealing the integral as merely the expectation
$$=\iint I(x+y\le a)f_{(X,Y)}(x,y)\,\mathrm{d}y\,\mathrm{d}x = \mathbb{E}(I(X+Y\le a))=\mathbb{P}(X+Y\le a),$$
as desired.

More generally, even when one or both of $X$ or $Y$ does not have a distribution function, we can still obtain
$$F_{X+Y}(a) = \mathbb{E}_X\left(F_Y(a-X)\right) = \mathbb{E}_Y\left(F_X(a-Y)\right)$$
directly from basic definitions, using the expectation of indicators to go back and forth between probabilities and expectations and exploiting the independence assumption to break the calculation into separate expectations with respect to $X$ and $Y$:
$$\eqalign{
\mathbb{P}(X+Y\le a) &= \mathbb{E}(I(X+Y\le a)) \\
&= \mathbb{E}_X\left(\mathbb{E}_Y(I(X+Y\le a)\right) \\
&= \mathbb{E}_X\left(\mathbb{P}_Y(Y\le a-X)\right) \\
&=\mathbb{E}_X(F_Y(a-X)).
}$$
This includes the usual formulas for discrete random variables, for instance, albeit in a slightly different form than usual (because it is stated in terms of the CDFs rather than the probability mass functions).
If you have a strong enough theorem about interchanging derivatives and integrals, you can differentiate both sides with respect to $a$ to obtain the density $f_{X+Y}$ in one stroke,
$$\eqalign{
f_{X+Y}(a) &= \frac{\mathrm{d}}{\mathrm{d}a} F_{X+Y}(a) =\mathbb{E}_X\left(\frac{\mathrm{d}}{\mathrm{d}a}  F_Y(a-X)\right) = \mathbb{E}_X \left(f_Y(a-X)\right) \\
&= \int f_X(x) f_Y(a-x) \,\mathrm{d} x.
}$$
