回归分析第3讲——线性回归的假设检验

一、模型检验

一旦我们完成了线性回归的参数拟合之后，就要对模型进行检验。那么为什么要进行检验？目的无非有两个：

1. 检测模型的显著性，即模型结构是否存在，在数学上也就是验证回归系数$\beta_1,\cdots,\beta_p$之中至少有一个系数不为0，这样模型才会有线性成分，模型结构才会成立；
2. 检验回归系数的显著性，即到底哪些回归系数可以认为不是0的，即真正对观测值$y$有贡献的是哪些项，这样我们才能精准地找到关键特征。

那么接下来我们就从一般线性回归开始，介绍模型检验的原理和方法。

1.1 一般线性回归的检验

为了不失一般性，并且将模型显著性和系数显著性检验统一起来，我们设模型为：

$$Y=X_{n\times p}\beta+e,e\sim N(0,\sigma^2I_n)$$

假定$rank(X)=r$，给定已知矩阵$H_{m\times p}$，$rank(H)=m$，那么我们的假设为：

• 原假设/零假设$H_0$：$H\beta=0$
• 备择假设$H_1$：$H\beta\ne0$

要完成对假设的检验，需要构造一些统计量。我们借助似然函数，构造两个假设的似然比：

$$\Lambda=\frac{\sup_{H_1}L(X,Y,\beta,\sigma^2)}{\sup_{H_0}L(X,Y,\beta,\sigma^2)}$$

如果模型不显著，似然比应该为1。如果似然比$\Lambda$越偏离1，说明假设$H_1$是更可能的，达到一定阈值我们就可以接受$H_1$，认为模型显著了。

对于分子来说，当$\hat\beta=(X^TX)^-X^TY$，$\sigma^2=\frac{||Y-X\hat\beta||^2}{n}$时取得最大值，即：

$$\sup_{H_1}L(X,Y,\beta,\sigma^2)=\left(\sqrt{\frac{2\pi e}{n}}||Y-X\hat\beta||\right)^{-n}$$

对于分母来说，当$\hat\beta$是$H\beta=0$约束下的最小二乘估计$\hat\beta_H$，$\sigma^2_H=\frac{||Y-X\hat\beta_H||^2}{n}$时取得最大值，即：

$$\sup_{H_0}L(X,Y,\beta,\sigma^2)=\left(\sqrt{\frac{2\pi e}{n}}||Y-X\hat\beta_H||\right)^{-n}$$

因此似然比可以写为：

$$\Lambda=\left(\frac{||Y-X\hat\beta_H||^2}{||Y-X\hat\beta||^2}\right)^{\frac{n}{2}}$$

接下来，我们想要将似然比和一个统计量关联起来。似然比中出现残差的平方，并且我们认为残差服从正态分布（联想到卡方分布），那么是不是可以让分子分母分别除以各自的自由度，就能得到一个$F$统计量呢？没错，线性回归中著名的$F$检验就是这样来的，因此我们的似然比可以这样转化为：

$$\Lambda=\left(\frac{||Y-X\hat\beta_H||^2}{||Y-X\hat\beta||^2}\right)^{\frac{n}{2}}=\left( 1+\frac{k}{n-r}F \right)^\frac{n}{2}$$

其中：

$$F=\frac{(ESS_H-ESS)/k}{ESS/(n-r)}\sim F_{k,n-r}$$

并且$ESS= ||Y-X\hat\beta||^2$，$ESS_H= ||Y-X\hat\beta_H||^2$，$k=r+m-rank\begin{pmatrix}X\\ H\end{pmatrix}$

只要计算了似然比，就能根据公式计算$F$，之后就能进行$F$检验了。

你可能注意到上面的表述中出现了残差(residual)一词。残差是指观测值和拟合值之间的差距，而之前我们提到的误差$e$则是观测值和真实值之间的差距。我们在构建模型时使用的是“真实值=线性预测子+误差”，但是在用实际数据拟合时，“真实值”是未知的，只能通过观测值去描述；“线性预测子”是未知的，只能通过拟合值去描述；那么观测值和拟合值之间的差距就是残差。

