目录
 
线性拟合
平面线性拟合
维度扩展
最小二乘参数估计
  
高斯噪声与最小二乘
探究高斯噪声与最小二乘的关系
参数估计
 

 
线性拟合
 
平面线性拟合
 
线性回归是最简单的数据拟合方法，平面空间的线性回归利于直观理解：
 
 平面上有很多样本点$(x,y)$，目标是找到一条直线$y=ax+b$拟合这些样本，在实际问题中，样本点的特征不局限于1维，而是任意的$p$维；
 
维度扩展
 
为了便于问题描述，引入符号：
 假设有一组样本
    
     
      
       
        D
       
       
        =
       
       
        
         {
        
        
         (
        
        
         
          x
         
         
          1
         
        
        
         ,
        
        
         
          y
         
         
          1
         
        
        
         )
        
        
         ,
        
        
         (
        
        
         
          x
         
         
          2
         
        
        
         ,
        
        
         
          y
         
         
          2
         
        
        
         )
        
        
         ,
        
        
         .
        
        
         .
        
        
         .
        
        
         ,
        
        
         (
        
        
         
          x
         
         
          n
         
        
        
         ,
        
        
         
          y
         
         
          n
         
        
        
         )
        
        
         }
        
       
      
      
       D=\left\{(x_{1},y_{1}),(x_{2},y_{2}),...,(x_{n},y_{n})\right\}
      
     
    D={(x1​,y1​),(x2​,y2​),...,(xn​,yn​)}，其中，$x_i$是一个$p$维向量，它可以表达第$i$个样本被观察的$p$个特征，$y_i$表示第$i$个样本的取值，是一个数值，所以拟合的直线表示样本$p$个特征到最终取值的线性映射关系；
 
因此，$N$个样本的集合写作：$X=[x_1,x_2,...,x_N{]}^T$；每个样本写作向量$x_i=[x_{i1},x_{i2},...,x_{ip}]$，最终目标是估计系数$w=[w_1,w_2,...,w_p]$和一个偏置$b$，从而建立映射：
 $$w^{T}x+b\Rightarrow y$$
 为了书写简便，可以将偏置$b$作为系数$w_0$，得到：
 
 即$w^{T}x+b\Rightarrow y$变成$w^{T}x\Rightarrow y$；
 
最小二乘参数估计
 
在最小二乘法估计中，定义目标函数如下：
 $$L(w)=\sum _{i=1}^N|w^T{x}_i-y_i{|}^2$$
 针对$N$个样本，找到一个系数向量$w$，使得拟合结果和真实值间误差的平方和最小；
 
高斯噪声与最小二乘
 
探究高斯噪声与最小二乘的关系
 
通过线性拟合，能否让直线精确通过每一个样本，使得拟合的误差为0；显然这是不可能的，因为样本本身携带噪声，带有随机性，所以可以得到另一种理解：拟合的直线代表样本分布的确定性，拟合值和真实值的误差（也可称为噪声），代表了随机性；
 
对于随机性的噪声$ϵ$，也许可以用高斯分布来描述：
 $$ϵ\sim N(0,{\sigma }^2)$$
 因此，可以得到：
 $$y=w^{T}x+ϵ$$
 显然，在$w$与$x$确定的情况下，$y$也将服从正态分布：
 $$y\sim N(w^{T}x,{\sigma }^2)$$
 其概率密度也是一个条件概率：
 $$p(y|x,w)=\frac{1}{\sqrt{2\pi {\sigma }^2}}exp(-\frac{(y-w^{T}x{)}^2}{2{\sigma }^2})$$
 此时，可以使用极大似然估计去获得参数$w$，对数似然函数为：
 $$L(w)=log\prod _{i=1}^{N}p(y_i|x_i,w)=\sum _{i=1}^{N}log(p(y_i|x_i,w))=\sum _{i=1}^N(log\frac{1}{\sqrt{2\pi {\sigma }^2}}-\frac{(y_i-w^T{x}_i{)}^2}{2{\sigma }^2})$$
 去除与$w$无关的项，目标化简为：
 $$w_{mle}=argma{x}_w\sum _{i=1}^N-(y_i-w^T{x}_i{)}^2=argmi{n}_w\sum _{i=1}^N(y_i-w^T{x}_i{)}^2$$
 可以看出，正好就是最小二乘法参数估计的目标函数$L(w)={\sum }_{i=1}^N|w^T{x}_i-y_i{|}^2$，所以，最小二乘法本身就隐含了噪声服从0均值正态分布的假设；
 
参数估计
 
