贝叶斯多元线性回归

统计学中，贝叶斯多元线性回归是一种多元线性回归（预测结果为相关随机变量向量，而非标量随机变量的线性回归）的贝叶斯推断方法。这种方法的更一般论述见最小均方误差。

细节
考虑一回归问题，其中需要预测的自变量不是实标量而是相关实数组成的m维向量。与标准回归设置一样，有n个观测值，其中每个观测i包含k−1个解释变量，归为k维向量\mathbf{x}_i（添加值为1的虚拟变量，以允许截距系数）。对每个观测i，可以视作m个相关回归问题：
\begin{align}
y_{i,1} &= \mathbf{x}_i^\mathsf{T}\boldsymbol\beta_{1} + \epsilon_{i,1} \\
&\;\;\vdots \\
y_{i,m} &= \mathbf{x}_i^\mathsf{T}\boldsymbol\beta_{m} + \epsilon_{i,m}
\end{align}
其中误差集\{ \epsilon_{i,1}, \ldots, \epsilon_{i,m}\}都是相关的。等价地，也可以视作单一回归问题，其中结果是行向量\mathbf{y}_i^\mathsf{T}，回归系数向量排在一起：
\mathbf{y}_i^\mathsf{T} = \mathbf{x}_i^\mathsf{T}\mathbf{B} + \boldsymbol\epsilon_{i}^\mathsf{T}.

系数矩阵B是k \times m矩阵，其中每个回归问题的系数向量\boldsymbol\beta_1,\ldots,\boldsymbol\beta_m垂直排列在一起：
\mathbf{B} =
\begin{bmatrix}
\begin{pmatrix} \\ \boldsymbol\beta_1 \\ \\ \end{pmatrix}
\cdots
\begin{pmatrix} \\ \boldsymbol\beta_m \\ \\ \end{pmatrix}
\end{bmatrix}
=
\begin{bmatrix}
\begin{pmatrix}
\beta_{1,1} \\ \vdots \\ \beta_{k,1}
\end{pmatrix}
\cdots
\begin{pmatrix}
\beta_{1,m} \\ \vdots \\ \beta_{k,m}
\end{pmatrix}
\end{bmatrix}
.

每个观测i的噪声向量\boldsymbol\epsilon_{i}服从联合正态分布，因此给定观测的结果是相关的：
\boldsymbol\epsilon_i \sim N(0, \boldsymbol\Sigma_{\epsilon}).

可以将整个回归问题写成矩阵形式：
\mathbf{Y} = \mathbf{X}\mathbf{B} + \mathbf{E},
其中Y、E是n \times m矩阵。设计矩阵X是n \times k矩阵，观测如标准线性回归垂直排列：

\mathbf{X} = \begin{bmatrix} \mathbf{x}^\mathsf{T}_1 \\ \mathbf{x}^\mathsf{T}_2 \\ \vdots \\ \mathbf{x}^\mathsf{T}_n \end{bmatrix}
= \begin{bmatrix} x_{1,1} & \cdots & x_{1,k} \\
x_{2,1} & \cdots & x_{2,k} \\
\vdots & \ddots & \vdots \\
x_{n,1} & \cdots & x_{n,k}
\end{bmatrix}.

经典频率学派解利用摩尔－彭若斯广义逆，简单地估计回归系数矩阵\hat{\mathbf{B}}：
\hat{\mathbf{B}} = (\mathbf{X}^\mathsf{T}\mathbf{X})^{-1}\mathbf{X}^\mathsf{T}\mathbf{Y}.

求贝叶斯解，要先指定条件似然，再找到适当的共轭先验。与线性贝叶斯回归不同，可以指定一个自然的条件共轭先验（与规模相关）。

把条件似然写成
\rho(\mathbf{E}|\boldsymbol\Sigma_{\epsilon}) \propto |\boldsymbol\Sigma_{\epsilon}|^{-n/2} \exp\left(-\tfrac{1}{2} \operatorname{tr}\left(\mathbf{E}^\mathsf{T} \mathbf{E} \boldsymbol\Sigma_{\epsilon}^{-1}\right) \right) ,
误差\mathbf{E}表为\mathbf{Y},\mathbf{X},\mathbf{B}，则有
\rho(\mathbf{Y}|\mathbf{X},\mathbf{B},\boldsymbol\Sigma_{\epsilon}) \propto |\boldsymbol\Sigma_{\epsilon}|^{-n/2} \exp(-\tfrac{1}{2} \operatorname{tr}((\mathbf{Y}-\mathbf{X} \mathbf{B})^\mathsf{T} (\mathbf{Y}-\mathbf{X} \mathbf{B}) \boldsymbol\Sigma_{\epsilon}^{-1} ) ) ,

