はじめに

正規分布の再生性は独立な場合が有名ですが、独立でない場合も実は正規分布に従い、再生性を持ちます。
独立な場合の証明は以下が丁寧にまとまっています。

結論

独立でない場合の結論は以下です。
相関がある場合は、分散に相関係数の項が追加されます。

$a,b$を定数とする。確率変数$X,Y$について
$ X\sim N(μ_{x}, σ_{x}^2)\quad Y\sim N(μ_{y}, σ_{y}^2)\quad$2分布の相関係数を$ρ$とする。
このとき、$aX+bY\sim N(\hspace{1pt}aμ_x+bμ_y,\hspace{4pt} a^2σ_x^2+b^2σ_y^2+2abρσ_xσ_y)$である

和$X+Y$の分布を考えたいときは$a=1, b=1$として

$X+Y\sim N(\hspace{1pt}μ_x+μ_y,\hspace{4pt} σ_x^2+σ_y^2+2ρσ_xσ_y)$

差$X-Y$の分布を考えたいときは$a=1, b=-1$として

$X-Y\sim N(\hspace{1pt}μ_x-μ_y,\hspace{4pt} σ_x^2+σ_y^2-2ρσ_xσ_y)$

のように分かります！

上記の結果を以下の2種類で証明します。
1. モーメント母関数を使った証明
2. 確率密度関数から導出する証明(現在制作中)

上で導入した文字を証明中そのまま使います。

1. モーメント母関数を使った証明

モーメント母関数と確率分布は1対1で対応します。
この証明は以下にあります。

https://academ-aid.com/statistics/unique
https://su-butsu-kikaigakusyuu.hatenablog.com/entry/2018/08/01/004535 (測度論から)

そのため以下の方針で証明します！

方針

確率変数$aX+bY$のモーメント母関数が$N(\hspace{1pt}aμ_x+bμ_y,\hspace{4pt} a^2σ_x^2+b^2σ_y^2+2abρσ_xσ_y)$に従う確率変数のモーメント母関数に一致することを示す

解答

確率変数$aX+bY$のモーメント母関数$M_{aX+bY}(t)$と書く。

  \begin{align}
    M_{aX+bY}(t) &= \text{E}[\text{exp}(t(aX+bY))] \\
                 &= \text{E}[\text{exp}(taX+tbY)] \\
                 &= \int_{-∞}^{∞}\int_{-∞}^{∞}\text{exp}(tax+tby)f(x,y)dxdy \\
                 &= \int_{-∞}^{∞}\int_{-∞}^{∞}\text{exp}(tax+tby)・\frac{1}{2 \pi \sigma_x \sigma_y \sqrt{1 - \rho^2}}\\
                 &  \exp\left( -\frac{1}{2(1 - \rho^2)} \left[ \frac{(x - \mu_x)^2}{\sigma_x^2} \\ - 2\rho \frac{(x - \mu_x)(y - \mu_y)}{\sigma_x \sigma_y} + \frac{(y - \mu_y)^2}{\sigma_y^2} \right] \right)  
  \end{align}

ここで重積分において$\displaystyle s=\frac{x-\mu_{x}}{\sigma_{x}}, \quad t=\displaystyle \frac{y-\mu_{y}}{\sigma_{y}}$と変数変換を行う。
$x=μ_x+sσ_x, y=μ_y+sσ_y$であることと、
|J| = 
\begin{vmatrix}
\frac{\partial x}{\partial s} & \frac{\partial x}{\partial t} \\
\frac{\partial y}{\partial s} & \frac{\partial y}{\partial t}
\end{vmatrix}
   = 
\begin{vmatrix}
σ_x & 0 \\
0 & σ_y
\end{vmatrix}
   = σ_xσ_y