现在需要求解$w_{mle}$，首先展开目标函数：
 
 对向量$[w^T{x}_1-y_1,w^T{x}_2-y_2,...,w^T{x}_N-y_N]$进行简单处理：
 $$[w^T{x}_1-y_1,w^T{x}_2-y_2,...,w^T{x}_N-y_N]=w^T[x_1,x_2,...,x_N]-[y_1,y_2,...,y_N]=w^T{X}^T-Y^T$$
 同样的，后一项作为转置有：
 $$(w^T{X}^T-Y^T{)}^T=Xw-Y$$
 得到：
 $$L(w)=(w^T{X}^T-Y^T)(Xw-Y)=w^T{X}^{T}Xw-2w^T{X}^{T}Y+Y^{T}Y$$
 计算关于$w$的偏导数，并令偏导数为0：
 $$\frac{\partial L(w)}{\partial w}=2X^{T}Xw-2X^{T}Y=0\Rightarrow w=(X^{T}X{)}^{-1}{X}^{T}Y$$

最小二乘估计的使用前提总是假设线性回归模型的误差是等方差且不相关的，即$Cov(e)={\sigma }^{2}I$，虽然在许多情况下，这个假定可以认为近似地成立，但有时我们的确要考虑假定不成立时的情况。
 为了讨论的简单，我们假定以下的的$\Sigma$（正常情况下是有参数的）是完全已知的。
 我们讨论的模型：$$\{\begin{array}{l}y=X\beta +e,(\ast )\\ E(e)=0,\\ Cov(e)={\sigma }^2\Sigma \end{array}$$
 (注：$\Sigma$是正定矩阵，故存在$P_{n\times n}$使得$\Sigma =P^{\prime }\Lambda P$,并记$({\Sigma }^{-\frac{1}{2}}{)}^2={\Sigma }^{-1}$)
 用${\Sigma }^{-\frac{1}{2}}$左乘$(\ast )$式，得$${\Sigma }^{-\frac{1}{2}}y={\Sigma }^{-\frac{1}{2}}X\beta +{\Sigma }^{-\frac{1}{2}}e$$记$Z={\Sigma }^{-\frac{1}{2}}y,U={\Sigma }^{-\frac{1}{2}}X\beta ,\epsilon ={\Sigma }^{-\frac{1}{2}}e$,即得到了以个满足基本假定的新模型(我们可以计算新模型期望和协方差阵，发现确实满高斯-马尔可夫假定计算协方差时用到公式$Cov(AX)=ACov(X)A^{\prime }$): 
    
     
      
       
      
      
       \quad
      
     
     $Z=U+\epsilon$, 因而我们可以得到新模型的最小二乘估计$${\beta }^{\ast }=(U^{\prime }U{)}^{-1}{U}^{\prime }Z=(X^{\prime }{\Sigma }^{-1}X{)}^{-1}{X}^{\prime }{\Sigma }^{-1}y$$一般地，我们就称${\beta }^{\ast }$为广义最小二乘估计。
 定理:
    
     
      
       
      
      
       \quad
      
     
    对于线性回归模型$(\ast )$的广义最小二乘估计${\beta }^{\ast }$，有以下性质：
 
$E({\beta }^{\ast })=\beta$
$Cov({\beta }^{\ast })={\sigma }^2(X^{\prime }{\Sigma }^{-1}X{)}^{-1}$
对于任意$c_{n\times 1}$向量，$c^{\prime }{\beta }^{\ast }$为$c^{\prime }\beta$的唯一的最小方差无偏估计。（说明对于一般线性回归模型$(\ast ),$广义最小二乘估计总是优于普通最小二乘估计的）

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

 
约束最小二乘估计的求解
 
在线性模型$y=X\beta +ϵ$中，我们考虑的约束也是线性的。假设系数$\beta$满足
 $$H\beta =d,H\in {\mathbb{R}}^{k\times p},rank(H)=k$$
 
并且$d$属于$H$的列空间（或者称为像空间），$d\in C(H)$，也就是说这个约束方程有界。假设$C(H^{\prime })\subset C(X^{\prime })$，即$H\beta$是$k$个线性无关的可估函数。
 
下面我们尝试用Lagrange乘子法求解带约束的最小二乘：
 $$\begin{gathered} \underset{\beta }{\min }Q={\Vert e\Vert }^2=(y-X\beta {)}^{\prime }(y-X\beta )=y^{\prime }y-2y^{\prime }X\beta +{\beta }^{\prime }{X}^{\prime }X\beta \\ s.t.H\beta =d \end{gathered}$$
 
用$L$表示Lagrange函数，$2\lambda \in {\mathbb{R}}^k$表示Lagrange乘子，则
 $$L(\beta ,\lambda )=y^{\prime }y-2y^{\prime }X\beta +{\beta }^{\prime }{X}^{\prime }X\beta +2{\lambda }^{\prime }(H\beta -d)$$
 
这里用$2\lambda$只是为了约掉2这个数值，让下面的正则方程形式上美观一点。计算Lagrange函数关于$\beta$的梯度可以得到正则方程：
 $$\begin{gathered} {\nabla }_{\beta }L=2X^{\prime }X\beta -2(X^{\prime }y-H^{\prime }\lambda )=0 \\ \Rightarrow X^{\prime }X\beta =X^{\prime }y-H^{\prime }\lambda \end{gathered}$$
 
数值计算的思路
 
记$\theta =[{\beta }^{\prime },{\lambda }^{\prime }{]}^{\prime }$, 约束方程可以写成
 $$\left[\begin{array}{cc}H& 0\end{array}\right]\theta =d$$
 
正则方程可以写成
 $$\left[\begin{array}{cc}{X}^{\prime }X& H^{\prime }\end{array}\right]\theta =X^{\prime }y$$
 
合并起来就是
 $$\left[\begin{array}{cc}{X}^{\prime }X& H^{\prime }\\ H& 0\end{array}\right]\theta =\left[\begin{array}{c}{X}^{\prime }y\\ d\end{array}\right]$$
 
求解$\theta$可以得到$\beta$与$\lambda$的估计值，
 $$\hat{\theta }={\left[\begin{array}{cc}{X}^{\prime }X& H^{\prime }\\ H& 0\end{array}\right]}^{-1}\left[\begin{array}{c}{X}^{\prime }y\\ d\end{array}\right]$$
 
系数估计量的解析式
 
数值上这样计算非常方便，但是我们想得到$\beta$估计量的解析式。考虑正则方程，
 $$\hat{\beta }=(X^{\prime }X{)}^{-1}(X^{\prime }y-H^{\prime }\hat{\lambda })={\hat{\beta }}_{OLS}-(X^{\prime }X{)}^{-1}{H}^{\prime }\hat{\lambda }$$
 
将这个结果代入约束方程中，
 $$\begin{gathered} H\hat{\beta }=H{\hat{\beta }}_{OLS}-H(X^{\prime }X{)}^{-1}{H}^{\prime }\hat{\lambda }=d \\ \Rightarrow \hat{\lambda }=[H(X^{\prime }X{)}^{-1}{H}^{\prime }{]}^{-1}(H{\hat{\beta }}_{OLS}-d) \end{gathered}$$
 
前面我们假设了$C(H^{\prime })\subset C(X^{\prime })$，并且$rank(H)=k$，因此$H(X^{\prime }X{)}^{-1}{H}^{\prime }$的逆与广义逆选取无关，这保证$\hat{\lambda }$形式的唯一性。由此我们得到系数的估计量为
 $$\hat{\beta }={\hat{\beta }}_{OLS}-(X^{\prime }X{)}^{-1}{H}^{\prime }[H(X^{\prime }X{)}^{-1}{H}^{\prime }{]}^{-1}(H{\hat{\beta }}_{OLS}-d)$$
 
约束最小二乘估计的统计性质
 
在约束参数空间
    
     
      
       
        {
       
       
        (
       
       
        β
       
       
        ,
       
       
        
         σ
        
        
         2
        
       
       
        )
       
       
        :
       
       
        H
       
       
        β
       
       
        =
       
       
        d
       
       
        }
       
      
      
       \{(\beta,\sigma^2):H\beta=d\}
      
     
    {(β,σ2):Hβ=d}中，${\hat{\sigma }}^2$是$\sigma$的无偏估计，其中
 $${\hat{\sigma }}^2=\frac{{\hat{e}}^{\prime }\hat{e}}{n-rank(X)+rank(H)},\hat{e}=y-X\hat{\beta }$$
 
与普通最小二乘法不同的是，约束最小二乘法的残差有更多自由度。普通最小二乘法总自由度为$n-1$，回归自由度（系数的自由度）为$rank(X)-1$；约束最小二乘法总自由度为$n+rank(H)-1$，回归自由度与普通最小二乘一样，所以多出来的自由度属于残差。
 
证明
 考虑${\hat{e}}^{\prime }\hat{e}={\Vert y-X\hat{\beta }\Vert }^2={\hat{e}}^{\prime }\hat{e}={\Vert y-X({\hat{\beta }}_{OLS}+\hat{\beta }-{\hat{\beta }}_{OLS})\Vert }^2$，进一步化简得到
 $${\Vert (y-X{\hat{\beta }}_{OLS})+X(\hat{\beta }-{\hat{\beta }}_{OLS})\Vert }^2$$
 
注意到$y-X{\hat{\beta }}_{OLS}$与$C(X^{\prime })$正交，因此上式等于
 $${\Vert y-X{\hat{\beta }}_{OLS}\Vert }^2+{\Vert X(\hat{\beta }-{\hat{\beta }}_{OLS})\Vert }^2$$
 
上一讲证明了
 $$E{\Vert y-X{\hat{\beta }}_{OLS}\Vert }^2=(n-rank(X)){\sigma }^2$$
 
并且证明了一个恒等式：如果$EX=\mu ,Cov(X)=\Sigma$，则
 $$E[X^{\prime }AX]={\mu }^{\prime }A\mu +tr(A\Sigma )$$
 
接下来我们基于这个恒等式计算$E{\Vert X(\hat{\beta }-{\hat{\beta }}_{OLS})\Vert }^2$,
 $$\begin{gathered} E{\Vert X(\hat{\beta }-{\hat{\beta }}_{OLS})\Vert }^2 \\ =E(H{\hat{\beta }}_{OLS}-d{)}^{\prime }[H(X^{\prime }X{)}^{-1}{H}^{\prime }{]}^{-1}(H{\hat{\beta }}_{OLS}-d) \\ =(H\beta -d{)}^{\prime }[H(X^{\prime }X{)}^{-1}{H}^{\prime }{]}^{-1}(H\beta -d) \\ +tr[[H(X^{\prime }X{)}^{-1}{H}^{\prime }{]}^{-1}Cov(H{\hat{\beta }}_{OLS})] \end{gathered}$$
 
在参数空间
    
     
      
       
        {
       
       
        (
       
       
        β
       
       
        ,
       
       
        
         σ
        
        
         2
        
       
       
        )
       
       
        :
       
       
        H
       
       
        β
       
       
        =
       
       
        d
       
       
        }
       
      
      
       \{(\beta,\sigma^2):H\beta=d\}
      
     
    {(β,σ2):Hβ=d}中，第一项$$(H\beta -d{)}^{\prime }[H(X^{\prime }X{)}^{-1}{H}^{\prime }{]}^{-1}(H\beta -d)=0$$
 
计算第二项，根据上一讲的最后一个定理，
 $$Cov(H{\hat{\beta }}_{OLS})={\sigma }^2{H}^{\prime }(X^{\prime }X{)}^{-1}H$$
 
因此
 $$\begin{gathered} [H(X^{\prime }X{)}^{-1}{H}^{\prime }{]}^{-1}Cov(H{\hat{\beta }}_{OLS})={\sigma }^2{I}_k \\ \Rightarrow tr[[H(X^{\prime }X{)}^{-1}{H}^{\prime }{]}^{-1}Cov(H{\hat{\beta }}_{OLS})]=tr({\sigma }^2{I}_k)=k{\sigma }^2 \end{gathered}$$
 
这里$k=rank(H)$，所以
 $$E{\Vert y-X\hat{\beta }\Vert }^2=(n-rank(X)+rank(H)){\sigma }^2$$
 
证毕

最小二乘估计矩阵形式的推导
 
最近写文章有用到一些算法，自己推一下，顺便mark下来。
 这么久我才发现csdn居然都能写Tex了(666)。
 
考虑一般线性回归模型(OLR)
 
考虑只含有一个指标的一般线性回归模型(ordinary linear regression model)有如下形式：
 $$y_i={\beta }_0+{\beta }_1{x}_{i1}+ϵ，i=1,2,\dots ,n$$
 显然这是基于$n$个观测数据或者叫样本的模型形式。其中${\beta }_0$称为截距项系数，${\beta }_1$称为$x_1$的回归系数，它们都是未知的常值参数。$ϵ$是不能被观测到的随机误差项，并且满足$E(ϵ)=0$,$\mathrm{V}\mathrm{a}\mathrm{r}(ϵ)={\sigma }^2>0$。其实是有$x_0$的,只是通常认为$x_0=1$。还有一个关键的假设就是$x$不是随机变量($x$要都随机了,这模型就没法玩了)。
 
实际上我们所研究的问题往往包含多个指标。那么这些指标$(x_1,x_2,...,x_p)$就对对应$({\beta }_0,{\beta }_1,...,{\beta }_p)$个回归系数，这个时候模型的形式就变成了多元线性回归模型：
 $$y_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_p{x}_{ip}+{ϵ}_i，i=1,2,\dots ,n$$
 所以为了简化计算和书写方便，我们可以把它写成矩阵的形式：
 $$Y=X\mathbf{\beta }+\mathbf{ϵ}$$
 
$$Y=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_n\end{array}\right]X=\left[\begin{array}{cccc}1& x_{11}& \cdots & x_{1p}\\ 1& x_{21}& \cdots & x_{2p}\\ ⋮& ⋮& ⋮& ⋮\\ 1& x_{n1}& \cdots & x_{np}\end{array}\right]\mathbf{\beta }=\left[\begin{array}{c}{\beta }_0\\ {\beta }_1\\ ⋮\\ {\beta }_p\end{array}\right]\mathbf{\epsilon }=\left[\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ ⋮\\ {\epsilon }_n\end{array}\right]$$
 其中$X$称为设计矩阵(只是习惯叫法)，$Y$就不多说了。同样也有一些前提：$X$必须是列满秩；随机误差向量$\mathbf{\epsilon }$要满足高斯-马尔科夫条件(1829年，高斯提供了最小二乘法的优化效果强于其他方法的证明，因此被称为高斯-马尔可夫定理)：
 (i) $E(\mathbf{\epsilon })=0$
 (ii)$\mathrm{V}\mathrm{a}\mathrm{r}(\mathbf{\epsilon })={\sigma }^2\mathbf{I}$
 