寻找一个自然共轭先验——联合密度\rho(\mathbf{B},\Sigma_{\epsilon})，其泛函形式与似然相同。由于似然在\mathbf{B}中是二次的，因此我们重写似然使其在(\mathbf{B}-\hat{\mathbf{B}})（与经典样本估计的差）是正态的。
用与相同的技术，可用矩阵形式的平方和分解指数项。不过此处还要用到矩阵微分（克罗内克积和向量化变换）。

首先，应用平方和得到新的似然表达式：
\rho(\mathbf{Y}|\mathbf{X},\mathbf{B},\boldsymbol\Sigma_{\epsilon}) \propto |\boldsymbol\Sigma_{\epsilon}|^{-(n-k)/2} \exp(-\operatorname{tr}(\tfrac{1}{2}\mathbf{S}^\mathsf{T} \mathbf{S} \boldsymbol\Sigma_{\epsilon}^{-1}))
|\boldsymbol\Sigma_{\epsilon}|^{-k/2} \exp(-\tfrac{1}{2} \operatorname{tr}((\mathbf{B}-\hat{\mathbf{B}})^\mathsf{T} \mathbf{X}^\mathsf{T} \mathbf{X}(\mathbf{B}-\hat{\mathbf{B}}) \boldsymbol\Sigma_{\epsilon}^{-1} ) )
,
\mathbf{S} = \mathbf{Y} - \mathbf{X}\hat{\mathbf{B}}

我们想开发一种先验的条件形式：
\rho(\mathbf{B},\boldsymbol\Sigma_{\epsilon}) = \rho(\boldsymbol\Sigma_{\epsilon})\rho(\mathbf{B}|\boldsymbol\Sigma_{\epsilon}),
其中\rho(\boldsymbol\Sigma_{\epsilon})服从逆威沙特分布，\rho(\mathbf{B}|\boldsymbol\Sigma_{\epsilon})是矩阵\mathbf{B}中某种形式的正态分布。这是通过向量化变换实现的，它将似然从矩阵\mathbf{B}, \hat{\mathbf{B}}的函数变换为向量\boldsymbol\beta = \operatorname{vec}(\mathbf{B}), \hat{\boldsymbol\beta} = \operatorname{vec}(\hat{\mathbf{B}})的函数。

\operatorname{tr}((\mathbf{B} - \hat{\mathbf{B}})^\mathsf{T}\mathbf{X}^\mathsf{T} \mathbf{X}(\mathbf{B} - \hat{\mathbf{B}}) \boldsymbol\Sigma_\epsilon^{-1}) = \operatorname{vec}(\mathbf{B} - \hat{\mathbf{B}})^\mathsf{T} \operatorname{vec}(\mathbf{X}^\mathsf{T} \mathbf{X}(\mathbf{B} - \hat{\mathbf{B}}) \boldsymbol\Sigma_{\epsilon}^{-1} )

令
\operatorname{vec}(\mathbf{X}^\mathsf{T} \mathbf{X}(\mathbf{B} - \hat{\mathbf{B}}) \boldsymbol\Sigma_{\epsilon}^{-1} ) = (\boldsymbol\Sigma_{\epsilon}^{-1} \otimes \mathbf{X}^\mathsf{T}\mathbf{X} )\operatorname{vec}(\mathbf{B} - \hat{\mathbf{B}}),
其中\mathbf{A} \otimes \mathbf{B}表示矩阵A、B的克罗内克积，其是外积的推广。

