7.2 线性回归与最小二乘

数据在我们的生活中随处可见. 我们想要通过建立模型来捕捉变量之间的关系, 并且用我们的模型进行预测. 我们假设获得的一组数据为

$$D_n := \{ (\vec x_1, y_1), (\vec x_2,y_2),\cdots,(\vec x_n,y_n) \},$$

而我们的目标便是尝试建立$\vec x_i$与$y_i$之间的联系. 我们将$\vec x_i= \begin{pmatrix} x_{i1} & x_{i2} & \cdots & x_{im} \end{pmatrix}^\top \in \mathbb{R}^m$称作是协变量向量 (Covariates) (也可叫做特征 (Features)), 其中$x_{ij}$代表数据$\vec x_i$的第$j$个协变量. 我们还定义$y_i \in \mathbb{R}$为响应变量 (Response Variable), 即我们感兴趣的研究对象. 协变量和响应变量的选择取决于我们的研究对象和研究目的. 在线性回归下, 我们希望建立响应变量$y_i$与协变量$\vec x_i$之间的线性关系. 即我们希望求出常数$\beta_0,\beta_1,\cdots,\beta_m \in \mathbb{R}$, 使得

$$y_i = \beta_0 + \beta_1 x_{i1} +\beta_2 x_{i2} +\cdots+\beta_mx_{im}+\epsilon_i,\tag{7.1}$$

其中$\beta_0$被称作该模型的截距 (Intercept), $\epsilon_i$被我们称作是模型的噪声 (Noise), 即无法用线性模型解释的波动因素, 它可以视作为误差, 噪声, 遗漏因素或随机扰动. 我们假设这些噪声为相互独立的随机变量, 且它们的均值为$0$, 方差为$\sigma^2_i$. 对于这些相互独立的噪声而言, 它们可以简要分为两种情况: 同方差性 (Homoskedasticity) 与 异方差性 (Heteroskedasticity). 对于同方差性而言, 我们假设所有的噪声均为独立同分布的样本, 即它们的方差为一常数$\sigma^2$. 而异方差性假设这些噪声的方差互不相同. 在接下来的讨论中, 我们默认模型的噪声满足同方差性.

如果我们考虑对所有的$(\vec x_i,y_i)$使用等式(7.1), 我们便可以得到一个线性方程组:

$$\begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} = \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1m} \\ 1 & x_{21} & x_{22} & \cdots & x_{2m} \\ \vdots & \vdots & \vdots & &\vdots\\ 1 & x_{n1} & x_{n2} & \cdots & x_{nm} \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_m \end{pmatrix}+\begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{pmatrix}.$$

我们也将上式简写为$\vec y = \vec X\vec\beta+\vec\epsilon$或$y_i = \vec x_i^\top\vec\beta+\epsilon_i$, 其中$\vec y$为响应向量 (Response Vector), $\vec X \in M_{ n\times(m+1)}(\mathbb{R})$为设计矩阵 (Design Matrix), $\vec\beta$为参数向量 (Coefficient Vector), $\vec\epsilon$为噪声向量 (Noise Vector). 我们的目标即为求解$\vec\beta$, 从而便能得到$y_i$与$\vec x_i$之间的关系. 在本节的学习中, 我们规定设计矩阵$\vec X$满足固定设计 (Fixed Design), 即$\vec X$中的元素均为已知量, 模型中的随机性仅来源于噪声;我们同时规定设计矩阵$\vec X$满列秩, 即$\vec X$的列向量彼此线性无关. 那么此时由于$\vec x_i$的维数为$m+1$, 它们之间彼此线性无关也就说明了$m+1\leq n$. 即样本点要足够多. 这对我们随后构建最小二乘有着至关重要的意义.

*{最小二乘估计量及其性质}

我们发现, 如果$m+1<n$, 那么便代表了在$\vec y = \vec X\vec\beta+\vec\epsilon$中方程的数量多于未知数的数量. 回顾第一章, 第一节的知识, 此时的方程很可能是无解的. 那么我们能不能尝试求出一个“近似解”呢? 即求出满足$\vec y$与$\vec X\beta$非常接近的参数向量. 在向量中, “非常接近”便对应了向量差值的模非常小. 因此我们希望$||\vec y - \vec X\beta||$非常小, 即噪声非常小 . 由此我们也就引出了最小二乘 (Least Squares): 即求出优化问题 $\arg\min_{\vec\beta \in \mathbb{R}^{m+1}}||\vec y - \vec X\vec\beta||^2$ 的解$\hat{\vec\beta}$.

