最小範數二次無偏估計法

在統計學中，最小範數二次不偏估計法（minimum norm quadratic unbiased estimation，MINQUE）理論 由卡利安普迪·拉達克里希納·拉奧提出。最小範數二次不偏估計法與矩估計法或最大似然估計等其他估計方法一起構成估計理論。與最佳線性不偏估計理論類似，最小範數二次不偏估計法是專門針對線性回歸的模型。該方法最初是為了估計多元線性回歸中的異方差誤差變異數而設計的。對於混合效應模型中的變異數分量，最小範數二次不偏估計法估計量也為最大似然估計量或限制極大似然估計量提供了一種替代方案。最小範數二次不偏估計法估計量是響應變量的二次型，用於估計變異數的線性函數。

## 原理

在混合效應模型中，隨機向量 $\mathbf{Y} \in \mathbb{R}^n$具有以下線性結構：

$$\mathbf{Y} = \mathbf{X}\boldsymbol\beta + \mathbf{U}_1 \boldsymbol\xi_1 + \cdots + \mathbf{U}_k \boldsymbol\xi_k$$

此處， $\mathbf{X} \in \mathbb{R}^{n\times m}$是固定效應的設計矩陣， $\boldsymbol\beta \in \mathbb{R}^m$代表未知的固定效應參數， $\mathbf{U}_i \in \mathbb{R}^{n\times c_i}$是設計矩陣第$i$個隨機效應分量，而$\boldsymbol\xi_i\in\mathbb{R}^{c_i}$是第$i$ 個隨機效應分量的隨機向量。隨機效應的均值被假設為零（$\mathbb{E}[\boldsymbol\xi_i]=\mathbf{0}$）且獨立（$\mathbb{V}[\boldsymbol\xi_i]=\sigma^2_i\mathbf{I}_{c_i}$），而任意兩個隨機效應向量也是獨立的（$\mathbb{V}[\boldsymbol\xi_i, \boldsymbol\xi_j]=\mathbf{0}\, \forall i\neq j$）。其中未知變異數$\sigma^2_1,\cdots,\sigma^2_k$代表模型的變異數分量。

這是一個涵蓋常用線性回歸模型的通用模型。

- 高斯-馬爾可夫模型：如果我們考慮一個單分量模型，其中$\mathbf{U}_1=\mathbf{I}_n$。那麼，該模型就等價於高斯-馬爾可夫模型中 $\mathbf{Y}=\mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon$和$\mathbb{E}[\boldsymbol\epsilon]=\mathbf{0}$和$\mathbb{V}[\boldsymbol\epsilon]=\sigma^2_1 \mathbf{I}_n$。
- 異方差模型：每組隨機變量$\mathbf{Y}$具有共同變異數的變量，並可以建模為具有適當係數的個體變異數分量： $\mathbf{U}_i$。

該模型的緊湊表示如下，其中$\mathbf{U} = \left[\begin{array}{c|c|c}\mathbf{U}_1&\cdots&\mathbf{U}_k\end{array}\right]$和$\boldsymbol\xi^\top = \left[\begin{array}{c|c|c} \boldsymbol\xi_1^\top&\cdots&\boldsymbol\xi_k^\top\end{array}\right]$：

$$\mathbf{Y}=\mathbf{X}\boldsymbol\beta+\mathbf{U}\boldsymbol\xi$$

以下模型沒有對除了第一動差和第二動差之外的$\mathbf{Y}$分布做任何假設：

$$\mathbb{E}[\mathbf{Y}] = \mathbf{X}\boldsymbol\beta$$

$$\mathbb{V}[\mathbf{Y}]=\sigma^2_1\mathbf{U}_1\mathbf{U}_1^\top + \cdots + \sigma^2_k \mathbf{U}_k \mathbf{U}_k^\top \equiv \sigma^2_1\mathbf{V}_1 + \cdots + \sigma^2_k \mathbf{V}_k$$

最小範數二次不偏估計法的目標透過是使用二次型$\hat{\theta}=\mathbf{Y}^\top \mathbf{A} \mathbf{Y}$估算$\theta = \sum_{i=1}^k p_i \sigma^2_i$。最小範數二次不偏估計法估計量是通過找出一個使得估計量具有一些理想的性質的矩陣$\mathbf{A}$而得到的。下面描述矩陣$\mathbf{A}$尋找方式：

### 約束最小範數二次不偏估計法的最佳估計量性質

#### 固定效應的平移不變性