除了通过代数层面去理解$F$统计量，我们还可以从几何层面去理解。我们已经知道，一般线性回归的回归系数的最小二乘解为：

$$\beta=(X^TX)^{-}X^TY$$

也就是说，我们通过拟合得到的预测值为：

$$\hat Y=X\hat\beta=X(X^TX)^-X^TY=P_XY$$

因此，线性回归拟合的本质就是把响应变量正交投影到由设计矩阵$X$的列向量所张成的空间$\mu(X)$当中。有了这个概念之后，那么就很容易理解$F$检验了，如图所示：

对于同一组响应变量$Y$，备择假设就是把$Y$投影到Large model space中，而原假设就是把$Y$投影到Small model space中，如图中两条虚线箭头所示。根据向量的运算规则，我们可以很容易的标注出两种假设对应的残差，以及两个模型的差异。可以想象，如果原假设和备择假设之间差异很小，那么图中的点状虚线的长度应当很短，反之则很长。因此，我们所构建的$F$统计量就是看这个差值究竟是不是显著的大。

进一步地，如果$\mu(H^T)\subset \mu(X^T)$，则说明$H\beta$的每个分量都是可估的，此时$rank(X)=rank\begin{pmatrix}X\\ H\end{pmatrix}=r$，那么$k=m$。此时$\hat\beta_H$有更简单的计算方法：

$$\hat\beta_H=\hat\beta-(X^TX)^-H^T[H(X^TX)^-H^T]^{-1}H\hat\beta$$

考虑完一般线性回归，接着我们把结论迁移到有约束线性回归模型当中。这里省略推导过程，如果原模型就已经有约束$L_{q\times p}\beta=0$，且$rank(L_{q\times p})=q$，那么$F$检验统计量转换为：

$$F=\frac{(ESS_{LH}-ESS_L)/(m_1-m_2)}{ESS_L/(n-m_1)}$$

其中：

$$\begin{aligned} ESS_L&= ||Y-X\hat\beta_L||^2\\ ESS_{LH}&= ||Y-X\hat\beta_{LH}||^2\\ m_1&=rank\begin{pmatrix}X\\ H\end{pmatrix}-rank(L)\\ m_2&=rank\begin{pmatrix}X\\ L\\ H \end{pmatrix}-rank\begin{pmatrix}L\\ H \end{pmatrix} \end{aligned}$$

最后，我们考虑带截距项的线性回归模型。我们对带截距项的线性回归模型进行中心化，并假定所有回归系数都是可检验的，那么其检验的假设是$H\begin{pmatrix}\alpha\\ \beta_I\end{pmatrix}=0$，其中$H=\begin{pmatrix}0\vdots I_{p-1}\end{pmatrix}$。仿造前面的构造方法，我们可以构造统计量：

$$F=\frac{\hat\beta_I^T\widetilde X_c^TY/(p-1)}{(Y^TY-n\overline y^2-\hat\beta_I^T\widetilde X_c^TY)/(n-p)}\sim F_{p-1,n-p}$$

如果拒绝了原假设，那么就需要逐个对$\beta_i$进行检验，即原假设变为$H_{0i}:\beta_i=0$。如果误差服从正态分布，那么有：

$$\hat\beta_i\sim N(\beta_i,\sigma^2c_{ii}),C=(c_{ii})_{p\times p}=(X^TX)^{-1}$$

此时构造$t$统计量求解：

$$t_i=\frac{\hat\beta_i}{\hat\sigma\sqrt{c_{ii}}}\sim t_{n-p}$$

注意，如果误差不服从正态分布，那么上述构造不会精确地服从$t$分布（因为$\hat\beta_i$不服从正态分布），但是在大样本情况下上述结论依然可以使用。

如果接受了$H_{0i}$，就要考虑剔除$\beta_i$。

1.2 线性回归检验的推广

1.1小节考虑的检验较为特殊（但已经满足日常应用需求），为了完整性，这一小节我们进行简单的推广。

