UA STAT687 线性模型II 最小二乘理论2 约束最小二乘估计
约束最小二乘估计的求解
数值计算的思路
系数估计量的解析式
约束最小二乘估计的统计性质

约束最小二乘估计的求解
在线性模型$y = X\beta+\epsilon$中，我们考虑的约束也是线性的。假设系数$\beta$满足

$H\beta = d, H \in \mathbb{R}^{k \times p},\ rank(H)=k$
并且$d$属于$H$的列空间（或者称为像空间），$d\in C(H)$，也就是说这个约束方程有界。假设$C(H') \subset C(X')$，即$H\beta$是$k$个线性无关的可估函数。
下面我们尝试用Lagrange乘子法求解带约束的最小二乘：

$\min_{\beta}\ \ Q = \left\| e \right\|^2 = (y-X\beta)'(y-X\beta)=y'y-2y'X\beta+\beta'X'X\beta \\ s.t.\ \ H\beta = d$
用$L$表示Lagrange函数，$2\lambda \in \mathbb{R}^k$表示Lagrange乘子，则

$L(\beta,\lambda) = y'y-2y'X\beta+\beta'X'X\beta+2\lambda'(H\beta-d)$
这里用$2\lambda$只是为了约掉2这个数值，让下面的正则方程形式上美观一点。计算Lagrange函数关于$\beta$的梯度可以得到正则方程：

$\nabla_{\beta} L = 2X'X\beta - 2(X'y-H'\lambda)=0 \\ \Rightarrow X'X\beta = X'y-H'\lambda$
数值计算的思路
记$\theta = [\beta',\ \lambda']'$, 约束方程可以写成

$\left[ \begin{matrix} H & 0 \end{matrix} \right]\theta = d$
正则方程可以写成

$\left[ \begin{matrix} X'X   & H' \end{matrix} \right]\theta = X'y$
合并起来就是