则
\begin{align}
&\operatorname{vec}(\mathbf{B} - \hat{\mathbf{B}})^\mathsf{T} (\boldsymbol\Sigma_{\epsilon}^{-1} \otimes \mathbf{X}^\mathsf{T}\mathbf{X} )\operatorname{vec}(\mathbf{B} - \hat{\mathbf{B}}) \\
&= (\boldsymbol\beta - \hat{\boldsymbol\beta})^\mathsf{T}(\boldsymbol\Sigma_{\epsilon}^{-1} \otimes \mathbf{X}^\mathsf{T}\mathbf{X} )(\boldsymbol\beta-\hat{\boldsymbol\beta})
\end{align}
产生的似然在(\boldsymbol\beta - \hat{\boldsymbol\beta})中正态。

有了更易理解的似然，就可以找到自然的（条件）共轭先验了。

共轭先验分布
由向量化的\boldsymbol\beta得到的自然共轭先验形式为
\rho(\boldsymbol\beta, \boldsymbol\Sigma_{\epsilon}) = \rho(\boldsymbol\Sigma_{\epsilon})\rho(\boldsymbol\beta|\boldsymbol\Sigma_{\epsilon}),
其中
\rho(\boldsymbol\Sigma_{\epsilon}) \sim \mathcal{W}^{-1}(\mathbf V_0,\boldsymbol\nu_0)

\rho(\boldsymbol\beta|\boldsymbol\Sigma_{\epsilon}) \sim N(\boldsymbol\beta_0, \boldsymbol\Sigma_{\epsilon} \otimes \boldsymbol\Lambda_0^{-1}).

后验分布
利用上述先验与似然，可得到后验
\begin{align}
\rho(\boldsymbol\beta,\boldsymbol\Sigma_{\epsilon}|\mathbf{Y},\mathbf{X})
\propto{}& |\boldsymbol\Sigma_{\epsilon}|^{-(\boldsymbol\nu_0 + m + 1)/2}\exp{(-\tfrac{1}{2}\operatorname{tr}(\mathbf V_0 \boldsymbol\Sigma_{\epsilon}^{-1}))} \\
&\times|\boldsymbol\Sigma_{\epsilon}|^{-k/2}\exp{(-\tfrac{1}{2} \operatorname{tr}((\mathbf{B}-\mathbf B_0)^\mathsf{T}\boldsymbol\Lambda_0(\mathbf{B}-\mathbf B_0)\boldsymbol\Sigma_{\epsilon}^{-1}))} \\
&\times|\boldsymbol\Sigma_{\epsilon}|^{-n/2}\exp{(-\tfrac{1}{2}\operatorname{tr}((\mathbf{Y}-\mathbf{XB})^\mathsf{T}(\mathbf{Y}-\mathbf{XB})\boldsymbol\Sigma_{\epsilon}^{-1}))},
\end{align}
其中\operatorname{vec}(\mathbf B_0) = \boldsymbol\beta_0。
记\boldsymbol\Lambda_0 = \mathbf{U}^\mathsf{T}\mathbf{U}，涉及\mathbf{B}的项可以分类为
\begin{align}
& \left(\mathbf{B} - \mathbf B_0\right)^\mathsf{T} \boldsymbol\Lambda_0 \left(\mathbf{B} - \mathbf B_0\right) + \left(\mathbf{Y} - \mathbf{XB}\right)^\mathsf{T} \left(\mathbf{Y} - \mathbf{XB}\right) \\
={}& \left(\begin{bmatrix}\mathbf Y \\ \mathbf U \mathbf B_0\end{bmatrix} - \begin{bmatrix}\mathbf{X}\\ \mathbf{U}\end{bmatrix}\mathbf{B}\right)^\mathsf{T} \left(\begin{bmatrix}\mathbf{Y}\\ \mathbf U \mathbf B_0\end{bmatrix}-\begin{bmatrix}\mathbf{X}\\ \mathbf{U}\end{bmatrix}\mathbf{B}\right) \\
={}& \left(\begin{bmatrix}\mathbf Y \\ \mathbf U \mathbf B_0\end{bmatrix} - \begin{bmatrix}\mathbf{X}\\ \mathbf{U}\end{bmatrix}\mathbf B_n\right)^\mathsf{T}\left(\begin{bmatrix}\mathbf{Y}\\ \mathbf U \mathbf B_0\end{bmatrix}-\begin{bmatrix}\mathbf{X}\\ \mathbf{U}\end{bmatrix}\mathbf B_n\right) + \left(\mathbf B - \mathbf B_n\right)^\mathsf{T} \left(\mathbf{X}^\mathsf{T} \mathbf{X} + \boldsymbol\Lambda_0\right) \left(\mathbf{B}-\mathbf B_n\right) \\
={}& \left(\mathbf{Y} - \mathbf X \mathbf B_n \right)^\mathsf{T} \left(\mathbf{Y} - \mathbf X \mathbf B_n\right) + \left(\mathbf B_0 - \mathbf B_n\right)^\mathsf{T} \boldsymbol\Lambda_0 \left(\mathbf B_0 - \mathbf B_n\right) + \left(\mathbf{B} - \mathbf B_n\right)^\mathsf{T} \left(\mathbf{X}^\mathsf{T} \mathbf{X} + \boldsymbol\Lambda_0\right)\left(\mathbf B - \mathbf B_n\right),
\end{align}
其中
\mathbf B_n = \left(\mathbf{X}^\mathsf{T}\mathbf{X} + \boldsymbol\Lambda_0\right)^{-1}\left(\mathbf{X}^\mathsf{T} \mathbf{X} \hat{\mathbf{B}} + \boldsymbol\Lambda_0\mathbf B_0\right) = \left(\mathbf{X}^\mathsf{T} \mathbf{X} + \boldsymbol\Lambda_0\right)^{-1}\left(\mathbf{X}^\mathsf{T} \mathbf{Y} + \boldsymbol\Lambda_0 \mathbf B_0\right).