より$dxdy = |J|dsdt = σ_xσ_ydsdt$をふまえて計算を進めると、

  \begin{align}
    (与式) &= \int_{-∞}^{∞}\int_{-∞}^{∞}\frac{1}{2\pi\sqrt{1-ρ^2}}\text{exp}\Bigg(ta(μ_x+σ_xs)+tb(μ_y+σ_yt) \\
           & -\frac{1}{2(1-ρ^2)}(s^2-2ρst+t^2)\Bigg)dsdt \\
           &= \text{exp}(taμ_x+tbμ_y) \int_{-∞}^{∞}\int_{-∞}^{∞}\frac{1}{2\pi\sqrt{1-ρ^2}}\text{exp}\Bigg(taσ_xs+tbσ_yt \\
           & -\frac{1}{2(1-ρ^2)}(s^2-2ρst+t^2)) \Bigg) dsdt \\
           % この行で行列表記に書き直す
           &= \text{exp}(taμ_x+tbμ_y) \int_{-∞}^{∞}\int_{-∞}^{∞}\frac{1}{2\pi\sqrt{1-ρ^2}} \\
           & \text{exp}\Bigg( -\frac{1}{2}
\begin{pmatrix}
 s & t
\end{pmatrix}
\begin{pmatrix}
 \frac{1}{1-ρ^2} & -\frac{ρ}{1-ρ^2} \\
 -\frac{ρ}{1-ρ^2} & \frac{1}{1-ρ^2} \\
\end{pmatrix}
\begin{pmatrix}
 s \\ t 
\end{pmatrix}
+
\begin{pmatrix}
 taμ_x & taμ_y 
\end{pmatrix}
\begin{pmatrix}
 s \\ t 
\end{pmatrix}
\Bigg) dsdt
  \end{align}

\mathbb{z}=
\begin{pmatrix}
s \\ t 
\end{pmatrix},\,\,
A=
\begin{pmatrix}
\frac{1}{1-ρ^2} & -\frac{ρ}{1-ρ^2} \\
-\frac{ρ}{1-ρ^2} & \frac{1}{1-ρ^2} \\
\end{pmatrix},\,\,
\mathbb{b}=
\begin{pmatrix}
taμ_x & taμ_y 
\end{pmatrix}とすると、
(被積分関数の指数部)=-\frac{1}{2}\mathbb{z}^\mathsf{T}A\mathbb{z}+\mathbb{b}^\mathsf{T}\mathbb{z}
ここで2次形式の平方完成を用いる。
-\frac{1}{2}\mathbb{z}^\mathsf{T}A\mathbb{z}+\mathbb{b}^\mathsf{T}\mathbb{z}=-\frac{1}{2} (\mathbb{z} - A^{-1} \mathbb{b})^\mathsf{T} A (\mathbb{z} - A^{-1} \mathbb{b}) + \frac{1}{2} \mathbb{b}^\mathsf{T} A^{-1} \mathbb{b}
この証明は係数比較によりできます。少し形が違いますが、以下のサイトで2次形式の平方完成の概要と証明を行っています。
https://blog.monophile.net/posts/20131227_quadratic_square.html
これを用いてさらに変形を進めます。

  \begin{align}
    (与式) &= \text{exp}(taμ_x+tbμ_y)\int_{-∞}^{∞}\int_{-∞}^{∞}\frac{1}{2\pi\sqrt{1-ρ^2}}\text{exp}\Bigg(-\frac{1}{2} (\mathbb{z} - A^{-1} \mathbb{b})^\mathsf{T} A  \\
           & (\mathbb{z} - A^{-1} \mathbb{b}) + \frac{1}{2} \mathbb{b}^\mathsf{T} A^{-1} \mathbb{b}\Bigg)dsdt \\
           &= \text{exp}(taμ_x+tbμ_y+\frac{1}{2} \mathbb{b}^\mathsf{T} A^{-1} \mathbb{b})\int_{-∞}^{∞}\int_{-∞}^{∞}\frac{1}{2\pi\sqrt{1-ρ^2}} \\
           & \text{exp}\Bigg(-\frac{1}{2} (\mathbb{z} - A^{-1} \mathbb{b})^\mathsf{T} A (\mathbb{z} - A^{-1} \mathbb{b})\Bigg)dsdt \\
           &= \text{exp}(t(aμ_x+bμ_y)+\frac{t^2}{2}(a^2μ_x^2+2abμ_xμ_y+b^2σ_y^2))
  \end{align}

$(\because$ 被積分関数は平均ベクトル$A^{-1} \mathbb{b}$, 分散$A$の2次元正規分布で、全面積なので1である$)$

これは正規分布$N(\hspace{1pt}aμ_x+bμ_y,\hspace{4pt} a^2σ_x^2+b^2σ_y^2+2abρσ_xσ_y)$のモーメント母関数に一致しています。

正規分布のモーメント母関数は以下に丁寧にまとめられています。
https://mathlandscape.com/normal-distrib-char/

以上により$aX+bY\sim N(\hspace{1pt}aμ_x+bμ_y,\hspace{4pt} a^2σ_x^2+b^2σ_y^2+2abρσ_xσ_y)$であることが示されました。

2. 確率密度関数から導出する証明

証明できることは確認しましたが、モーメント母関数による解答より骨が折れる証明です。
現在執筆中です。少々お待ちくださいませ。