之前我们都在检验$H\beta=0$，那如果原假设是$H\beta=d$该如何处理？

我们先考虑$H\beta$的每个分量都可估，即$\mu(H^T)\subset\mu(X^T)$，此时对于线性回归$Y=X\beta+e,e\sim N(0,\sigma^2I_n)$，取$H\beta=d$的一个特解$\beta_0$，令$Z=Y-X\beta_0,\theta=\beta-\beta_0$，那么模型可以转换为：

$$Z=X\theta+e,e\sim N(0,\sigma^2I_n)$$

这样原假设就又重新变为$H\theta=0$了，此时对应的$F$统计量为：

$$F=\frac{(H\hat\beta-d)^T[H(X^TX)^-H^T]^{-1}(H\hat\beta-d)/m}{||Y-X\hat\beta||^2/(n-r)}\sim F_{m,n-r}$$

这样就能按照之前的方法检验了。

接着我们考虑$\mu(H^T)\not\subset\mu(X^T)$，也就是说假设当中含有不可估的成分，那么我们可以把$H$进行划分：

$$H=\begin{pmatrix}H_1\\H_2\end{pmatrix}$$

其中$H_1\beta$是可估的部分，$H_2\beta$是不可估的部分。那么此时存在两个原假设：

• $H_0: H\beta=0$
• $H_{01}:H_1\beta=0$

于是线性回归可以写为：

$$\begin{aligned} &Y=\theta+e,\theta\in\mu(X),e\sim N(0,\sigma^2I_n)\\ &H_0:\theta\in S=\left\{X\beta|H\beta=0,\beta\in \mathbb R^p\right\}\\ &H_{01}:\theta\in S_1=\left\{X\beta|H_1\beta=0,\beta\in \mathbb R^p\right\} \end{aligned}$$

显然有$S\subset S_1$，并且根据两个空间的维数：

$$\begin{aligned} &dim\ S_1=rank\begin{pmatrix}X\\H_1\end{pmatrix}-rank(H_1)=rank(X)-rank(H_1)\\ &dim\ S=rank\begin{pmatrix}X\\H\end{pmatrix}-rank(H)=rank\begin{pmatrix}X\\H_2\end{pmatrix}-rank\begin{pmatrix}H_1\\H_2\end{pmatrix}=rank(X)-rank(H_1) \end{aligned}$$

可以发现，其实$S=S_1$，也就是说$H_0$和$H_{01}$等价，因此对于不可估的假设是无法检验的，称为不可检验的假设(Non-testable hypothesis)。

二、参数的置信区间

前面我们讲解了如何对线性回归中的回归系数进行检验，那么接下来就需要对计算回归系数的置信区间，以把握系数的波动范围。

2.1 置信区间的构造

设$\Phi=H_{m\times p}\beta=(h_1^T\beta,\cdots,h_m^T\beta)^T$是$m$个独立可估函数，那么$rank(H)=m$，$\mu(H^T)=\mu(X^T)$。在没有初始约束下，由于$\hat\beta=(X^TX)^-X^TY$，那么：

$$\hat\Phi=H\hat\beta\sim N_m(\Phi,\sigma^2H(X^TX)^-H^T)$$

要构造置信区间，首先要构造一个枢轴量（pivot）服从一定的分布。由于$\hat\Phi$服从正态分布，那么标准化后再平方，一定服从卡方分布，即：

$$(\hat\Phi-\Phi)^T[H(X^TX)^-H^T]^{-1}(\hat\Phi-\Phi)/\sigma^2\sim \chi^2_m$$

由于方差的估计值为$\hat\sigma^2=||Y-X\hat\beta||^2/(n-r)$，且由于假设误差服从正态分布，所以根据之前的结论，方差估计值与$\hat\Phi$相互独立，并且：

$$(n-r)\hat\sigma^2/\sigma^2\sim\chi_{n-r}^2$$

所以两个式子分别除以自由度，再相除，就可以构造一个$F$统计量作为枢轴量：

