[proofplan] We prove multivariate normality by checking all one-dimensional linear projections of $Y=AX+b$. For an arbitrary vector $c \in \mathbb{R}^q$, the scalar projection $c^\top Y$ is rewritten as $(A^\top c)^\top X+c^\top b$, so the defining projection property of $X \sim \mathcal{N}_p(\mu,\Sigma)$ applies. We then compute the resulting mean and variance and identify them with the projection mean and variance associated to $A\mu+b$ and $A\Sigma A^\top$. [/proofplan] [step:Declare the affine image and reduce to scalar projections] Define \begin{align*} Y:(\Omega,\mathcal{F}) &\to (\mathbb{R}^q,\mathcal{B}(\mathbb{R}^q)) \\ \omega &\mapsto A X(\omega)+b. \end{align*} Since $X$ is measurable and the map $T:\mathbb{R}^p \to \mathbb{R}^q$ given by $T(x)=Ax+b$ is continuous, $Y=T\circ X$ is a random vector. To prove that $Y$ is $q$-variate normal, it suffices by the defining projection characterization of the multivariate normal distribution to show that, for every $c \in \mathbb{R}^q$, the real-valued random variable $c^\top Y$ is normally distributed with mean $c^\top(A\mu+b)$ and variance $c^\top A\Sigma A^\top c$. [guided] We first name the random vector whose distribution we want to identify: \begin{align*} Y:(\Omega,\mathcal{F}) &\to (\mathbb{R}^q,\mathcal{B}(\mathbb{R}^q)) \\ \omega &\mapsto A X(\omega)+b. \end{align*} The map $T:\mathbb{R}^p \to \mathbb{R}^q$ defined by $T(x)=Ax+b$ is continuous, hence Borel measurable. Because $X:(\Omega,\mathcal{F}) \to (\mathbb{R}^p,\mathcal{B}(\mathbb{R}^p))$ is measurable, the composition $Y=T\circ X$ is measurable. Thus $Y$ is indeed a $\mathbb{R}^q$-valued random vector. The definition of a $q$-variate normal random vector is tested through scalar projections: $Y$ is $\mathcal{N}_q(m,\Gamma)$ exactly when every linear projection $c^\top Y$, with $c \in \mathbb{R}^q$, is a univariate normal random variable with mean $c^\top m$ and variance $c^\top \Gamma c$. Therefore, to prove the theorem, we fix an arbitrary $c \in \mathbb{R}^q$ and compute the distribution of $c^\top Y$. [/guided] [/step] [step:Rewrite each projection as a projection of $X$ plus a constant] Fix $c \in \mathbb{R}^q$. Define $u \in \mathbb{R}^p$ by \begin{align*} u:=A^\top c. \end{align*} Then, for every $\omega \in \Omega$, \begin{align*} c^\top Y(\omega) &=c^\top(AX(\omega)+b) \\ &=c^\top A X(\omega)+c^\top b \\ &=(A^\top c)^\top X(\omega)+c^\top b \\ &=u^\top X(\omega)+c^\top b. \end{align*} Since $X \sim \mathcal{N}_p(\mu,\Sigma)$, the scalar random variable $u^\top X$ is normally distributed with mean $u^\top \mu$ and variance $u^\top \Sigma u$. Adding the deterministic constant $c^\top b$ preserves normality and shifts the mean by $c^\top b$, so $c^\top Y$ is normally distributed with mean $u^\top\mu+c^\top b$ and variance $u^\top\Sigma u$. [guided] Fix $c \in \mathbb{R}^q$. We want to express the scalar projection $c^\top Y$ in a form where the hypothesis on $X$ can be used. Define \begin{align*} u:=A^\top c \in \mathbb{R}^p. \end{align*} Then, for each $\omega \in \Omega$, \begin{align*} c^\top Y(\omega) &=c^\top(AX(\omega)+b) \\ &=c^\top A X(\omega)+c^\top b \\ &=(A^\top c)^\top X(\omega)+c^\top b \\ &=u^\top X(\omega)+c^\top b. \end{align*} This is the key reduction: $c^\top Y$ is not a new kind of random variable; it is a linear projection of $X$, followed by addition of the deterministic scalar $c^\top b$. Since $X \sim \mathcal{N}_p(\mu,\Sigma)$, the defining projection property gives \begin{align*} u^\top X \sim \mathcal{N}_1(u^\top\mu,u^\top\Sigma u). \end{align*} Adding the deterministic scalar $c^\top b$ shifts the mean and leaves the variance unchanged. Hence \begin{align*} c^\top Y \sim \mathcal{N}_1(u^\top\mu+c^\top b,u^\top\Sigma u). \end{align*} [/guided] [/step] [step:Identify the projected mean and variance] Using $u=A^\top c$, the mean computed above satisfies \begin{align*} u^\top\mu+c^\top b &=(A^\top c)^\top\mu+c^\top b \\ &=c^\top A\mu+c^\top b \\ &=c^\top(A\mu+b). \end{align*} Similarly, the variance satisfies \begin{align*} u^\top\Sigma u &=(A^\top c)^\top\Sigma(A^\top c) \\ &=c^\top A\Sigma A^\top c. \end{align*} Therefore, for every $c \in \mathbb{R}^q$, \begin{align*} c^\top Y \sim \mathcal{N}_1(c^\top(A\mu+b),c^\top A\Sigma A^\top c). \end{align*} By the projection characterization of the multivariate normal distribution, this proves \begin{align*} Y \sim \mathcal{N}_q(A\mu+b,A\Sigma A^\top). \end{align*} Since $Y=AX+b$, the desired conclusion follows. [/step]