我们做一个大致证明: 严格的证明依赖于严格凸函数的性质.

此时的$\hat{\vec\beta}$被我们称作普通最小二乘估计量 (Ordinary Least Squares Estimator), 或简称为OLS估计量. 我们有时候将其记作$\hat{\vec\beta}_{OLS}$, 或$\hat{\vec\beta}$. 求解完$\hat{\vec\beta}$, 我们便可以用得到的参数来对我们的模型进行检验和预测.

该定理的证明涉及到一些没有在本章第一节提及的概率论知识, 感兴趣的读者欢迎阅读.

实际上, 我们可以根据条件无偏性得到无偏性. 这是因为$\E[\hat{\vec\beta}] = \E[\E[\hat{\vec\beta} | \vec X]]$, 因此我们也称$\hat{\vec\beta}$为一无偏估计量. 对于方差而言, 根据公式

$$\Var(\hat{\vec\beta}) = \E[\Var(\hat{\vec\beta}|\vec X)] + \Var(\E[\hat{\vec\beta} | \vec X]) = \sigma^2 \E\Big[ \Big(\vec X^\top \vec X\Big)^{-1} \Big] = \sigma^2\cdot(\vec X^\top\vec X)^{-1},$$

我们可以得到相似的结果. 在所有的无偏估计量中, 我们的目标是寻找方差最小的一个估计量. 这一估计量其实就是最小二乘估计量. 该定理被称作Gauss-Markov定理 (Gauss-Markov Theorem):

该定理的证明略.

在模型$\vec y = \vec X\vec\beta+\vec\epsilon$中, 数据的残差 (Residual) 为$\epsilon_i = y_i - \vec x_i^\top\vec\beta$. 对于最小二乘估计量$\hat{\vec\beta}$, 我们称$\hat{y}_i = \vec x_i^\top \hat{\vec\beta}$为数据点$(\vec x_i,y_i)$的拟合值 (Fitted Value). 观测值$y_i$和拟合值$\hat{y}_i$的差$\hat{\epsilon_i} = y_i - \vec x_i^\top \hat{\vec\beta}$被称作该数据点的残差估计量 (Estimated Residual). 那么如果我们有一个新的数据点$\vec x_0$, 我们便可以用公式$\hat{y}_0 = \vec x_0^\top \vec\beta$求出此时响应变量的拟合值. 随后我们便可以构建置信区间来检验拟合的精确性. 由于本书中缺乏对概率论与数理统计的深入讨论, 因此有关置信区间的构建我们在此不作涉及.

当我们初步了解完最小二乘的性质之后, 我们给出两个矩阵的定义: 它们分别为帽子矩阵 (Hat Matrix) 和残差生成矩阵 (Annihilator Matrix).

对于帽子矩阵而言, 我们不难发现

$$\vec H\vec y = \vec X(\vec X^\top\vec X)^{-1}\vec X^\top\vec y = \vec X\hat{\vec\beta} = \hat{\vec y},$$

而对于残差生成矩阵, 同样地我们有

$$\vec M\vec y = \vec y - \hat{\vec y} = \hat{\vec \epsilon}.$$

本节的课后练习题中会有更多关于帽子矩阵和残差生成矩阵的性质等待读者去探索. 在普通最小二乘中，我们通过最小化 $\sum_{i=1}^n (y_i-x_i^\top\vec\beta)^2$ 来选择参数$\vec\beta$. 这相当于认为每一个观测点的误差$\epsilon_i$都具有相同的重要性. 然而在实际问题中, 不同数据点的可靠程度可能并不相同. 例如:有些观测值测量误差较小,因此相对更值得信任; 有些观测值噪声较大,因此不应对模型产生过强影响. 为了反映这种差异, 我们可以引入一个权函数 (Weight Function), 使得在观测点$\vec x_i$处的权重为$w_i>0$. 权重越大, 则表示该观测点在拟合中越重要;权重越小,则表示该观测点对最终模型的影响较弱. 因此我们此时的最小二乘优化的便是 $\sum_{i=1}^n w_i(y_i - \vec x_i^\top\vec\beta)^2$ 的最小值.