最小二乘估计
 
最小二乘估计法$(LSE)$，它和机器学习领域的梯度下降法还是有一定的区别的(后者没有这么多假设，实用性更广泛)，准确的来讲$LES$只是一种算法，因为随机误差向量$\mathbf{ϵ}$并不能被观测，所以回归方程不存在解，我们只能尽可能的去接近真实值从而解出全局最优解，即确定一个$\hat{\mathbf{\beta }}$使得$\mathbf{\epsilon }=Y-X\mathbf{\beta }$各元素的平方和达到最小，可以记为：
 $$\begin{array}{rl}Q(\mathbf{\beta })& =\sum _{i=1}^n{\epsilon }_i^2\\ & ={\mathbf{\epsilon }}^T\mathbf{\epsilon }\\ & =(Y-X\mathbf{\beta }{)}^T(Y-X\mathbf{\beta })\\ & =(Y^{T}Y-2{\mathbf{\beta }}^T{X}^{T}Y+{\mathbf{\beta }}^T{X}^{T}X\mathbf{\beta })\end{array}$$
 
令：
 $$\frac{\partial Q(\mathbf{\beta })}{\partial \beta }=-2X^{T}Y+2X^{T}X\mathbf{\beta }=0$$
 这里需要一些矩阵求导的概念，接下来我们就可以得到一个叫做正规方程 的东西：
 $$X^{T}X\mathbf{\beta }=X^{T}Y$$
 由$\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(X^{T}X)=\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(X)=p+1$知$X^{T}X$是正定矩阵，所以$X^X$存在逆矩阵，那么正规方法就有唯一解了：
 $$\hat{\mathbf{\beta }}=({\hat{\beta }}_0,{\hat{\beta }}_1,\cdots ,{\hat{\beta }}_p{)}^T=(X^{T}X{)}^{-1}{X}^{T}Y$$
 此时$\mathbf{\beta }$的估计就得到了，如果再把它带回到模型中去就有：
 $$\hat{Y}=X\hat{\mathbf{\beta }}=X(X^{T}X{)}^{-1}{X}^{T}Y=SY$$
 一般统计学上称$S$是$Y$的帽子矩阵，这个称呼是因为有$S$的存在使$Y$带上了帽子（总感觉怪怪的？）接下来看残差：
 $$\hat{\mathbf{\epsilon }}=Y-\hat{Y}=(I-H)Y$$
 $I$是$n$阶的单位矩阵，显然残差的总和为0，是因为$Q(\mathbf{\beta })$对截距项求偏导数等于0时：
 $$-2\sum _{i=1}^n[y_i-({\beta }_0+\sum _{i=1}^p{\beta }_i{x}_i)]=0$$
 这个式子很明显表达了当存在截距项时，残差和必然为0，这也是为什么200年前拉普拉斯放弃了最小一乘法。也可以证明最小二乘法得到的估计和最大似然估计的结果是相同的，都是无偏估计。关于最小二乘法的BLUE性质不是本文的重点不再赘述。
 
补充几个推导过程中用到的矩阵求偏导法则
 
$$\frac{\partial x^{T}a}{\partial x}=\frac{\partial a^{T}x}{\partial x}=a$$
 $$\frac{\partial x^{T}Ax}{\partial x}=Ax+A^{T}x$$
 如果$A$是对称的:$Ax+A^{T}x=2Ax$.
 至此推导过程完毕。
 
参考文献：梅长林,王宁《近代回归分析方法》[M],科学出版社，2012.