考慮對原始固定效應的轉換，產生一個新的固定效應參數$\boldsymbol\gamma=\boldsymbol\beta - \boldsymbol\beta_0$，則新的等效模型如下所示：

$$\mathbf{Y} - \mathbf{X}\boldsymbol\beta_0 = \mathbf{X}\boldsymbol\gamma + \mathbf{U}\boldsymbol\xi$$

在這個等效模型下，最小範數二次不偏估計法估計量現在是$(\mathbf{Y} - \mathbf{X}\boldsymbol\beta_0)^\top \mathbf{A} (\mathbf{Y} - \mathbf{X}\boldsymbol\beta_0)$。

卡利安普迪·拉達克里希納·拉奧認為，由於底層模型是等價的，因此該估計量應該等於$\mathbf{Y}^\top \mathbf{A} \mathbf{Y}$ 此估計量可以通過尋找符合條件的$\mathbf{A}$（能使得$\mathbf{A}\mathbf{X} = \mathbf{0}$以確保二次型的展開式中除$\mathbf{Y}^\top \mathbf{A} \mathbf{Y}$係數之外的所有係數均為零）來尋得。

#### 不偏估計

假設如上文所述限制$\mathbf{A}\mathbf{X} = \mathbf{0}$，最小範數二次不偏估計法估計量具有以下形式：

$$\begin{align} \hat{\theta} &= \mathbf{Y}^\top \mathbf{A} \mathbf{Y}\\ &= (\mathbf{X}\boldsymbol\beta + \mathbf{U}\boldsymbol\xi)^\top \mathbf{A} (\mathbf{X}\boldsymbol\beta + \mathbf{U}\boldsymbol\xi)\\ &= \boldsymbol\xi^\top\mathbf{U}^\top\mathbf{A}\mathbf{U}\boldsymbol\xi \end{align}$$

為了確保該估計量是不偏的，需要確保估計量期望值$\mathbb{E}[\hat{\theta}]$必須等於所關注的參數$\theta$，在下方公式中，由於各分量彼此獨立，估計量的期望值可以分解為每個分量的期望值之和。此外，跡的循環性質被用於評估與 $\boldsymbol\xi_i$有關的期望值。

$$\begin{align} \mathbb{E}[\hat{\theta}] &= \mathbb{E}[\boldsymbol\xi^\top \mathbf{U}^\top \mathbf{A} \mathbf{U} \boldsymbol\xi]\\ &= \sum_{i=1}^k \mathbb{E}[\boldsymbol\xi_i^\top\mathbf{U}_i^\top\mathbf{A}\mathbf{U}_i\boldsymbol\xi_i]\\ &= \sum_{i=1}^k \sigma_i^2 \mathrm{Tr}[\mathbf{U}_i^\top \mathbf{A} \mathbf{U}_i] \end{align}$$

為了確保該估計量不偏，卡利安普迪·拉達克里希納·拉奧建議設定$\sum_{i=1}^k \sigma_i^2 \mathrm{Tr}[\mathbf{U}_i^\top \mathbf{A} \mathbf{U}_i] = \sum_{i=1}^k p_i \sigma_i^2$，以藉約束$\mathbf{A}$以對所有組件使得$\mathrm{Tr}[\mathbf{U}_i^\top \mathbf{A} \mathbf{U}_i] = \mathrm{Tr}[\mathbf{A}\mathbf{V}_i] = p_i$。

#### 最小範數

拉奧認為，如果$\boldsymbol\xi$被觀測到，因為$\mathbb{E}[\boldsymbol\xi_i^\top\boldsymbol\xi_i]=c_i \sigma_i^2$，一個對$\theta$的「自然」估計量將如以下公式所示。其中， $\boldsymbol\Delta$被定義為對角矩陣：

$$\frac{p_1}{c_1}\boldsymbol\xi_1^\top\boldsymbol\xi_1 + \cdots + \frac{p_k}{c_k}\boldsymbol\xi_k^\top\boldsymbol\xi_k = \boldsymbol\xi^\top\left[\mathrm{diag}\left(\frac{p_1}{c_i},\cdots,\frac{p_k}{c_k}\right)\right]\boldsymbol\xi \equiv \boldsymbol\xi^\top\boldsymbol\Delta\boldsymbol\xi$$

所提出的估計量與自然估計量之間的差是$\boldsymbol\xi^\top (\mathbf{U}^\top \mathbf{A} \mathbf{U} - \boldsymbol\Delta)\boldsymbol\xi$。此差可以通過最小化矩陣範數$\lVert \mathbf{U}^\top\mathbf{A}\mathbf{U}-\boldsymbol\Delta \rVert$被最小化。