现在可以用更有用的形式来写后验：
\begin{align}
\rho(\boldsymbol\beta,\boldsymbol\Sigma_{\epsilon}|\mathbf{Y},\mathbf{X})
\propto{}&|\boldsymbol\Sigma_{\epsilon}|^{-(\boldsymbol\nu_0 + m + n + 1)/2}\exp{(-\tfrac{1}{2}\operatorname{tr}((\mathbf V_0 + (\mathbf{Y}-\mathbf{XB_n})^\mathsf{T} (\mathbf{Y}-\mathbf{XB_n}) + (\mathbf B_n-\mathbf B_0)^\mathsf{T}\boldsymbol\Lambda_0(\mathbf B_n-\mathbf B_0))\boldsymbol\Sigma_{\epsilon}^{-1}))} \\
&\times|\boldsymbol\Sigma_{\epsilon}|^{-k/2}\exp{(-\tfrac{1}{2}\operatorname{tr}((\mathbf{B}-\mathbf B_n)^\mathsf{T} (\mathbf{X}^T\mathbf{X} + \boldsymbol\Lambda_0) (\mathbf{B}-\mathbf B_n)\boldsymbol\Sigma_{\epsilon}^{-1}))}.
\end{align}

其形式为逆威沙特分布乘以矩阵正态分布：
\rho(\boldsymbol\Sigma_{\epsilon}|\mathbf{Y},\mathbf{X}) \sim \mathcal{W}^{-1}(\mathbf V_n,\boldsymbol\nu_n)
\rho(\mathbf{B}|\mathbf{Y},\mathbf{X},\boldsymbol\Sigma_{\epsilon}) \sim \mathcal{MN}_{k,m}(\mathbf B_n, \boldsymbol\Lambda_n^{-1}, \boldsymbol\Sigma_{\epsilon}).

此后验的参数由下式给出
\mathbf V_n = \mathbf V_0 + (\mathbf{Y}-\mathbf{XB_n})^\mathsf{T}(\mathbf{Y}-\mathbf{XB_n}) + (\mathbf B_n - \mathbf B_0)^\mathsf{T}\boldsymbol\Lambda_0(\mathbf B_n-\mathbf B_0)
\boldsymbol\nu_n = \boldsymbol\nu_0 + n
\mathbf B_n = (\mathbf{X}^\mathsf{T}\mathbf{X} + \boldsymbol\Lambda_0)^{-1}(\mathbf{X}^\mathsf{T} \mathbf{Y} + \boldsymbol\Lambda_0\mathbf B_0)
\boldsymbol\Lambda_n = \mathbf{X}^\mathsf{T} \mathbf{X} + \boldsymbol\Lambda_0

另见
*

• 矩阵正态分布

参考文献
*
*
*