1. 模型求解
假设输入数据为\vec{x}=\left[\begin{array}{c}1 & x_1 & x_2 & \cdots & x_k \end{array}\right]，共具有k个特征，线性回归模型可写为
\LARGE{y=\vec{x}\cdot\vec{\beta}+\epsilon}其中\vec{\beta}=\left[\begin{array}{c}\beta_{0} \\ \beta_{1} \\ \vdots \\ \beta_{k} \end{array}\right]为待求系数，\epsilon\sim{N(0,\sigma^2)}表示误差。从统计学上来说，这个模型可以表示为\Large{y|\vec{x}\text{ }\sim{N(\vec{x}\cdot\vec{\beta}, \sigma^2)}}其中\vec{\beta}以及\sigma^2为统计分布的参数。在文章这一部分中主要解决以下问题：
假定共有n个训练数据，则输入数据集可记为X=\left[\begin{array}{c} \vec{x}^{(1)} \\ \vec{x}^{(2)} \\ \vdots \\ \vec{x}^{(n)}\end{array}\right]=\begin{bmatrix} 1 & x^{(1)}_1 & x^{(1)}_2 & \cdots & x^{(1)}_k \\ 1 & x^{(2)}_1 & x^{(2)}_2 & \cdots & x^{(2)}_k \\ & & \vdots \\ 1 & x^{(n)}_1 & x^{(n)}_2 & \cdots & x^{(n)}_k \end{bmatrix}对应的n个输出值可记为\vec{y}=\left[\begin{array}{c}y_1 \\ y_2 \\ \vdots \\ y_n\end{array}\right]=X\vec{\beta}+\vec{\epsilon}=X\vec{\beta}+\left[\begin{array}{c}\epsilon_1 \\ \epsilon_2 \\ \vdots \\ \epsilon_n\end{array}\right]其中\epsilon_i相互独立并且\epsilon_i\sim{N(0,\sigma^2)}。接下来使用最大似然估计来计算\vec{\beta}，由于误差相互独立且为正态分布，概率公式可以表示为\large{p(y_1,\cdots,y_n|\vec{x}^{(1)},\cdots,\vec{x}^{(n)})=\prod_{i=1}^n{p(y_i|\vec{x}^{(i)})}\propto{e}^{-\frac{1}{2\sigma^2}\color{orange}{\sum_{i=1}^n{(y_i-\vec{x}^{(i)}\cdot\vec{\beta})^2}}}}若使该式最大，需使黄色部分最小，这里为方便后续计算写成矩阵形式：\Large{\argmin_{\vec \beta}[(\vec{y}-X\vec{\beta})^T(\vec{y}-X\vec{\beta})]}使式子的一阶导为0，即X^T(\vec{y}-X\vec{\beta})=0。\vec{\beta}的估计值为\Large{\hat{\vec{\beta}}=(X^TX)^{-1}X^T\vec{y}}
获得参数\vec{\beta}的估计值\hat{\vec{\beta}}后，\vec{y}的估计值可写为\large{\hat{\vec{y}}=X\hat{\vec{\beta}}=X(X^TX)^{-1}X^T\vec{y}=H\vec{y}}其中H又被称为Hat Matrix，在后续的分析中会经常用到这一矩阵，因此在这里先介绍一下它的几点性质：
- 是一个对称矩阵，即H^T=H
- H^2=H
- 矩阵迹tr(H)=\sum_{i=1}^nh_{ii}=tr(X^TX(X^TX)^{-1})=k+1，h_{ii}为对角线上的元素
由第一条性质可以对矩阵进行特征值分解\large{H=V\Lambda{V^T}\text{, }V=\left[\begin{array}{c} \vec{v}_1 & \vec{v}_2 & \cdots & \vec{v}_n \end{array}\right]}其中V为正交方阵，满足V^TV=VV^T=I，它的列向量\vec{v}_i为H的特征向量，\Lambda为对角阵，对角线上的元素\lambda_i为H的特征值。由第二条性质可得\large{\vec{v}_i^TH\vec{v}_i=\lambda_i\|\vec{v}_i\|_2^2=\vec{v}_i^TH^2\vec{v}_i=\|H\vec{v}_i\|_2^2=\lambda_i^2\|\vec{v}_i\|_2^2}容易看出H的特征值只能为0或1。由第三条性质可得\large{tr(H)=tr(V\Lambda{V^T})=tr(V^TV\Lambda)=tr(\Lambda)=\sum_{i=1}^n\lambda_i=k+1}综合以上三条性质可以将H的特征值分解写为\large{H=\left[\begin{array}{c} V_1 & V_2 \end{array}\right]\begin{bmatrix}I_{k+1} & O \\ O & O \end{bmatrix}\left[\begin{array}{c} V_1^T \\ V_2^T \end{array}\right]=V_1V_1^T}其中V_1=\left[\begin{array}{c} \vec{v}_1 & \vec{v}_2 & \cdots & \vec{v}_{k+1} \end{array}\right]\text{, }V_2=\left[\begin{array}{c} \vec{v}_{k+2} & \vec{v}_{k+3} & \cdots & \vec{v}_n \end{array}\right]，O为零矩阵。容易看出V_1和V_2满足以下性质：\large{\begin{cases}V_1^TV_1=I_{k+1} \\ V_2^TV_2=I_{n-k-1} \\ V_1^TV_2=O \\ V_2^TV_1=O\end{cases}}（注：这些零矩阵的行数和列数并不相同，只是为了方便统一记为O，但并不影响后续分析）
2. 参数估计的统计分布
文章这一部分主要解决以下问题：
- 参数\vec{\beta}的估计值\hat{\vec{\beta}}的统计分布的均值和协方差矩阵
- 如何计算参数\sigma^2的估计值s^2以及它的统计分布
- 参数估计值\hat{\vec{\beta}}与s^2的相关性
首先计算参数估计值\hat{\vec{\beta}}的均值：\large{E(\hat{\vec{\beta}})=[(X^TX)^{-1}X^T]E(\vec{y})=[(X^TX)^{-1}X^T]X\vec{\beta}=\vec{\beta}}
接着计算参数估计值\hat{\vec{\beta}}的协方差矩阵：\large{Cov(\hat{\vec{\beta}})=[(X^TX)^{-1}X^T]Cov(\vec{y})[(X^TX)^{-1}X^T]^T=\sigma^2(X^TX)^{-1}}由正态分布的性质以及上述两式容易看出\Large{\hat{\vec{\beta}}\text{ }\sim{N(\vec{\beta}, \sigma^2(X^TX)^{-1})}}
现在对我们来说并不知道\sigma^2的值，因此还需要对这个量进行估计。定义残差向量为\large{\vec{e}=\left[\begin{array}{c} {e}_1 \\ {e}_2 \\ \vdots \\ {e}_n\end{array}\right]=\left[\begin{array}{c} y_1-\hat{y}_1 \\ y_2-\hat{y}_2 \\ \vdots \\ y_n-\hat{y}_n\end{array}\right]=(I-H)\vec{y}}根据第一部分中H的特征值分解形式，并且V为正交方阵，可以得到\large{I-H=V(I-\Lambda)V^T=\left[\begin{array}{c} V_1 & V_2 \end{array}\right]\begin{bmatrix}O & O \\ O & I_{n-k-1} \end{bmatrix}\left[\begin{array}{c} V_1^T \\ V_2^T \end{array}\right]=V_2V_2^T}利用公式(I-H)X\vec{\beta}=X\vec{\beta}-X(X^TX)^{-1}X^TX\vec{\beta}=X\vec{\beta}-X\vec{\beta}=O，总残差RSS可写为RSS=\sum_{i=1}^n{e_i^2}=\|(I-H)\vec{y}\|_2^2=\|(I-H)(X\vec{\beta}+\vec{\epsilon})\|_2^2=\|(I-H)\vec{\epsilon}\|_2^2=\vec{\epsilon}^TV_2\underbrace{V_2^TV_2}_{I_{n-k-1}}V_2^T\vec{\epsilon}=\|V_2^T\vec{\epsilon}\|_2^2令\widetilde{\vec{e}}=V_2^T\vec{\epsilon}，容易看出\large{\begin{cases}E(\widetilde{\vec{e}})=V_2^TE(\vec{\epsilon})=0 \\ Cov(\widetilde{\vec{e}})=V_2^TCov(\vec{\epsilon})V_2=\sigma^2I_{n-k-1}\end{cases}}也就是说\widetilde{\vec{e}}为n-k-1个相互独立的分布为N(0,\sigma^2)的随机变量。RSS的均值为\large{E(RSS)=E(\|\widetilde{\vec{e}}\|_2^2)=\sum_{i=1}^{n-k-1}E(\widetilde{e}_i^2)=(n-k-1)\sigma^2}因此方差\sigma^2的无偏估计值可以写为\large{s^2=\frac{RSS}{n-k-1}=\frac{1}{n-k-1}\sum_{i=1}^n{(y_i-\hat{y}_i)^2}}并且满足分布\large{\frac{(n-k-1)s^2}{\sigma^2}\sim{\chi^2_{n-k-1}}}
将\hat{\vec{\beta}}改写成如下形式：\hat{\vec{\beta}}=(X^TX)^{-1}X^T\vec{y}=\vec{\beta}+(X^TX)^{-1}X^T\vec{\epsilon}=\vec{\beta}+(X^TX)^{-1}X^TX(X^TX)^{-1}X^T\vec{\epsilon}=\vec{\beta}+(X^TX)^{-1}X^TH\vec{\epsilon}根据公式\large{Cov(H\vec{\epsilon}, V_2^T\vec{\epsilon})=E(H\vec{\epsilon}\vec{\epsilon}^TV_2)=HE(\vec{\epsilon}\vec{\epsilon}^T)V_2=\sigma^2HV_2=\sigma^2V_1\underbrace{V_1^TV_2}_O=O}容易看出参数估计值\hat{\vec{\beta}}与s^2是相互独立的。