在阅读 The Elements of Statistical Learning 第三章的时候，有一个式子我没有弄明白：设 $Y = X\beta + \epsilon$，其中 $\epsilon$ 的均值为 0，方差为 $\sigma^2$；再设 $X$ 是 $N \times (p+1)$ 的矩阵（每条训练样本含常数项 1），那么对 $\sigma^2$ 的无偏估计是 $$\hat{\sigma}^2 = \frac{1}{N-p-1}\sum_{i=1}^N(y_i-\hat{y}_i)^2$$ 其中 $\hat{Y} = X\hat{\beta}$，$\hat{\beta}$ 是用 least square 得到的参数。
这个式子最奇怪的就是前面的系数 $1/(N-p-1)$。我们一般计算方差的时候都是以 $1/N$ 作为系数的，样本方差是以 $1/(N-1)$ 作为系数的。那么样本方差的 $1/(N-1)$ 和上面那个式子的 $1/(N-p-1)$ 是怎么来的呢？
我们先来看一下在一般情况下计算方差时，为什么以 $1/N$ 为系数的估计是有偏的。设 $x_1, x_2, \dots, x_n$ 是我们观察到的数据，它们的均值 $\bar{x} = (x_1 + x_2 + \dots + x_n)/n$，我们估计的方差是 $$\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2$$ 设数据真实的期望为 $\mu$，真实的方差为 $\sigma^2$，则我们估计的方差的期望为 $$\mathbb{E}(\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2) \\ = \mathbb{E}(\frac{1}{n}\sum_{i=1}^n((x_i-\mu)-(\bar{x}-\mu))^2) \\ = \mathbb{E}(\frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2 - \frac{2}{n}(\bar{x}-\mu)\sum_{i=1}^n(x_i-\mu) + (\bar{x}-\mu)^2)$$ 注意到 $$\sum_{i=1}^n(x_i-\mu) = n\bar{x}-n\mu$$ 所以 $$\mathbb{E}(\frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2 - \frac{2}{n}(\bar{x}-\mu)\sum_{i=1}^n(x_i-\mu) + (\bar{x}-\mu)^2) \\ = \mathbb{E}(\frac{1}{n}\sum_{i=1}^n(x_i-\mu)^2 - (\bar{x}-\mu)^2) \\ = \sigma^2 - \mathbb{E}((\bar{x}-\mu)^2)$$ 注意到 $$\mathbb{E}(\bar{x}) = \mu$$ 展开第二项有 $$\mathbb{E}((\bar{x}-\mu)^2) = \mathbb{E}(\bar{x}^2) - 2\mu\mathbb{E}(\bar{x}) + \mu^2 \\ = \mathbb{E}(\bar{x}^2) - \mathbb{E}^2(\bar{x}) = \text{Var}(\bar{x})$$ $\bar{x}$ 是 $n$ 个相互独立且方差均为 $\sigma^2$ 的变量的均值，所以 $$\text{Var}(\bar{x}) = \frac{\sigma^2}{n}$$ 所以 $$\mathbb{E}(\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2) = \sigma^2 - \mathbb{E}((\bar{x}-\mu)^2) \\ = \frac{n-1}{n}\sigma^2 \ne \sigma^2$$ 这就是以 $1/N$ 为系数的方差是有偏估计的原因。相应的 $$\mathbb{E}(\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})^2) \\ = \frac{n}{n-1}\mathbb{E}(\frac{1}{n}\sum_{i=1}^n(x_i-\bar{x})^2) \\ = \frac{n}{n-1}\frac{n-1}{n}\sigma^2 = \sigma^2$$ 所以以 $1/(N-1)$ 为系数的方差才是无偏估计。
回到最开始的问题，为了说明 $\hat{\sigma}^2$ 是对 $\sigma^2$ 的无偏估计，我们需要证明 $\mathbb{E}(\hat{\sigma}^2) = \sigma^2$。这里设 $\epsilon$ 是一个协方差矩阵为 $\sigma^2I$ 的向量。$$\mathbb{E}(\sum_{i=1}^N(y_i-\hat{y}_i)^2) = \mathbb{E}(|Y-X\hat{\beta}|^2) \\ = \mathbb{E}(|X\beta + \epsilon - X(X^TX)^{-1}X^T(X\beta+\epsilon)|^2) \\ = \mathbb{E}(|X\beta + \epsilon - X\beta - X(X^TX)^{-1}X^T\epsilon|^2) \\ = \mathbb{E}(|(I-X(X^TX)^{-1}X^T)\epsilon|^2)$$ 令 $X(X^TX)^{-1}X^T = H$，容易验证 $H^T = H^2 = H$。我们有 $$\mathbb{E}(|(I-X(X^TX)^{-1}X^T)\epsilon|^2) \\ = \mathbb{E}(\epsilon^T(I-H)^T(I-H)\epsilon) \\ = \mathbb{E}(\epsilon^T(I-H)\epsilon) = \sum_{i=1}^N\sum_{j=1}^N(I-H)_{ij}\mathbb{E}(\epsilon_i\epsilon_j)$$ 注意到除非 $i = j$，否则 $\epsilon_i$ 与 $\epsilon_j$ 互相独立，且 $\mathbb{E}(\epsilon) = 0$，所以 $$\sum_{i=1}^N\sum_{j=1}^N(I-H)_{ij}\mathbb{E}(\epsilon_i\epsilon_j) \\ = \sum_{i=1}^N(I-H)_{ii}(\mathbb{E}(\epsilon_i^2) - \mathbb{E}^2(\epsilon_i)) \\ = \sigma^2(N-\text{tr}(H)) \\ = \sigma^2(N-\text{tr}(X(X^TX)^{-1}X^T)) \\ = \sigma^2(N-\text{tr}(X^TX(X^TX)^{-1})) = (N-p-1)\sigma^2$$ 所以 $$\mathbb{E}(\hat{\sigma}^2) = \frac{1}{N-p-1}\mathbb{E}(\sum_{i=1}^N(y_i-\hat{y}_i)^2) = \sigma^2$$