不难发现, 一般情况下的线性回归即对应了权函数$w_i\equiv 1$. 最小二乘有什么几何意义呢? 由于我们知道$n>m+1$, 也就是说$\vec y$的维数大于$\vec X$的列空间的位数. 此时可以将最小二乘其视作响应变量$\vec y$在设计矩阵$\vec X$的列空间上的正交投影. 此时$\vec y$满足正交分解$\vec y = \vec X\vec\beta + (\vec y - \vec X\vec\beta)$, $\vec y - \vec X\vec\beta \in \image(\vec X)^\perp$. 回顾我们在内积空间中提出的定理: $\image(\vec X)^\perp = \Ker(\vec X^\top)$. 因此$\vec y - \vec X\vec\beta \in \Ker(\vec X^\top)$, 即

$$\vec X^\top (\vec y - \vec X\vec\beta)=\vec 0.$$

将其展开我们便得到了系统的正规方程: $\vec X^\top\vec y = \vec X^\top\vec X\vec\beta$. 当设计矩阵满列秩时我们便得到了前面推出的$\hat{\vec\beta} = (\vec X^\top\vec X)^{-1}\vec X^\top\vec y$. 下面的图便很好地展示了这一几何关系:

*{多项式回归}

值得我们注意的是, 线性回归仅代表我们的模型关于参数向量$\vec\beta$是线性的, 而关于协变量$\vec x$则不需要满足线性关系. 比如我们完全可以设$\vec x = \begin{pmatrix} x & x^2 & \cdots & x^m \end{pmatrix}^\top$, 这样一来我们得到的模型便为

$$y = \beta_0+\beta_1x+\beta_2 x^2+\cdots+\beta_m x^m,$$

其中$y$是关于$x$的多项式. 我们也把这类问题称作是多项式回归 (Polynomial Regression). 给出数据点$D_n = \{ (x_1, y_1),\cdots,(x_n,y_n)\}$, 假设我们将使用关于$x$的$m$次多项式去估计$y$, 那么此时我们有