$$(\hat\Phi-\Phi)^T[H(X^TX)^-H^T]^{-1}(\hat\Phi-\Phi)/m\hat\sigma^2\sim F_{m,n-r}$$

根据这个枢轴量，我们就可以构造置信区间：

$$D=\left\{ \Phi|(\hat\Phi-\Phi)^T[H(X^TX)^-H^T]^{-1}(\hat\Phi-\Phi)\le m\hat\sigma^2F_{m,n-r}(\alpha) \right\}$$

可以看出，$D$是以$\hat\Phi$为中心的一个置信椭球，且有$P(\Phi\in D)=1-\alpha$。

另外，如果只有一个可估函数，即$m=1$，我们就可以用更准确的$t$分布来构造，即：

$$\left| \frac{h^T\hat\beta-h^T\beta}{\hat\sigma\sqrt{h^T(X^TX)^-h}} \right|<t_{n-r}(\alpha/2)$$

于是有：

$$h^T\beta\in h^T\hat\beta\pm t_{n-r}(\alpha/2)\hat\sigma\sqrt{h^T(X^TX)^-h}$$

2.2 同时置信区间

不知道你发现没有，上面这种置信椭球的构造有点问题。如果我们对每一个可估函数都要计算一次置信椭球，那么犯Ⅰ类错误的概率会增大，假阳性的概率会提高。因此，对每一个可估函数应当同时给出区间估计，由此得到的置信区间称为同时置信区间。

这里我们主要介绍三种同时置信区间：

1. Bonferroni-t区间
2. 最大模-t区间
3. Scheffe置信带

我们首先介绍Bonferroni-t区间。Bonferroni提出了一个不等式，即对于一系列事件$A$：

$$P\left(\bigcup_iA_i\right)\le \sum_iP(A_i)$$

也就是说，并事件的概率不会超过事件独立时的概率之和。因此，对于$m$个可估函数和置信度$\alpha$，应当有：

$$P\left(\bigcap_{i=1}^mA_i\right)\ge 1-m\alpha$$

因此，为了取得名义上的$1-\alpha$，只需要用$\frac{\alpha}{m}$来代替原来的置信度$\alpha$即可，此时：

$$h_i^T\hat\beta\pm t_{n-r}(\alpha/2m)\hat\sigma\sqrt{h_i^T(X^TX)^-h_i}$$

这样得到的区间称为Bonferroni-t区间，这个区间的好处就是简单易行，且不需要任何额外的条件，但是缺点就是太严格，一刀切，假阳性的控制过强。

接下来就介绍最大模-t区间。最大模-t区间依赖的是多元非中心化$t$分布，如果我们令$V=(v_{ij})_{m\times m}=H(X^TX)^-H^T$，$x_i=(h_i^T\hat\beta-h_i^T\beta)/\sqrt{v_{ii}}$，因此：

$$X=(x_1,\cdots,x_m)^T\sim N_m(0,\sigma^2 R),R=(r_{ij})_{m\times n},r_{ij}=\frac{v_{ij}}{\sqrt{v_{ii}v_{jj}}}$$

那么对于$t_i=x_i/\hat\sigma$，$t\sim t_m(0,R,n-r)$，于是：

$$P\left\{\max_{1\le i\le m}|t_i|\le t_m^{\alpha/2}(0,R,n-r)\right\}=1-\alpha$$

所以构造$m$个同时置信区间为：

$$h_i^T\hat\beta\pm \hat\sigma\sqrt{v_{ii}}t_m^{\alpha/2}(0,R,n-r)$$

最大模-t区间放宽了限制，但是有个缺点，即$t_m^{\alpha/2}(0,R,n-r)$的计算非常困难，同时由于$R$是依赖于设计矩阵$X$的，每次$X$变化时都要重新计算，不利于查表。当然，如果$R=I_m$，那么还是比较容易计算的，不过计算得到的区间会被放大。

最后我们介绍scheffe区间。scheffe认为，对于任意可估的$l^T\beta\ (l\in\mu(H^T))$，其同时置信区间为：

