Let $X \sim N(\mu_1, \sigma_1^2)$ and $Y \sim N(\mu_2, \sigma_2^2)$ be independent. It can easily be shown via moment generating functions that $Z=X+Y$ is also normally distributed with $Z\sim N(\mu_1+\mu_2, \sigma_1^2+\sigma_2^2)$.
I was wondering if one could instead use the change-of-variables formula to derive the distribution of $Z$. The obvious way seems to be to go from the vector $(X,Y)$ to the vector $(Z,Y) = (X+Y, Y)$. The Jacobian of this is 1 and the result seems at first to be straightforward.
However the result of course is the joint distribution of $Z$ and $Y$, so to get the distribution of $Z$ one has to then marginalise out $y$ by integration.
It is easy to see from the result that in the general case, i.e. Gaussian or not, the sum of independent random variables is distributed as the convolution (see also Wikipedia). In the Gaussian case, how would one go about computing the required integral? That is, the integral
$$\int \frac{1}{\sqrt{2\pi}} \frac{1}{\sigma_1} \exp(-\frac{1}{2} \frac{((z-y)-\mu_1)^2}{\sigma_1^2} ) \frac{1}{\sqrt{2\pi}} \frac{1}{\sigma_2} \exp(-\frac{1}{2} \frac{(y-\mu_2)^2}{\sigma_2^2} d y \,?$$
Or, instead of $(X,Y) \rightarrow (X+Y,Y)$, is there another change or variables that would result in an easier integral?
Best Answer
The technique used to evaluate that integral is called completing the square in the exponent which is pretty close to what you learned (and undoubtedly soon forgot) in high school when you derived the formula for the roots of a quadratic polynomial. Multiply the exponential functions that are in the integrand by using the identity $\exp(a+b) = \exp(a)\exp(b)$ bass ackwards to get $\exp\left(-\frac{1}{2\alpha} (y^2 - 2\beta y +\gamma)\right)$ where $\beta$ is a linear function of $z$, $\gamma$ is a quadratic function of $z$ and $\alpha$ does not depend on $z$ at all. Then, massage the argument of the exponent to write it as \begin{align} y^2 - 2\beta y +\gamma &= y^2 - 2\beta y + \beta^2 + \gamma - \beta^2\\ &= (y-\beta)^2 + (\gamma - \beta^2) \end{align} Consequently, the integral becomes of the form $$\frac{1}{2\pi \sigma_1\sigma_2} \exp\left(-\frac{1}{2\alpha} (\gamma-\beta^2)\right) \int_{-\infty}^\infty \exp\left(-\frac{1}{2\alpha} (y-\beta^2)\right) \, \mathrm dy$$ where the integrand needs to be recognized as almost the density function of a $N(\beta, \alpha)$ random variable giving $$f_Z(z) = A \exp\left(-\frac{1}{2\alpha} (\gamma-\beta^2)\right) $$ for specific value of $A$ that I will leave it for you to work out. Notice that the argument of the exponential is a quadratic function of $z$ and so $Z$ has a Gaussian density. As Huber remarks in a comment, somewhere along such a calculation it might dawn on you that the density of $Z$ is going to end up being a Gaussian density at which point you can save yourself a lot of algebra by working out the mean and variance of $Z$ and plugging in results into the Gaussian density formula. And no, I don't mean plugging in $E[Z]$ etc into the integrand of the convolution integral; the integrand of the convolution integral has nothing to do with $Z$.
If one is not all set on using the convolution integral formula and nothing else must be used, the easiest way of directly determining the pdf of $Z=X+Y$ where $X$ and $Y$ are independent Gaussian random variables is to find the CDF first, recognize it as a Gaussian CDF and hence write down the density of $Z$. The argument is as follows.
Consider independent random variables $A, B \sim N(0,1)$ and set $$Z = (\sigma_1 A + \mu_1) + (\sigma_2 B + \mu_2) = X + Y$$ where $X\sim N(\mu_1, \sigma_1^2), Y \sim N(\mu_2, \sigma_2^2)$ are independent random variables. Then, $$F_Z(z) = P\{\sigma_1 A + \sigma_2 B + \mu_1 + \mu_2)\leq z\} =\iint_{\sigma_1a + \sigma_2b + (\mu_1 + \mu_2)\leq z} \phi(a)\phi(b) \,\mathrm da \mathrm db$$ where $\phi(\cdot)$ is the standard Gaussian density. But the integral has circular symmetry and so the probability depends only on the distance of the line $\sigma_1 a + \sigma_2 b + (\mu_1 + \mu_2) = z$. Indeed, by a rotation of coordinates we can make this line vertical (say), and if it crosses the $a$ axis at $d$, then $$F_Z(z) = P\{A \leq d\} = \Phi(d)$$ showing that $Z$ is a Gaussian random variable. At this point, as whuber says, we can compute $$E[Z] = E[X+Y] = \mu_1 + \mu_2, \operatorname{var}(Z) =\operatorname{var}(X+Y) = \operatorname{var}(X) + \operatorname{var}(Y) = \sigma_1^2+\sigma_2^2$$ and write down the density of $Z$. Or, we can bull ahead and note that standard results in analytical geometry tell us that $$d = \left(\frac{z - (\mu_1+\mu_2)}{\sqrt{\sigma_1^2 + \sigma_2^2}}\right)$$ and so \begin{align} F_Z(z) &= \Phi\left(\frac{z - (\mu_1+\mu_2)}{\sqrt{\sigma_1^2 + \sigma_2^2}}\right)\\ f_Z(z) &= \frac{1}{\sqrt{\sigma_1^2 + \sigma_2^2}}\phi\left(\frac{z - (\mu_1+\mu_2)}{\sqrt{\sigma_1^2 + \sigma_2^2}}\right) \end{align} where we can recognize this last expression as the density of a $N(\mu_1+\mu_2,\sigma_1^2+\sigma_2^2)$ random variable.