$$\begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} = \begin{pmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^m \\ 1 & x_2 & x_2^2 & \cdots & x_2^m\\ \vdots & \vdots & \vdots & &\vdots \\ 1 & x_n & x_n^2 & \cdots & x_n^m \end{pmatrix}\begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_m \end{pmatrix}+\begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{pmatrix}.$$

此时的的最小二乘为

$$\sum_{i=1}^n |y_i - (\beta_0+\beta_1x+\beta_2 x^2+\cdots+\beta_mx^m)|.$$

类似地, 如果我们假设的模型为

$$y= \beta_0 f_0(x) + \beta_1 f_1(x)+\cdots+\beta_m f_m(x),$$

那么我们有

$$\begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix} = \begin{pmatrix} f_0(x_1) & f_1(x_1) & f_2(x_1) & \cdots & f_m(x_1) \\ f_0(x_2) & f_1(x_2) & f_2(x_2) & \cdots & f_m(x_2)\\ \vdots & \vdots & \vdots & &\vdots \\ f_0(x_n) & f_1(x_n) & f_2(x_n) & \cdots & f_m(x_n) \end{pmatrix}\begin{pmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_m \end{pmatrix}+\begin{pmatrix} \epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n \end{pmatrix}.$$

此时我们考虑多项式回归中的一个特殊情况: 假设$m=1$, 此时的最小二乘可以写作

$$\begin{aligned} (\hat{\beta}_0,\hat{\beta}_1)&=\arg\min_{\beta_0,\beta_1} \sum_{i=1}^n\{ y_i - (\beta_0+\beta_1x_i)\}^2\\ &=\arg\min_{\beta_0,\beta_1}\frac{1}{n} \sum_{i=1}^n\{y_i^2 +\beta_0^2+2\beta_0\beta_1 x_i+\beta_1^2x_i^2 -2y_i\beta_0-2\beta_1y_ix_i \}. \end{aligned}$$

我们令

$$\frac{1}{n}\sum_{i=1}^n x_i := \bar{x}_n \hspace{1cm} \frac{1}{n} \sum_{i=1}^n y_i :=\bar{y}_n,$$

由此我们进一步将原式化简为

$$(\hat{\beta}_0,\hat{\beta}_1) := \arg\min_{\beta_0,\beta_1} \left\{ \beta_0^2+\frac{1}{n}\sum_{i=1}^n x_i^2 \beta_1^2 +2\bar{x}_n\beta_0\beta_1-2\bar{y}_n \beta_0 - \frac{2}{n}\sum_{i=1}^n x_i y_i \beta_1 \right\}.$$

我们定义关于$\beta_0,\beta_1$的二元函数

$$f(\beta_0,\beta_1)=\beta_0^2+\frac{1}{n}\sum_{i=1}^n x_i^2 \beta_1^2 +2\bar{x}_n\beta_0\beta_1-2\bar{y}_n \beta_0 - \frac{2}{n}\sum_{i=1}^n x_i y_i \beta_1,$$

因此$f$的梯度为

$$\nabla f =\left( \frac{\partial f}{\partial \beta_0}, \frac{\partial f}{\partial \beta_1} \right)= \left( 2\beta_0+2\bar{x}_n \beta_1-2\bar{y}_n\hspace{0.2cm}, \hspace{0.2cm}\frac{2}{n}\sum_{i=1}^n x_i^2 \beta_1 +2\bar{x}_n\beta_0-\frac{2}{n}\sum_{i=1}^n x_iy_i\right).$$

令$\nabla f = 0$, 则我们有

$$\begin{cases} 2\beta_0+2\bar{x}_n \beta_1-2\bar{y}_n &=0\\ \frac{2}{n}\sum_{i=1}^n x_i^2 \beta_1 +2\bar{x}_n\beta_0-\frac{2}{n}\sum_{i=1}^n x_iy_i&=0 \end{cases}.$$

其中不难发现

$$\hat{\beta}_0 = \bar{y}_n - \bar{x}_n\hat{\beta}_1. \tag{7.2}$$

然后利用(7.2)中$\hat{\beta}_0$的取值带入另外一式求出$\beta_1$:

$$\frac{1}{n}\sum_{i=1}^n x_i^2 \hat{\beta}_1 + \bar{x}_n (\bar{y}_n - \bar{x}_n\hat{\beta}_1) - \frac{1}{n}\sum_{i=1}^n x_i y_i = 0,$$

即

$$\hat{\beta}_1 = \frac{\sum_{i=1}^n x_i y_i - n\bar{x}_n \bar{y}_n}{\sum_{i=1}^n (x_i - \bar{x}_n)^2}.$$

我们将这一发现总结成下面的推论:

此时读者会想: 如果我们有$n$组数据$(x_1,y_1),\cdots,(x_n,y_n)$, 且这些$x_i$互不相同, 那么根据Lagrange多项式, 我们可以找到一个$(n-1)$次多项式$p(x)$, 使得$p(x_i)=y_i$. 也就是说, 理论上完全存在这样的函数, 使得这个函数能够与数据完美吻合, 那么我们还干嘛费尽心思地去用最小二乘求解呢? 我们不妨来看下图:

在回归分析中, 我们关注的不仅仅是模型在以观测到的数据集上的表现, 我们还要求模型要有精准的预测能力. 在Lagrange多项式拟合的结果中, 虽然拟合函数精确地穿过了每一个数据点, 做到了零误差, 但我们看到得到的函数在数据点之间会来回大幅度地波动, 并且在$x$坐标相距很近的情况下对应的响应变量$y$却有天壤之别. 因此Lagrange多项式的拟合结果显然没有捕捉到数据之间的大致变动趋势. 相反, 红色的线性函数虽然没有精确地穿过每一个数据点, 但它在一定程度上反映出了$x,y$之间的变化趋势. 这个例子说明: 复杂的模型并不一定意味着预测效果越好. 若模型过度追逐训练数据中的偶然波动或噪声, 就可能在已有数据上表现很好, 却在新的数据上表现较差. 这种现象称为过拟合 (Overfitting).

为避免过拟合的发生, 我们通常把数据集分成训练集 (Training Set)和检验集 (Testing Set). 即我们利用训练集中的数据构造设计矩阵, 进行回归分析. 然后再用得到的参数估计去拟合检验集中的数据, 然后比较误差. 不过, 当数据较少时, 单纯划分一次训练集和检验集会产生一定程度的数据浪费, 因此我们也可以使用交叉验证 (Cross Vaildation)的方法. 交叉验证的基本逻辑是通过重复实验, 多次改变训练集和检验集, 使得每一个数据都能被最大程度地利用. 常见的交叉验证方法有K折交叉验证 (K-fold Cross Validation), 即将数据分成$K$份, 其中$K-1$份作为训练集, 剩下的一份作为检验集, 然后重复实验$K$次, 取模型的平均检验误差. 这样一来每一份数据都能够作为检验集, 通过平均, 这样得到的误差通常比单次划分更加稳定.

此时, 如果我们的数据如果不再是离散的数据点$(x_i,y_i)$, 而是一整个连续的区间$(a,b)$, 那么在这个区间上我们能否也类比最小二乘的知识呢? 此时我们会联想到积分的知识. 即当我们在区间$(a,b)$上用多项式$p(x) = \beta_0+\beta_1x+\cdots+\beta_mx^m$去近似$y=f(x)$时, 我们所要优化的便是积分

$$f(\beta_0,\cdots,\beta_m) = \int_a^b | f(x) - p(x)|^2dx = \int_a^b \left( f(x) - \sum_{i=0}^m \beta_i x^i \right)^2 dx.$$

我们此时求出$f$关于$\beta_i$的偏导数, 得到

$$\begin{aligned} \frac{\pl f}{\pl \beta_i} &= \frac{\pl}{\pl\beta_i}\int_a^b\left( f(x) - \sum_{i=0}^m \beta_i x^i \right)^2 dx\\ &= \int_a^b \frac{d}{d\beta_i}\left[ f^2(x) - 2f(x)\sum_{i=0}^m \beta_i x^i + \left( \sum_{i=0}^m \beta_i x^i \right)^2 \right]dx\\ &=\int_a^b -2f(x)x^i+2x^i\sum_{j=0}^m\beta_jx_j dx. \end{aligned}$$

在极小值处, 我们令$\frac{\pl f}{\pl \beta_i}=0$, 我们因此得到此时的正规方程:

$$\int_a^b x^if(x) dx = \sum_{j=0}^m \beta_j \int x^{i+j} dx, \quad i=0,1,\cdots,m.$$

由此我们便得到了如下的矩阵形式:

$$\begin{pmatrix} \int_a^b f(x)\\ \int_a^b x f(x) dx \\ \vdots \\ \int_a^b x^mf(x)dx \end{pmatrix} = \begin{pmatrix} \int_a^b 1 dx & \int_a^b xdx & \int_a^b x^2 dx & \cdots & \int_a^b x^m dx\\ \int_a^b xdx & \int_a^b x^2 dx & \int_a^b x^3 dx & \cdots & \int_a^b x^{m+1}dx\\ \vdots & \vdots & \vdots & &\vdots\\ \int_a^b x^m dx & \int_a^b x^{m+1}dx & \int_a^b x^{m+2}dx & \cdots & \int_a^b x^{2m} dx \end{pmatrix} \begin{pmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \\ \vdots \\ \beta_m \end{pmatrix},$$

其中等式右边的$(m+1)\times(m+1)$矩阵被称作Gram 矩阵. Gram矩阵的列向量彼此线性无关(这一点读者不妨自行验证), 因此上面的系统便存在唯一解. 此时对于确定的$f(x),m,a,b$, 我们便可以求出每一个积分的值. 随后我们便可以运用求解线性方程组的方法去求解$\beta_0,\cdots,\beta_m$了.

{7.2 练习}

1. 设$\vec X \in M_{n(m+1)}(\mathbb{R})$为设计矩阵, 帽子矩阵为$\vec H = \vec X(\vec X^\top\vec X)^{-1}\vec X^\top$; 残差生成矩阵为$\vec M = \vec I - \vec H$.

(i) 证明: $\vec H, \vec M$均为对称矩阵;

(ii) 证明: $\vec H, \vec M$均为投影矩阵, 即$\vec H^2 = \vec H, \vec M^2 = \vec M$;

(iii) 证明: $\tr(\vec H)=m+1$, $\tr(\vec M) = n-m-1$.

(iv) 对于样本$\vec x_i$而言, 我们定义其杠杆值 (Leverage)为$h_{ii} = \vec x_i^\top (\vec X^\top \vec X)^{-1}\vec x_i$. 证明: $0 \leq h_{ii} \leq 1$.

(v) 我们假设样本中的噪声满足同方差性, 证明$\Var(\hat{\epsilon}_i) = \sigma^2(1-h_{ii})$.

2. 设$\vec X \in M_{n(m+1)} (\mathbb{R})$为设计矩阵, 我们假设$\vec X$从左往右的的第一列为截距列， 即第一列所有元素为$1$. 样本$\vec x_i$的第一个协变量也为截距项$1$.

(i) 对于任意的样本$\vec x_i$, 证明$\sum_{j=1}^n \vec x_i^\top (\vec X^\top \vec X)^{-1} \vec x_j =1$.

(ii) 证明

$$\sum_{i=1}^n \hat{\epsilon}_i = 0.$$

3. 设$\vec X \in M_{n(m+1)}$为设计矩阵, $\hat{\vec \epsilon} = \vec y - \vec X\hat{\vec\beta}$, $\vec\epsilon = \vec y - \vec X\vec\beta$. $\vec H$为帽子矩阵, $\vec M$为残差生成矩阵. 证明$\hat{\vec\epsilon}^\top\hat{\vec\epsilon} = \vec\epsilon^\top\vec M\vec\epsilon$.

4. 当我们假设线性模型$y_i = \vec x_i^\top\vec\beta + \epsilon_i$时, 若噪声满足同方差性且服从正态分布$\epsilon_i \sim N(0,\sigma^2)$, 我们可以将样本的似然函数写作

$$L(\vec\beta, \sigma^2) = \frac{1}{(2\pi\sigma^2)^{n/2}} \cdot\exp\left( -\sum_{i=1}^n \frac{(y_i - \vec x_i^\top \vec\beta)^2}{2\sigma^2} \right).$$

求出在此情况下$\vec\beta$和$\sigma^2$的最大似然估计.

5. 在以下的回归分析问题中, 我们假设有且仅有一个分类项作为协变量. 分类项中我们假设含有$K$个不同的种类. 定义$\mathbb{1}_{\vec x_i = g}$为样本$\vec x_i$是否属于种类$g$的示性函数 (Indicator Function), 即

$$\mathbb{1}_{\vec x_i=g} = \begin{cases} 1 & \text{ 若样本} \vec x_i \text{属于种类} g \\ 0 & \text{ 若样本} \vec x_i \text{不属于种类} g \end{cases}.$$

我们定义$\vec z_i = \begin{pmatrix} \mathbb{1}_{\vec x_i=1} \\ \mathbb{1}_{\vec x_i=2} \\ \vdots \\ \mathbb{1}_{\vec x_i = K} \end{pmatrix}$, 设计矩阵$\vec Z=\begin{pmatrix} \vec z_1^\top \\ \vec z_2^\top \\ \vdots \\ \vec z_n^\top \end{pmatrix}$. 设$n_g$为属于种类$g$的样本数, 证明此时在模型$\vec y = \vec Z\vec\beta+\vec\epsilon$中$\vec\beta$的最小二乘估计量$\hat{\vec\beta}$为

$$\hat{\vec\beta} = \begin{pmatrix} \frac{1}{n_1}\sum_{y \in n_1} y \\ \frac{1}{n_2} \sum_{y \in n_2} y \\ \vdots \\ \frac{1}{n_K} \sum_{y \in n_K} y \end{pmatrix}.$$

6. 在加权线性回归中, 给出设计矩阵$\vec X$和权函数矩阵$\vec W$我们得到的关于$\vec\beta$的加权最小二乘估计量为

$$\hat{\vec\beta}_{GLS} = \Big( \vec X^\top \vec W \vec X\Big)^{-1}\vec X^\top \vec W\vec y.$$

假设噪声相互独立,现给出两个不同条件: 噪声满足同方差性 (M1); 噪声满足异方差性 (M2).

(i) 证明: 在(M1), (M2)条件下均有 $\E[\hat{\vec\beta}_{GLS}|\vec X] = \vec \beta$;

(ii) 在(M1)条件下证明

$$\Var(\hat{\vec\beta}_{GLS} | \vec X) = \sigma^2 (\vec X^\top\vec X)^{-1};$$

(iii) 在(M2)条件下, 假设$\Var(\epsilon_i|\vec X) = \sigma_i^2$, 证明

$$\Var(\hat{\vec\beta}_{GLS}|\vec X) = (\vec X^\top\vec X)^{-1} \vec X^\top \diag\begin{pmatrix} \sigma_1^2 & \sigma_2^2 & \cdots & \sigma_n^2 \end{pmatrix}.$$

(iv) 在(M2)条件下, 我们假设$\Var(\vec\epsilon | \vec X) = \sigma^2 \vec\Sigma$, 其中$\vec\Sigma$为对角矩阵. 设模型

$$\widetilde{\vec y} = \widetilde{\vec X} \vec\beta + \widetilde{\vec\epsilon},$$

其中$\widetilde{\vec y} = \vec\Sigma^{-1/2} \vec y$, $\widetilde{\vec X} = \vec\Sigma^{-1/2} \vec X$, $\widetilde{\vec\epsilon} = \vec\Sigma^{-1/2} \vec\epsilon$. 此时的最小二乘估计为$\widetilde{\vec\beta} = (\widetilde{\vec X}^\top \widetilde{\vec X})^{-1} \widetilde{\vec X}^\top \widetilde{\vec y}$. 求出$\E[\widetilde{\vec\beta}|\widetilde{\vec X}]$与$\Var(\widetilde{\vec\beta} | \widetilde{\vec X})$.

7. 在线性模型$y=\vec x^\top\vec\beta+\epsilon$中我们假设噪声满足同方差性, 设计矩阵为固定设计且满列秩. 现有观测数据

$$D=\{ (\vec x_1, y_1), (\vec x_2, y_2),\cdots, (\vec x_n,y_n)\},$$

我们将采用留一交叉验证 (Leave One Out Cross Validation)对最小二乘估计量进行分析. 在留一交叉验证中, 我们将原数据复制$n$份, 在每一份中采用$n-1$个数据作为训练集, 随后剩余的一个数据作为检验集. 这样一来每一个数据都可以作为检验集. 我们定义集合$E_i \subset D$为除去原数据中除去第$i$个样本之后得到的新数据集, 即$E_i = D\setminus \{ (\vec x_i,y_i)\}$. 那么在数据集$E_i$上关于$\vec\beta$的最小二乘估计量为

$$\hat{\vec\beta}^{(i)} = \arg\min_{\vec\beta} \frac{1}{n-1}\sum_{j \neq i} (y_i - \vec x_i^\top\vec\beta)^2.$$

(i) 利用Sherman-Morrison公式:

$$(\vec A+\vec u \vec v^\top )^{-1} = \vec A^{-1} - \frac{\vec A^{-1} \vec u \vec v^\top \vec A^{-1}}{1+\vec v^\top\vec A^{-1}\vec u},$$

证明此时的最小二乘估计量$\hat{\vec\beta}^{(i)}$满足

$$\hat{\vec\beta}^{(i)} = \hat{\vec\beta} - \frac{(\vec X^\top\vec X)^{-1}\vec x_i \hat{\epsilon}_i}{1-h_{ii}}.$$

其中$\displaystyle{\hat{\vec\beta} = (\vec X^\top\vec X)^{-1}\vec X^\top\vec y = \left(\sum_{i=1}^n \vec x_i \vec x_i^\top \right)^{-1} \cdot \left( \sum_{i=1}^n \vec x_i y_i \right)}$为普通最小二乘估计量, $\hat{\epsilon}_i = y_i - \vec x_i^\top \hat{\vec\beta}$.

(ii) 证明: $\displaystyle{\hat{\epsilon}_i^{(i)} = y_i - \vec x_i^\top \hat{\vec\beta}^{(i)} = \frac{\hat{\epsilon}_i}{1- h_{ii}}}$.

8. 在推论$7.1$中, 证明$\hat{\beta}_1$的估计还可以写作

$$\hat{\beta}_1 = \frac{\sum_{i=1}^n(x_i - \bar{x}_n)(y_i - \bar{y}_n)}{\sum_{i=1}^n(x_i - \bar{x}_n)^2}.$$ 9.

(i) 设函数$f(x)=\frac{1}{x}$, $1\leq x\leq 3$. 利用Gram矩阵, 求出在该区间上满足最小二乘条件的二次多项式$p(x)=\beta_0+\beta_1x+\beta_2x^2$;

(ii) 设$A=(-1,0),B=(0,1),C=(1,4)$. 以$y = \beta_0 x + \beta_1 2^x$为模型, 利用最小二乘求出此时的参数$\beta_0,\beta_1$.