$\left[ \begin{matrix} X'X   & H' \\ H & 0\end{matrix} \right]\theta = \left[ \begin{matrix} X'y \\ d \end{matrix} \right]$
求解$\theta$可以得到$\beta$与$\lambda$的估计值，

$\hat{\theta} = \left[ \begin{matrix} X'X   & H' \\ H & 0\end{matrix} \right]^{-1}\left[ \begin{matrix} X'y \\ d \end{matrix} \right]$
系数估计量的解析式
数值上这样计算非常方便，但是我们想得到$\beta$估计量的解析式。考虑正则方程，

$\hat\beta = (X'X)^{-1}(X'y-H'\hat\lambda) = \hat\beta_{OLS}-(X'X)^{-1}H'\hat\lambda$
将这个结果代入约束方程中，

$H\hat\beta = H\hat\beta_{OLS}-H(X'X)^{-1}H'\hat\lambda=d \\ \Rightarrow \hat\lambda = [H(X'X)^{-1}H']^{-1}(H\hat\beta_{OLS}-d)$
前面我们假设了$C(H')\subset C(X')$，并且$rank(H)=k$，因此$H(X'X)^{-1}H'$的逆与广义逆选取无关，这保证$\hat{\lambda}$形式的唯一性。由此我们得到系数的估计量为

$\hat\beta = \hat\beta_{OLS}-(X'X)^{-1}H'[H(X'X)^{-1}H']^{-1}(H\hat\beta_{OLS}-d)$
约束最小二乘估计的统计性质
在约束参数空间$\{(\beta,\sigma^2):H\beta=d\}$中，$\hat{\sigma}^2$是$\sigma$的无偏估计，其中

$\hat{\sigma}^2 = \frac{\hat{e}'\hat{e}}{n-rank(X)+rank(H)},\ \hat{e}=y-X\hat{\beta}$
与普通最小二乘法不同的是，约束最小二乘法的残差有更多自由度。普通最小二乘法总自由度为$n-1$，回归自由度（系数的自由度）为$rank(X)-1$；约束最小二乘法总自由度为$n+rank(H)-1$，回归自由度与普通最小二乘一样，所以多出来的自由度属于残差。
证明

考虑$\hat{e}'\hat{e} = \left\| y-X\hat{\beta}\right\|^2 = \hat{e}'\hat{e} = \left\| y-X(\hat{\beta}_{OLS}+\hat\beta-\hat{\beta}_{OLS})\right\|^2$，进一步化简得到

$\left\| (y-X\hat{\beta}_{OLS})+X(\hat\beta-\hat{\beta}_{OLS})\right\|^2$
注意到$y-X\hat{\beta}_{OLS}$与$C(X')$正交，因此上式等于

$\left\| y-X\hat{\beta}_{OLS}\right\|^2+\left\| X(\hat\beta-\hat{\beta}_{OLS})\right\|^2$
上一讲证明了

$E\left\| y-X\hat{\beta}_{OLS}\right\|^2=(n-rank(X))\sigma^2$
并且证明了一个恒等式：如果$EX=\mu,Cov(X)=\Sigma$，则

$E[X'AX]=\mu'A\mu+tr(A\Sigma)$
接下来我们基于这个恒等式计算$E\left\| X(\hat\beta-\hat{\beta}_{OLS})\right\|^2$,

$E\left\| X(\hat\beta-\hat{\beta}_{OLS})\right\|^2\\=E(H\hat\beta_{OLS}-d)'[H(X'X)^{-1}H']^{-1}(H\hat\beta_{OLS}-d) \\ = (H\beta-d)'[H(X'X)^{-1}H']^{-1}(H\beta-d) \\+tr[[H(X'X)^{-1}H']^{-1}Cov(H\hat\beta_{OLS})]$
在参数空间$\{(\beta,\sigma^2):H\beta=d\}$中，第一项$(H\beta-d)'[H(X'X)^{-1}H']^{-1}(H\beta-d)=0$
计算第二项，根据上一讲的最后一个定理，

$Cov(H\hat\beta_{OLS})=\sigma^2H'(X'X)^{-1}H$
因此

$[H(X'X)^{-1}H']^{-1}Cov(H\hat\beta_{OLS})=\sigma^2I_k \\ \Rightarrow tr[[H(X'X)^{-1}H']^{-1}Cov(H\hat\beta_{OLS})] = tr(\sigma^2I_k)=k\sigma^2$
这里$k=rank(H)$，所以

$E \left\| y-X\hat{\beta}\right\|^2 = (n-rank(X)+rank(H))\sigma^2$
证毕

设分别为的最小二乘估计值,
 于是的观测值

　　　　, ,
 (2.1)

其中为误差的估计值,
 称为残差或剩余。令为的估计值,
 则有

　　　　,
 (2.2)

　　　　, ,
 (2.3)

(2.3)式表示实际值与估计值的偏离程度。欲使估计值与实际值拟合的最好,
 则应使残差平方和

　　　　

达到最小, 为此, 我们可以应用微分求极值原理确定,
 即解下列方程组

　　　　,
 (2.4)

即

　　　　,
 (2.5)

整理并化简则得以下正规方程组:

　　　　,
 (2.6)

如果记(2.6)式的系数矩阵为,
 右端常数项矩阵记为, 则有

　　　　,
 (2.7)

　　　　,
 (2.8)

因此正规方程(2.6)的矩阵形式为

　　　　,
 (2.9)

或

　　　　,
 (2.10)

其中为正规方程中待定的未知实数向量,
 如果系数矩阵满秩, 则存在,
 此时有

　　　　,
 (2.11)

(2.11)式即为多元线性回归模型(1.2)式中参数的最小二乘估计。

　　正规方程组(2.6)亦可表达为下述另一种形式, 如果记

　　　　, ,

　　　　,

则由(2.6)式中第一等式可解出

　　　　,
 (2.12)

再将(2.12)代入到(2.6)其它各式中并经化简整理可得

　　　　,
 (2.13)

又由

　　　　, ,

　　　　, ,

如果记

　　　　, ,
 (2.14)

　　　　, ,
 (2.15)

则(2.13)式可以表示为

　　　　,
 (2.16)

(2.16)式称为正规方程组, 解此方程组可得,
 再代入到(2.12)式中则得, 于是得回归方程

　　　　,
 (2.17)

(2.17)式称为回归超平面方程。

　　如果记(2.16)式的系数矩阵为,
 右端常数项向量为, 则

　　　　,

　　　　,

且记,
 则正规方程组(2.16)的矩阵形式为

　　　　,
 (2.18)

解(2.18)得

　　　　,
 (2.19)

再代回到(2.12), 则得到。

　

　　以下是一对多线性回归分析的两个例子。

　

例2.1 某养猪场估算猪的毛重, 测得14头猪的体长(cm)、胸围(cm)与体重(kg)数据如表１,
 试建立与及的预测方程。

表2.1




序号



体长()



胸围()



体重()



 1



 41



49



28



 2



 45



58



39



 3



 51



62



41



 4



 52



71



44



 5



 59



62



43



 6



 62



74



50



 7



 69



71



51



 8



 72



74



57



 9



 78



79



63



10



 80



84



66



11



 90



85



70



12



 92



94



76



13



 98



91



80



14



103



95



81




　　经计算: , , , ,

　　　　　　,

　　　　　　, 

　　　　　　,

　　　　　　,

　　　　　　,

　　于是正规方程组为

　　　　　　,

　　解此方程组得

　　　　　　, ,

　　又

　　　　　　,

　　因此所求预测回归方程为

　　　　　　

　　回归方程中系数与的含义是体长每增加1cm,
 则猪体重毛重平均增加0.522kg, 胸围每增加1cm, 则猪体重毛重平均增加0.475kg。

　

例2.2 某地区二化螟的第一代成虫发生量与四个因素有关,
 这四个因素分别如下, 已知原始观测数据如表2.2, 试建立二化螟发生总量的回归方程。

　　:
 冬季积雪期限(单位为周),

　　:
 每年化雪日期(以2月1日为1),

　　:
 二月份平均气温(℃),

　　:
 三月份平均气温(℃),

　　:
 二化螟发生总量(头),

　　经计算:

　　　　,

　　　　,

表2.2





序号























 1



10



26



 0.2



3.6



 9



 2



12



26



-1.4



4.4



17



 3



14



40



-0.8



1.7



34



 4



16



32



 0.2



1.4



42



 5



19



51



-1.4



0.9



40



 6



16



33



 0.2



2.1



27



 7



 7



26



 2.7



2.7



 4



 8



 7



25



 1.0



4.0



27



 9



12



17



 2.2



3.7



13



10



11



24



-0.8



3.0



56



11



12



16



-0.5



4.9



15



12



 7



16



 2.0



4.1



 8



13



11



15



 1.1



4.7



20







154



347



 4.7



41.2



312







11.8462



26.6923



 0.3615



 3.1692



24




　　　　,

　　于是

　　　　,

　　又

　　　　

　　　　＝24 + 0.99742×11.8462 + 1.62581×26.6923 + 11.19263×0.3615 + 16.95291×3.1692 ＝ 136.98554,

　　因此所求二化螟发生总量的预测回归方程为

　　　　。




本系列转自：http://hutangao.blog.163.com/blog/static/4888314200982852442975/

已知变量X和Y为线性关系（这里XY均为nx1的列向量），为了得知X和Y到底具有怎样的线性关系（也即求解X的系数），如果这是一个工程问题，我们解决这一问题的方法就是对X和Y进行采样，获得很多组样本，然后就能求解出系数了，按照线代的理论，系数矩阵为nxn方阵，且秩为n时，方程具有唯一解，如果采样点过多，也即方程的数目多于未知数的数目，则方程组无解，这时只能求出一个近似解，以不同的目的获得的近似解是不同的，如果为了使方程左右两边的误差的平方和最小，而获得的近似解，就是最小二乘解(所谓“二乘”，就是“平方”的意思，最小二乘就是最小平方和)。这个问题的证明在研究生的矩阵分析引论数学课上学过，现在也忘光了，只记得结论表达式了。

好歹还记得原理：对于一个超定方程组 AX=Y，求X。

解：假设X可能的解为Z，那么把Z代入X可得到Y1=AZ，那么建立目标函数J=e^T * e，其中e=(Y - Y1);必定有一个Z能够使得J最小，这个Z就是X的最小二乘解。

e^T * e就是Y的误差的平方和，能使得J取的最小值的解就是最小二乘解，求X的过程就是令J的一阶导为0。这都是废话，严格的求解方法如下：

 

 

例子：假设变量y是n个变量xi的线性组合，求系数。

设方程为AX=y，也即



 

为了计算出系数a1、a2、···an的值，我们至少需要n次X、Y的采样值，形如：





 

这样就能求解出系数a1、a2、···an的值了，如果采样样本不止n个，而是多于n个，也不要紧，虽然会会造成方程组无解，但是可以求出最小二乘解。

把方程组写成矩阵形式为：

XA=Y     （式2）



至此，就求出了所有的系数ai

 

 

 

 

-----------------------------------------------------后记---------------------------------------------------------------------

 

又翻了电子版课本，把最小二乘的证明过程也贴上来：

证明过程在《矩阵分析引论》第30页。证明过程需要用到子空间的概念，这一概念的定义在第6页。






 

 

 

最后简单整理一下证明过程：



 

多元加权最小二乘估计

选择“分析”-“回归”-“线性”

选入自变量与因变量

选择“保存”，在“残差”中将“未标准化”选中

选择“转换”-“计算变量”

输入如下公式，点击“确定”

选择“分析”-“相关”-“双变量”



将如下变量选入框中，点击“斯皮尔曼”，点击“确定”

得到等级相关系数图表如下

选择“分析”-“回归”-“权重估计”，将以下变量选入框中，因为x2与ABSE的等级相关系数为0.721，大于X1与ABSE的等级相关系数为0.443，因此权重变量选择x2



得到结果如下



加权最小二乘的回归方程为y=-266.962+1.696x1+0.47x2