$$l^T\hat\beta\pm[mF_{m,n-r}(\alpha)]^{\frac{1}{2}}\hat\sigma[l^T(X^TX)^-l]^{\frac{1}{2}},m=rank(H)$$

现证明上述置信区间的正确性。

要完成证明，首先给出引理，引理证明过程请查看附录。

引理：设矩阵$A_{n\times n}\gt0$，则$\sup_{b\ne0}\frac{(a^Tb)^2}{b^TAb}=a^TA^{-1}a$。

有了引理之后，由于$l\in\mu(H^T)$，说明存在$b$使得$l=H^Tb$，也就是$l^T=b^TH$；

令$a=H(\hat\beta-\beta)$，那么：

$$\begin{aligned}&P\left(l^T\hat\beta\le [mF_{m,n-r}(\alpha)]^{\frac{1}{2}}\hat\sigma[l^T(X^TX)^-l]^{\frac{1}{2}}\right)\\ &\Leftrightarrow P\left(\forall b,(a^Tb)^2\le mF_{m,n-r}(\alpha)\hat\sigma^2[b^TH(X^TX)^-Hb]\right)\\ &=P\left(\forall b\ne0,\frac{(a^Tb)^2}{b^TH(X^TX)^-H^Tb}\le mF_{m,n-r}(\alpha)\hat\sigma^2\right)\\ &=P\left(\sup_{b\ne0}\frac{(a^Tb)^2}{b^TH(X^TX)^-H^Tb}\le mF_{m,n-r}(\alpha)\hat\sigma^2\right)\\ &=P\left(a^T[H(X^TX)^-H^Tb]^{-1}a\le mF_{m,n-r}(\alpha)\hat\sigma^2\right)\\ &=P\{(H\hat\beta-H\beta)^T[H(X^TX)^-H^T]^{-1}(H\hat\beta-H\beta)\\ &\le mF_{m,n-r}(\alpha)\hat\sigma^2\}\\ &=1-\alpha\end{aligned}$$

原命题得证。

这里还需要补充额外的说明：

1. 对于有限的可估函数，scheffe置信带会偏长，只有当$m=r$时所有可估。
2. 如果$X$是列满秩的，即$rank(X_{n\times p})=p$，那么所有的$l^T\beta$可估。此时由于$l$的区间变化，所以会构成置信带。
3. 对于线性模型$y=\beta_0+\beta_1x$，其$1-\alpha$置信带为：

$$(\hat\beta_0+\hat\beta_1x)\pm [2F_{2,n-2}(\alpha)]^{\frac{1}{2}}\hat\sigma\sqrt{\frac{1}{n}+\frac{(x-\overline x)^2}{\sum_{i=1}^n(x_i-\overline x)^2}}$$

该置信带关于经验直线$\hat y=\hat\beta_0+\hat\beta_1x$对称，并且在$x=\overline x$处带宽最小。

三、线性回归的预测

前面我们对回归系数的显著性进行了检验，并讨论了系数的置信区间，最后我们要考虑在得到可信的回归系数后如何利用我们的模型进行预测。

3.1 点预测

对于线性模型$y_i=x_i^T\beta+e_i,i=1,\cdots,n$，$Ee=0$，$Cov(e)=\sigma^2\Sigma$，$\Sigma>0$且已知，$rank(X_{n\times p})=r$。现要预测$m$个点$x_i=(x_{i1},\cdots,x_{ip})^T,i=n+1,\cdots,n+m$的预测值$y_{n+1},\cdots,y_{n+m}$。

要完成预测，那么令$X_0=(X_{n+1}^T,\cdots,X_{n+m}^T)^T$，$Y_0=(y_{n+1},\cdots,y_{n+m})$，$e_0=(e_{n+1},\cdots,e_{n+m})$，于是所有的预测点也可以构成线性回归模型：

$$Y_0=X_0\beta+e_0,Ee_0=0,Cov(e_0)=\sigma^2\Sigma_0$$

假定$\mu(X_0^T)\subset\mu(X^T)$，分两种情况讨论：