### 運算流程

鑑於上述最佳性質推導出的約束條件和優化策略，最小範數二次不偏估計法估計量$\hat{\theta}$為了$\theta=\sum_{i=1}^k p_i\sigma_i^2$是通過選擇一個在符合$\mathbf{A}\mathbf{X}=\mathbf{0}$且$\mathrm{Tr}[\mathbf{A}\mathbf{V}_i]=p_i$條件下能最小化$\lVert \mathbf{U}^\top\mathbf{A}\mathbf{U}-\boldsymbol\Delta \rVert$矩陣$\mathbf{A}$推導出來的。

## 估計量範例

### 同方差性誤差的標準估計量

在高斯-馬爾可夫模型中, 變異數 $\sigma^2$ 由以下公式推導：

$$s^2 = \frac{1}{n-m}(\mathbf{Y}-\mathbf{X}\hat{\boldsymbol\beta})^\top(\mathbf{Y}-\mathbf{X}\hat{\boldsymbol\beta})$$

該估計量是不偏的，並且可以證明它能最小化如下形式的歐幾里得範數： $\lVert \mathbf{U}^\top\mathbf{A}\mathbf{U}-\boldsymbol\Delta \rVert$ 因此，高斯-馬爾可夫模型中誤差變異數的標準估計量是最小範數二次不偏估計法估計量。

### 具有共同均值和異方差性誤差的隨機變量

對於隨機變量$Y_1,\cdots,Y_n$均值相同但變異數不同的$\sigma^2_1,\cdots,\sigma^2_n$，最小範數二次不偏估計法估計量$\sigma^2_i$是$\frac{n}{n-2}(Y_i - \overline{Y})^2 - \frac{s^2}{n - 2}$，其中$\overline{Y} = \frac{1}{n} \sum_{i=1}^n Y_i$且$s^2 = \frac{1}{n-1} \sum_{i=1}^n (Y_i - \overline{Y})^2$。

### 變異數分量估計量

卡利安普迪·拉達克里希納·拉奧提出了一種基於最小化歐幾里得範數的變異數分量模型之最小範數二次不偏估計法估計量。歐氏範數$\lVert \cdot \rVert_2$是矩陣中所有元素平方和的平方根。在評估以下範數時， $\mathbf{V}=\mathbf{V}_1+\cdots+\mathbf{V}_k = \mathbf{U} \mathbf{U}^\top$此外，利用跡的循環性質，可得： $\mathrm{Tr}[\mathbf{U}^\top\mathbf{A}\mathbf{U}\boldsymbol\Delta] = \mathrm{Tr}[\mathbf{A}\mathbf{U}\boldsymbol\Delta\mathbf{U}^\top] = \mathrm{Tr}\left[\sum_{i=1}^k \frac{p_i}{c_i} \mathbf{A}\mathbf{V}_i \right] = \mathrm{Tr}[\boldsymbol\Delta\boldsymbol\Delta]$。

$$\begin{align} \lVert \mathbf{U}^\top\mathbf{A}\mathbf{U} - \boldsymbol\Delta \rVert^2_2 &= (\mathbf{U}^\top\mathbf{A}\mathbf{U} - \boldsymbol\Delta)^\top (\mathbf{U}^\top\mathbf{A}\mathbf{U} - \boldsymbol\Delta)\\ &= \mathrm{Tr}[\mathbf{U}^\top\mathbf{A}\mathbf{U}\mathbf{U}\mathbf{A}\mathbf{U}^\top] - \mathrm{Tr}[2\mathbf{U}^\top\mathbf{A}\mathbf{U}\boldsymbol\Delta] + \mathrm{Tr}[\boldsymbol\Delta\boldsymbol\Delta]\\ &= \mathrm{Tr}[\mathbf{A}\mathbf{V}\mathbf{A}\mathbf{V}] - \mathrm{Tr}[\boldsymbol\Delta\boldsymbol\Delta] \end{align}$$

由於$\mathrm{Tr}[\boldsymbol\Delta\boldsymbol\Delta]$不取決於$\mathbf{A}$，具有歐幾里得範數的最小範數二次不偏估計法矩陣可通過能最大限度地減少$\mathrm{Tr}[\mathbf{A}\mathbf{V}\mathbf{A}\mathbf{V}]$且受前面所提最小範數二次不偏估計法限制的識別矩陣$\mathbf{A}$取得。

拉奧指出滿足此最佳化問題的矩陣$\mathbf{A}$是：

