Linear Transformation of Gaussian Random Variable I've been trying to prove that if $\mathbf{x}$ is a random variable with multivariable normal distribution $Pr(\mathbf{x}) = Norm_\mathbf{x}[\mathbf{\mu}, \mathbf{\Sigma}]$ and $\mathbf{y}$ is a random variable where  $\mathbf{y} = \mathbf{Ax} + \mathbf{b}$, then $\mathbf{y}$ has distribution
$$Pr(\mathbf{y}) = Norm_\mathbf{y}[\mathbf{A\mu} + \mathbf{b}, \mathbf{A \Sigma A}^T]$$

Algebraically, I've proven that
$$Norm_\mathbf{x}[\mathbf{\mu}, \mathbf{\Sigma}] = \frac{\mid \mathbf{A \Sigma A}^T \mid}{\mid \mathbf{A} \mid}Norm_\mathbf{y}[\mathbf{A\mu} + \mathbf{b}, \mathbf{A \Sigma A}^T]$$
and if the right side is normalized with respect to $\mathbf{y}$, it results in
$$Norm_\mathbf{y}[\mathbf{A\mu} + \mathbf{b}, \mathbf{A \Sigma A}^T]$$
Intiuitively, I understand this as taking my normal distribution over $\mathbf{x}$, getting all of the relative frequencies of $\mathbf{y}$ from those, then renormalizing. Nonetheless, I would like to turn this into a more formal proof, rather than arguing "relative frequencies" - is there a way to formalize my proof, or is this as formal as it can get? I like this proof here, but the issue is that we have to assume $\mathbf{y}$ is be normally distributed. 
 A: I deprecate your dreadful notation, it will not help you achieve