• $Y_0$与$Y$不相关，即预测点与历史值无关，此时$Cov(e_0,e)=0$
• $Y_0$与$Y$相关，给定$V_{m\times n}$，设$Cov(e_0,e)=\sigma^2V$

对于不相关的情况，可以直接用$EY_0=X_0\beta$去预测，于是：

$$Y_0^*=X_0\beta^*=X_0(X^T\Sigma^{-1}X)^-X^T\Sigma^{-1}Y$$

设预测偏差$z=Y_0^*-Y_0$，由于$Ez=0$，所以这是无偏估计，那么：

$$Cov(z)=\sigma^2\left[ \Sigma_0+X_0(X^T\Sigma^{-1}X)^-X_0^T \right]$$

对于相关的情况，令$Y_0^*=C_{m\times n}Y$是$Y_0$的一个线性预测。要得到一个好的预测，我们给定评价指标为广义预测均方误差(Generalized Prediction Mean Squared Error,PMSE)，其计算方法为：

$$Give\ A\gt0,PMSE(Y_0^*)=E(Y_0^*-Y_0)^TA(Y_0^*-Y_0)$$

如果预测无偏且$PMSE$最小，则称为最优线性无偏预测(BLUP)。

对于$Cov(e_0,e)=\sigma^2V$，假定$X_0\beta$可估，那么$Y_0$的BLUP为：

$$Y_0^*=X_0\beta^*+V\Sigma^{-1}(Y-X\beta^*)$$

可以看到，如果$V=0$，那么BLUP转化为$Y_0^*=X_0\beta^*$。

上述$Y_0$的BLUP表达式的求解过程见附录。

3.2 区间预测

除了预测单独的点，我们还可以给出预测点所在的预测区间。为了简单起见，我们假设误差服从正态分布，即：

$$e\sim N_n(0,\sigma^2\Sigma),e_0\sim N_n(0,\sigma^2\Sigma_0)$$

并且认为$Y$与$Y_0$不相关，即$V=0$。因此，偏差$Z=Y_0^*-Y_0$服从正态分布，即：

$$Z=Y_0^*-Y_0\sim N(0,\sigma^2[\Sigma_0+X_0(X^T\Sigma^{-1}X)^-X_0^T])$$

同样假定$\mu(X_0^T)\subset\mu(X^T)$，$rank(X_{n\times p})=r$，并记：

$$\sigma^{*^2}=\frac{(Y-X\beta^*)^T\Sigma^{-1}(Y-X\beta^*)}{n-r},\Sigma_0=\left(\sigma_{ij}^{(0)}\right)_{n+1\le i,j\le n+m}$$

那么类似于2.1小节，对于第$i=n+1,\cdots,n+m$个点，$y_i$的$1-\alpha$预测区间为：

$$x_i^T\beta^*\pm t_{n-r}(\alpha/2)\sigma^*\left[\sigma_{ii}^{(0)}+x_i^T(X^T\Sigma^{-1}X)^-x_i\right]^{1/2}$$

进一步地，$y_{n+1},\cdots,y_{n+m}$的一个Bonferroni至少$1-\alpha$同时预测区间为：

$$x_i^T\beta^*\pm t_{n-r}(\alpha/2m)\sigma^*\left[\sigma_{ii}^{(0)}+x_i^T(X^T\Sigma^{-1}X)^-x_i\right]^{1/2}$$

同理，$y_{n+1},\cdots,y_{n+m}$的一个Scheffe至少$1-\alpha$同时预测区间为：

$$x_i^T\beta^*\pm \left[mF_{m,n-r}(\alpha)\right]^{1/2}\sigma^*\left[\sigma_{ii}^{(0)}+x_i^T(X^T\Sigma^{-1}X)^-x_i\right]^{1/2}$$

附录

2.2节引理：设矩阵$A_{n\times n}\gt0$，证明$\sup_{b\ne0}\frac{(a^Tb)^2}{b^TAb}=a^TA^{-1}a$。