$\mathbf{A}_\star=\sum_{i=1}^k \lambda_i \mathbf{R}\mathbf{V}_i\mathbf{R}$,

其中$\mathbf{R} = \mathbf{V}^{-1}(\mathbf{I}-\mathbf{P})$， $\mathbf{P}=\mathbf{X}(\mathbf{X}^\top\mathbf{V}^{-1}\mathbf{X})^{-}\mathbf{X}^\top\mathbf{V}^{-1}$兩者是投影到$\mathbf{X}$列空間的矩陣，而$(\cdot)^{-}$表示一矩陣的廣義逆陣。

因此，最小範數二次不偏估計法估計量如下所示，其中向量 $\boldsymbol\lambda$和$\mathbf{Q}$是根據總和定義的：

$$\begin{align} \hat{\theta} &= \mathbf{Y}^\top \mathbf{A}_\star\mathbf{Y}\\ &= \sum_{i=1}^k \lambda_i \mathbf{Y}^\top\mathbf{R}\mathbf{V}_i\mathbf{R}\mathbf{Y}\\ &\equiv\sum_{i=1}^k \lambda_i Q_i\\ &\equiv \boldsymbol\lambda^\top \mathbf{Q} \end{align}$$

向量$\boldsymbol\lambda$是通過使用約束條件$\mathrm{Tr}[\mathbf{A}_\star\mathbf{V}_i]=p_i$獲得的。也就是說，該向量$\boldsymbol\lambda$表示以下方程組的解。$\forall j\in\{1,\cdots,k\}$：

$$\begin{align} \mathrm{Tr}[\mathbf{A}_\star\mathbf{V}_j] &= p_j\\ \mathrm{Tr}\left[ \sum_{i=1}^k \lambda_i \mathbf{R}\mathbf{V}_i\mathbf{R}\mathbf{V}_j \right] &= p_j\\ \sum_{i=1}^k \lambda_i \mathrm{Tr}[\mathbf{R}\mathbf{V}_i\mathbf{R}\mathbf{V}_j] &= p_j \end{align}$$

這可以寫成矩陣乘積的形式 $\mathbf{S}\boldsymbol\lambda=\mathbf{p}$，其中$\mathbf{p}=[p_1\,\cdots\,p_k]^\top$而$\mathbf{S}$如下所示：

$$\mathbf{S}=\begin{bmatrix} \mathrm{Tr}[\mathbf{R}\mathbf{V}_1\mathbf{R}\mathbf{V}_1] & \cdots & \mathrm{Tr}[\mathbf{R}\mathbf{V}_k\mathbf{R}\mathbf{V}_1]\\ \vdots & \ddots & \vdots\\ \mathrm{Tr}[\mathbf{R}\mathbf{V}_1\mathbf{R}\mathbf{V}_k] & \cdots & \mathrm{Tr}[\mathbf{R}\mathbf{V}_k\mathbf{R}\mathbf{V}_k] \end{bmatrix}$$

因此， $\boldsymbol\lambda=\mathbf{S}^{-}\mathbf{p}$，代表 最小範數二次不偏估計法中$\hat{\theta}=\boldsymbol\lambda^\top\mathbf{Q}=\mathbf{p}^\top(\mathbf{S}^{-})^\top\mathbf{Q}=\mathbf{p}^\top\mathbf{S}^{-}\mathbf{Q}$。此外，$\theta=\sum_{i=1}^k p_i \sigma_i^2 = \mathbf{p}^\top\boldsymbol\sigma$，其中$\boldsymbol\sigma = [\sigma^2_1\,\cdots\,\sigma^2_k]^\top$。因此，變異數分量的估計量為$\hat{\boldsymbol\sigma}=\mathbf{S}^{-}\mathbf{Q}$。

## 擴展

最小範數二次不偏估計法估計量可以在不滿足不變性準則的情況下獲得，在這種情況下，估計量只是不偏的，並且最小化範數。這類估計量對最小化問題有略微不同的限制。

此模型可以擴展到估計協變異數分量。在這種模型中，分量的隨機效應被假設具有共同的協方差結構$\mathbb{V}[\boldsymbol\xi_i]=\boldsymbol\Sigma$。此外，一種用於混合方差和協方差分量的最小範數二次不偏估計法估計量之模型也被提出。在該模型中，在$i\in \{1,\cdots,s\}$下 $\mathbb{V}[\boldsymbol\xi_i]=\boldsymbol\Sigma$，而$i\in\{s+1,\cdots,k\}$下 $\mathbb{V}[\boldsymbol\xi_i]= \sigma_i^2\mathbf{I}_{c_i}$。