any kind of rigor at all.
If $\mathbf X$ is a (jointly continuous) vector random variable with density function $f_{\mathbf X}({\mathbf x})$, then with $\mathbf G$ 
being an invertible matrix, the density function of 
$\mathbf{Y} =\mathbf{XG}$ is
$$f_{\mathbf Y}({\mathbf y}) 
= \frac{f_{\mathbf X}\left({\mathbf y\mathbf G^{-1}}\right)}{\vert\det{\mathbf G}\vert}.\tag{1}$$
(Note that I am using row vectors instead of the column vectors
that have been claimed to be more common usage in statistical circles).
This is a general result that holds regardless of whether multivariate
normality is assumed or not.  Also, the mean vector $E[\mathbf Y]$
of $\mathbf Y$ is related to the mean vector $E[\mathbf X]$ as 
$$E[\mathbf Y] = E[\mathbf{XG}] = E[\mathbf{X}]\mathbf{G}\tag{2}$$ while
the covariance matrix is
\begin{align}\Sigma_{\mathbf Y} 
&= E\left[\left(\mathbf Y - E[\mathbf Y]\right)^T\left(\mathbf Y - E[\mathbf Y]\right)\right]\\
&= E\left[\left(\mathbf {XG} - E[\mathbf X]\mathbf G\right)^T\left(\mathbf {XG} - E[\mathbf X]\mathbf G\right)\right]\\
&= \mathbf G^T E\left[\left(\mathbf {X} - E[\mathbf X]\right)^T
\left(\mathbf {X} - E[\mathbf X]\right)\right]\mathbf G\\
&= \mathbf G^T\Sigma_{\mathbf X}\mathbf G. \tag{3}
\end{align}
For the special case of $\mathbf X$ having a $n$-variate 
normal distribution and the covariance matrix
$\Sigma_{\mathbf X}$ being a strictly positive invertible
$n\times n$ matrix, the density function
$f_{\mathbf X}({\mathbf x})$ is of the form
$$f_{\mathbf X}({\mathbf x})
= \frac{1}{(2\pi)^{n/2}\sqrt{\det\left(\Sigma_{\mathbf X}\right)}}
\exp\left[-\frac 12 \left(\mathbf {x} - E[\mathbf X]\right)
\Sigma_{\mathbf X}^{-1}\left(\mathbf {x} - E[\mathbf X]\right)^T\right]
\tag{4}$$
and so from $(1)$ we get
\begin{align}
f_{\mathbf Y}({\mathbf y}) 
&= \frac{f_{\mathbf X}\left({\mathbf y\mathbf G^{-1}}\right)}{\vert\det{\mathbf G}\vert}\\
&= \frac{1}{(2\pi)^{n/2}\sqrt{\det\left(\Sigma_{\mathbf X}\right)}
\vert\det{\mathbf G}\vert}
\exp\left[-\frac 12 \left(\mathbf{yG}^{-1} - E[\mathbf X]\right)
\Sigma_{\mathbf X}^{-1}\left(\mathbf{yG}^{-1} - E[\mathbf X]\right)^T\right]\\
&= \frac{1}{(2\pi)^{n/2}\sqrt{\det\left(\Sigma_{\mathbf X}\right)}
\vert\det{\mathbf G}\vert^{1/2}\vert\det{\mathbf G}^T\vert^{1/2}}
\exp\left[-\frac 12 \left(\mathbf{yG}^{-1} - E[\mathbf X]\right)
\Sigma_{\mathbf X}^{-1}\left(\mathbf{yG}^{-1} - E[\mathbf X]\right)^T\right]\\
&= \frac{1}{(2\pi)^{n/2}\sqrt{\left|\det\left(\mathbf G^T\Sigma_{\mathbf X}\mathbf G\right)\right|}}
\exp\left[-\frac 12 \left(\mathbf{y} - E[\mathbf X]\mathbf G\right)
\mathbf G^{-1}
\Sigma_{\mathbf X}^{-1}\left(\left(\mathbf{y} 
- E[\mathbf X]\mathbf G\right)\mathbf G^{-1}\right)^T\right]\\
&= \frac{1}{(2\pi)^{n/2}\sqrt{\left|\det\left(\mathbf G^T\Sigma_{\mathbf X}\mathbf G\right)\right|}}
\exp\left[-\frac 12 \left(\mathbf{y} - E[\mathbf Y]\right)
\mathbf G^{-1}
\Sigma_{\mathbf X}^{-1}\left(\mathbf G^{-1}\right)^T\left(\mathbf{y} 
- E[\mathbf Y]\right)^T\right]\\
&= \frac{1}{(2\pi)^{n/2}\sqrt{\left|\det\Sigma_{\mathbf Y}\right|}}
\exp\left[-\frac 12 \left(\mathbf{y} - E[\mathbf Y]\right)
\Sigma_{\mathbf Y}^{-1}\left(\mathbf{y} 
- E[\mathbf Y]\right)^T\right]
\end{align}
which shows that $\mathbf Y = \mathbf{XG}$ also has a $n$-variate normal density with mean vector and covariance matrix as given in (2) and (3).
A: This can be shown very succinctly by using the characteristic function of distributions. Let $\phi_X(t) = E[ \exp(i t ^ \mathsf T X) ]$ be the characteristic function of a random variable $X \in \mathbb R^n$.
If $x$ is normally distributed $x \sim \mathcal N(\mu, \Sigma)$, then we have $\phi_x(t) = \exp \Big (i t ^ \mathsf T \mu - \tfrac {1}{2} t ^ \mathsf T \Sigma t \Big )$.
If $y = A x + b$, then
\begin{align}
\phi_y(t) &= E[ \exp \Big(i t ^ \mathsf T (A x + b) \Big) ] \\
&= E[ \exp \Big(i t ^ \mathsf T b \Big) \exp \Big(i t ^ \mathsf T A x \Big) ] \\
&= \exp \Big(i t ^ \mathsf T b \Big) E[ \exp \Big(i (A ^ \mathsf T t) ^ \mathsf T x \Big) ] \\
&= \exp \Big(i t ^ \mathsf T b \Big) \phi_x (A ^ \mathsf T t) \\
&= \exp \Big(i t ^ \mathsf T b \Big) \exp \Big (i (A ^ \mathsf T t) ^ \mathsf T \mu - \tfrac {1}{2} (A ^ \mathsf T t) ^ \mathsf T \Sigma (A ^ \mathsf T t) \Big ) \\
&= \exp \Big (i t ^ \mathsf T (A \mu + b) - \tfrac {1}{2} t ^ \mathsf T A \Sigma A ^ \mathsf T t \Big ) \\
\end{align}
Since the characteristic function uniquely defines the distribution, we have $y \sim \mathcal N(A \mu + b, A \Sigma A ^ \mathsf T)$ as wanted.
