貝葉斯多元線性回歸

統計學中，貝葉斯多元線性回歸是一種多元線性回歸（預測結果為相關隨機變量向量，而非標量隨機變量的線性回歸）的貝葉斯推斷方法。這種方法的更一般論述見最小均方誤差。

## 細節

考慮一回歸問題，其中需要預測的自變量不是實標量而是相關實數組成的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$矩陣，觀測如標準線性回歸垂直排列：

$$<pre> \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}. </pre>$$

經典頻率學派線性最小二乘解利用摩爾－彭若斯廣義逆，簡單地估計回歸係數矩陣$\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$$