证明：对于矩阵$A$，一定存在正交阵$Q$使得$A=Q^T\Lambda Q$，其中$\Lambda$是由特征值构成的对角阵，那么如果令$\widetilde a=Qa$，$\widetilde b=Qb$，则有：

$$\begin{aligned} \sup_{b\ne 0}\frac{(a^Tb)^2}{b^TAb}&=\sup_{b\ne0}\frac{(a^TQ^TQb)^2}{b^TQ^T\Lambda Qb}\\ &=\sup_{\widetilde b\ne0}\frac{(\widetilde a^T\widetilde b)^2}{\widetilde b^T\Lambda \widetilde b}\\ &=\sup_{\widetilde b\ne0}\frac{(\sum_{i=1}^n\widetilde{a_i}\widetilde{b_i})^2}{\sum_{i=1}^n\lambda_i\widetilde{b_i}^2}\\ &=\sup_{\widetilde b\ne0}\frac{\sum_{i=1}^n\widetilde{a_i}^2\frac{1}{\sqrt{\lambda_i}}\sqrt{\lambda_i}\widetilde{b_i}^2}{\sum_{i=1}^n\lambda_i\widetilde{b_i}^2}\\ &\le\sup_{\widetilde b\ne0}\frac{\sum_{i=1}^n\widetilde{a_i}^2\frac{1}{\lambda_i}\sum_{i=1}^n\lambda_i\widetilde{b_i}^2}{\sum_{i=1}^n\lambda_i\widetilde{b_i}^2}\\ &=\sum_{i=1}^n\widetilde{a_i}^2\frac{1}{\lambda_i}=\widetilde a^T\Lambda^{-1}\widetilde a\\ &=a^TQ^T\Lambda^{-1}Qa\\ &=a^TA^{-1}a \end{aligned}$$

证毕。

3.1节BLUP求解。当前的优化目标是：

$$\min E\left[(Y_0^*-Y_0)^TA(Y_0^*-Y_0)\right]$$

约束条件是：

$$Y_0^*=CY,E(Y_0^*-Y_0)=0$$

由约束条件可以得到：

$$CX\beta-X_0\beta=0\Rightarrow CX=X_0$$

又因为：

$$\begin{aligned}Var(Y_0^*-Y_0)&=Var(Y_0^*)+Var(Y_0)-2Cov(Y_0^*,Y_0)\\ &=\sigma^2\left[C\Sigma C^T+\Sigma_0-2CV^T\right]\end{aligned}$$

于是原命题转化为：

$$\min tr[A(Y_0^*-Y_0)]=\min tr[AC\Sigma C^T-2ACV^T]$$

构造拉格朗日乘子：

$$L(C,\Lambda)=tr[AC\Sigma C^T-2ACV^T]+2tr[\Lambda^T(CX-X_0)]$$

不难证明以下两个求导规则：

$$\begin{aligned}\frac{\partial tr(AC\Sigma C^T)}{\partial C}&=A^TC\Sigma^T+AC\Sigma\\ \frac{\partial tr(ACB)}{\partial C}&=A^TB^T\end{aligned}$$

于是：

$$\frac{\partial L(C,\Lambda)}{\partial C}=2AC\Sigma-2AV+2\Lambda X^T=0$$

再结合约束条件$CX=X_0$，求解得到：

$$\begin{aligned}\Lambda&=A(V\Sigma^{-1}X-X_0)(X^T\Sigma^{-1}X)^-\\ C&=V\Sigma^{-1}-(V\Sigma^{-1}X-X_0)(X^T\Sigma^{-1}X)^-X^T\Sigma^{-1}\end{aligned}$$

那么最终：

$$\begin{aligned}Y_0^*&=CY\\ &=V\Sigma^{-1}Y-(V\Sigma^{-1}X-X_0)(X^T\Sigma^{-1}X)^-X^T\Sigma^{-1}Y\\ &=V\Sigma^{-1}-(V\Sigma^{-1}X-X_0)\beta^*\\ &=X_0\beta^*+V\Sigma^{-1}(Y-X\beta^*)\end{aligned}$$

于是原命题得证。