平差和统计的概念区别

平差认为,或者说测量学的书籍认为,测量获得的值有一定的“误差”,需要进行“改正”,所以:
真值=观测值+改正值
而统计中和物理学中认为观测值是真值加上一些误差(或者称“噪声”)引起的,所以:
观测值=真值+误差
在摄影测量和数据处理的相关教材中沿用此规定。本文以统计的概念入手,所以沿用此规定,希望读者注意区别。

改正数和误差互为相反数。

线性模型

在数学中,变量关系有两种基本类型:函数关系和相关关系(dependent relationship),此种关系没有密切到如上所述的确定关系。

假设因变量y和k个自变量x1,x2,,xkx_1,x_2,\cdots,x_k之间存在简单的线性关系:
y=β0+β1x1++βkxk+εy=\beta_0+\beta_1x_1+\cdots+\beta_kx_k+\varepsilon
其中ε\varepsilon是一个随机变量进一步假定对自变量的n组不同取值,得到因变量y的n次观测,则有:
Y=Xβ+εY=X\beta+\varepsilon

其中
Y=(y1y2yn),X=(1x11x1k1x21x2k1xn1xnk),β=(β1β2βn),ε=(ε1ε2εn),Y= \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\y_n \end{pmatrix}, X=\begin{pmatrix} 1 & x_{11} & \ldots & x_{1k}\\ 1 & x_{21} & \ldots & x_{2k}\\ \vdots & \vdots & & \vdots \\ 1 & x_{n1} & \ldots & x_{nk} \end{pmatrix}, \beta= \begin{pmatrix} \beta_1 \\ \beta_2 \\ \vdots \\\beta_n \end{pmatrix}, \varepsilon= \begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \vdots \\\varepsilon_n \end{pmatrix},
这里ε\varepsilon表示随机误差向量,满足E(ε)=0,cov(ε,ε)=ΣE(\varepsilon)=\mathbf 0,cov(\varepsilon,\varepsilon)=\Sigma
称上述模型为线性模型,记为(Y,Xβ,Σ)(Y,X\beta,\Sigma),其中,Y称为观测向量,X是k个自变量在n此观测中的取值,因为选择观测的量是可以控制的,这是试验设计问题,称X为设计矩阵(design matirx),β\beta 是未知的参数向量,一般假定Σ\Sigma是已知的,在许多问题中还假定n此观测相互独立,有公共方差,此时Σ=σ2In\Sigma=\sigma^2I_n,这里σ2\sigma^2是未知参数,称为误差方差。

Gauss-Markov条件:cov(ε,ε)=σ2In,E(ε)=0cov(\varepsilon,\varepsilon)=\sigma^2I_n,E(\varepsilon)=\mathbf{0}

有时还需假定正态条件:εN(0,σ2In)\varepsilon\sim N(\mathbf{0},\sigma^2I_n)
(注意到正态分布的性质Y=Xβ+εN(Xβ,σ2In)Y=X\beta+\varepsilon\sim N(X\beta,\sigma^2I_n)
为了对未知参数进行估计,总假定试验次数n不小于线性回归模型包含的未知参数个数,且设计矩阵X是列满秩的,即:
rank(X)=k+1rank(X)=k+1

使偏离平方和i=1n(yiy~i)2\sum_{i=1}^n(y_i-\tilde y_i)^2取最小值的β\beta称为它的最小二乘估计(Least Squares Estimate),简记为LS估计。这种求估计量的方法称为最小二乘法(Method of Least Squares),始于C.F.Gauss(1809),H后来A.A.Markov(1900)做了重要工作,奠定了这方面基础。

一元线性回归模型

y=β0+β1x+εy=\beta_0+\beta_1x+\varepsilon

假定E(ε)=0,var(ε)=σ2E(\varepsilon)=0,var(\varepsilon)=\sigma^2,该式称为一元线性回归模型(Simple Linear Regression Model)。如果加上正态条件,称其为一元正态线性回归模型。因为相信大家已有概率统计知识,故直接讲一般情况——多元线性回归模型。

多元线性回归模型

我们在这里讨论多元线性回归模型(Y,Xβ,σ2In)(Y,X\beta,\sigma^2I_n)

参数β\beta的估计

使总偏离平方和
Q(β)=i=1n[yi(β0+j=1kxijβj)]2=(YXβ)T(YXβ)=YXβ2\begin{aligned} Q(\beta)&=\sum_{i=1}^n[y_i-(\beta_0+\sum_{j=1}^kx_{ij}\beta_j)]^2 \\ &= (Y-X\beta)^T(Y-X\beta)=||Y-X\beta||^2 \end{aligned}
达到最小的β\beta记为β^\hat\beta即:
YXβ^2=minβYXβ2||Y-X\hat\beta||^2=\underset{\beta}{min}||Y-X\beta||^2
(2)
β^\hat\betaβ\beta的最小二乘估计。

Q=Q(β)=(YXβ)T(YXβ)=YTY2βTXTY+βTXTXβ\begin{aligned} Q=&Q(\beta)=(Y-X\beta)^T(Y-X\beta)\\ =& Y^TY-2\beta^TX^TY+\beta^TX^TX\beta \end{aligned}
β\beta求导,并令其为0
Qβ=2XTY+2XTXβ=0\frac{\partial Q}{\partial \beta}=-2X^TY+2X^TX\beta=0
整理后即有:
XTXβ=XTYX^TX\beta=X^TY
称该式为正规方程组(system of normal equations),记N为正规方程组的系数阵
N=XTXN=X^TX
因为N为列满秩矩阵,可逆,则:
β^=(XTX)1XTY\hat\beta=(X^TX)^{-1}X^TY
(1)
是正规方程组的解。

因为存在偏导数为0的点不是极值点的情况,下面的定理告诉我们正规方程组的解就是β\beta的最小二乘估计

定理

  1. 正规方程组的解(1)必是β\beta的最小二乘估计
  2. β\beta的最小二乘估计必为正规方程组的解。

proof.proof.

  1. β~\tilde\beta是正规方程组的解。即β~\tilde\beta满足:
    (XX)β~=XY(X'X)\tilde\beta=X'Y
    β,\forall \beta,
    Q(β)=YXβ2=(YXβ)(YXβ)=[(YXβ~)+X(β~β)][(YXβ~)+X(β~β)]=(YXβ~)(YXβ~)+(β~β)XX(β~β)+2(ββ~)X(YXβ~)\begin{aligned} Q(\beta)=& ||Y-X\beta||^2=(Y-X\beta)'(Y-X\beta)\\ =& [(Y-X\tilde\beta)+X(\tilde\beta-\beta)]'[(Y-X\tilde\beta)+X(\tilde\beta-\beta)]\\ =& (Y-X\tilde\beta)' (Y-X\tilde\beta)+(\tilde\beta-\beta)'X'X(\tilde\beta-\beta)+2(\beta-\tilde\beta)'X'(Y-X\tilde\beta) \end{aligned}
    注意到:
    (ββ~)X(YXβ~)=(ββ~)(XYXXβ~)=0(\beta-\tilde\beta)'X'(Y-X\tilde\beta)=(\beta-\tilde\beta)'(X'Y-X'X\tilde\beta)=0
    (β~β)XX(β~β)=X(β~β)20(\tilde\beta-\beta)'X'X(\tilde\beta-\beta)=||X(\tilde\beta-\beta)||^2\ge 0
    Q(β)(YXβ~)(YXβ~)=Q(β~)\therefore Q(\beta)\ge (Y-X\tilde\beta)' (Y-X\tilde\beta)=Q(\tilde\beta)(3)
    其中等号成立的充要条件为X(β~β)2=0||X(\tilde\beta-\beta)||^2= 0,由β\beta的任意性可知β~\tilde\beta满足(2)式,因此β~\tilde\betaβ\beta的最小二乘估计。
  2. β^\hat\betaβ\beta的最小二乘估计,β~\tilde\beta是正规方程组的解,由(2)式:
    Q(β^)Q(β~)Q(\hat\beta)\le Q(\tilde\beta)由(3)式:
    Q(β~)Q(β^)Q(\tilde\beta)\le Q(\hat\beta)所以
    Q(β^)=Q(β~)Q(\hat\beta)=Q(\tilde\beta)这意味着
    (β~β)XX(β~β)=X(β~β)2=0(\tilde\beta-\beta)'X'X(\tilde\beta-\beta)=||X(\tilde\beta-\beta)||^2= 0因此Xβ~=Xβ^X\tilde\beta=X\hat\beta因为β~\tilde\beta满足正规方程,则
    XXβ^=XXβ~=XYX'X\hat\beta=X'X\tilde\beta=X'Yβ^\hat\beta是正规方程组的解。

最小二乘估计的性质

设Y,Z为随机向量,A,B为常数矩阵,有下面两个常用结论:

  1. E(AY)=AE(Y)E(AY)=AE(Y)
  2. cov(AY,BZ)=Acov(Y,Z)BTcov(AY,BZ)=Acov(Y,Z)B^T

性质一 β^\hat\betaβ\beta线性无偏估计
proof.proof.
由(1),β^\hat\beta是Y的线性函数,故为线性估计。
E(β^)=E((XX)1XY)=(XX)1XE(Y)=(XX)1XXβ=βE(\hat\beta)=E((X'X)^{-1}X'Y)=(X'X)^{-1}X'E(Y)=(X'X)^{-1}X'X\beta=\beta

性质二 cov(β^,β^)=σ2(XX)1=σ2N1cov(\hat\beta,\hat\beta)=\sigma^2(X'X)^{-1}=\sigma^2N^{-1}
proof.proof.
cov(β^,β^)=cov((XX)1XY,(XX)1XY)=(XX)1Xcov(Y,Y)X(XX)1=σ2(XX)1XX(XX)1=σ2(XX)1\begin{aligned} cov(\hat\beta,\hat\beta)&=cov((X'X)^{-1}X'Y,(X'X)^{-1}X'Y)\\ &=(X'X)^{-1}X'cov(Y,Y)X(X'X)^{-1}\\ &=\sigma^2(X'X)^{-1}X'X(X'X)^{-1}\\ &=\sigma^2(X'X)^{-1} \end{aligned}
可见β^\hat\beta的各分量在一般情况下并不独立

对任一k+1维向量C=(c0,c1,,ck)C=(c_0,c_1,\ldots,c_k)',若存在n维向量L,使E(LY)=CβE(L'Y)=C'\beta,则称CβC'\beta为可估函数,而可估函数CβC'\beta的最小方差线性无偏估计,称为它的最好线性无偏估计(Best Linear Unbiased Estimate,BLUE)

性质三 (Gauss-Markov定理):Cβ^C'\hat\betaCβC'\beta的最好线性无偏估计,其中β^\hat\betaβ\beta的最小二乘估计,
proof.proof.
由性质1易见Cβ^C'\hat\betaCβC'\beta的无偏估计,它显然是Y的线性函数,只需证明对CβC'\beta的任一线性无偏估计T=LYT=L'Y,有var(T)var(Cβ^)var(T)\ge var(C'\hat\beta).
T=LYT=L'YCβC'\beta的无偏估计,则对任意β\beta有:
E(T)=E(LY)=LE(Y)=LXβ=CβE(T)=E(L'Y)=L'E(Y)=L'X\beta=C'\beta
LX=C\therefore L'X=C'
注意到性质二,有
var(T)=var(LY)=Lcov(Y,Y)L=σ2LLvar(T)=var(L'Y)=L'cov(Y,Y)L=\sigma^2L'L
var(Cβ)=Ccov(β^,β^)C=σ2C(XX)1Cvar(C'\beta)=C'cov(\hat\beta,\hat\beta)C=\sigma^2C(X'X)^{-1}C

0LX(XX)1C2=LLC(XX)1Cvar(Cβ^)var(T)\begin{aligned} 0\le&||L-X(X'X)^{-1}C||^2\\ =&L'L-C'(X'X)^{-1}C\\ \therefore& var(C'\hat\beta)\le var(T) \end{aligned}
由于T是CβC'\beta的线性无偏估计,所以Cβ^C'\hat\betaCβC'\beta的BLUE。

Gauss-Markov定理指出,Cβ^C'\hat\betaCβC'\beta的一切线性无偏估计中是方差最小的,但在CβC'\beta的一切无偏估计中不一定方差最小。如果在正态性条件下,,Cβ^C'\hat\beta是UMVUE。

性质四 在正态条件下,Cβ^C'\hat\betaCβC'\beta的一致最小方差无偏估计(UMVUE)。(UMVUE定义参看数理统计基础一文)
该性质证明需要用到的结论较多,故不证。

σ2\sigma^2的估计

由最小二乘估计原理知道,在线性模型(Y,Xβ,σ2In)(Y,X\beta,\sigma^2I_n)中,β\beta用它的最小二乘估计β^\hat\beta代替时,Q达到最小,记
Y^=Xβ^\hat Y=X\hat\beta
表示n个试验点处Y的回归值,
ε^=YY^=YXβ^\hat\varepsilon=Y-\hat Y=Y-X\hat\beta
表示实际观测值Y与它的回归值之差,称为残差,关于残差有如下性质

性质五

  1. E(ε^)=0E(\hat\varepsilon)=\mathbf{0}
  2. cov(ε^,ε^)=σ2(InX(XX)1X)cov(\hat\varepsilon,\hat\varepsilon)=\sigma^2(I_n-X(X'X)^{-1}X')
  3. cov(β^,ε^)=0cov(\hat\beta,\hat\varepsilon)=0

1.式显然成立,下证另外两式。

ε^=YXβ^=YX(XX)1XY=(InX(XX)1X)Y=AY\begin{aligned} \hat\varepsilon=&Y-X\hat\beta=Y-X(X'X)^{-1}X'Y\\ =&(I_n-X(X'X)^{-1}X')Y\overset{\triangle}{=}AY \end{aligned}
这里A=InX(XX)1XA=I_n-X(X'X)^{-1}X'
不难验证A是对称幂等阵
A=A,A2=AA'=A,A^2=A
cov(ε^,ε^)=cov(AY,AY)=Acov(Y,Y)A=σ2AA=σ2A=σ2(InX(XX)1X)cov(β^,ε^)=cov((XX)1XY,AY)=(XX)1Xcov(Y,Y)A=σ2(XX)1XA=0\begin{aligned} cov(\hat\varepsilon,\hat\varepsilon)=&cov(AY,AY)=Acov(Y,Y)A'\\ =&\sigma^2AA'=\sigma^2A=\sigma^2(I_n-X(X'X)^{-1}X')\\ cov(\hat\beta,\hat\varepsilon)=&cov((X'X)^{-1}X'Y,AY)\\ =&(X'X)^{-1}X'cov(Y,Y)A'\\ =&\sigma^2(X'X)^{-1}X'A=0 \end{aligned}

几何意义

下面给出最小二乘估计几β^\hat\beta与残差ε^\hat\varepsilon的几何意义。如果把某个随机变量的n个观测值看成n维欧氏空间的一个向量,在此空间中,向量Y的长度定义为Y=YY||Y||=\sqrt{Y'Y},两个向量的距离定义为Y1Y2||Y_1-Y_2||

记设计矩阵的列向量为X0,X1,,XkX_0,X_1,\cdots,X_k,是n维欧氏空间中的k+1个向量,它们的线性组合构成了n维空间的一个线性子空间,记为L(X)\mathscr L(X),对任一向量β^,Xβ^L(X)\hat\beta,X\hat\beta\in\mathscr L(X),因此β的最小二乘估计β^\hat\beta就是在L(X)\mathscr L(X)中寻找一个向量Xβ^X\hat\beta使得相应的ε^\hat\varepsilon长度最短,这仅当Xβ^X\hat\betaYYL(X)\mathscr L(X)中的投影时才能达到,如图由Xβ^=X(XX)1XYX\hat\beta=X(X'X)^{-1}X'Y,可以称:
P=X(XX)1XP=X(X'X)^{-1}X'

为空间L(X)\mathscr L(X)上的投影阵(projective matrix),或帽子阵(hat matirx),容易看出,投影阵具有对称性,幂等性。(投影阵性质特殊,可以将(n+1)维欧氏空间的向量投影到n维超平面上,是向量投影运算的一般化,在其他领域如CV等应用广泛,请特别留意)
β^\hat\betaβ\beta的最小二乘估计时,ε^\hat\varepsilon表示Y到L(X)\mathscr L(X)的垂线,性质五的3式表示Y^\hat Yε^\hat\varepsilon相互垂直,

记残差向量的长度平方,即
Qe=ε^2=ε^ε^Q_e=||\hat\varepsilon||^2=\hat\varepsilon'\hat\varepsilon
为残差平方和
Qe=ε^ε^=(YXβ^)(YXβ^)=(AY)(AY)=YAY=YYYX(XX)1(XX)(XX)1XY=YYβ^XXβ^=YYY^Xβ^=YYY^Y^\begin{aligned} Q_e=&\hat\varepsilon'\hat\varepsilon=(Y-X\hat\beta)'(Y-X\hat\beta)=(AY)'(AY)\\ =&Y'AY=Y'Y-Y'X(X'X)^{-1}(X'X)(X'X)^{-1}X'Y\\ =&Y'Y-\hat\beta'X'X\hat\beta=Y'Y-\hat Y'X\hat\beta\\ =& Y'Y-\hat Y'\hat Y \end{aligned}
(4)
上式说明残差向量ε^\hat\varepsilon与估计量Y^\hat Y的长度平方和等于观测向量YY的长度平方,同时也给出了QeQ_e的不同表达式。

残差ε^\hat\varepsilon与随机误差σ2\sigma^2有关,所以用Qe=ε^2Q_e=||\hat\varepsilon||^2来估计σ2\sigma^2是合理的。

证明残差的平方和与σ2\sigma^2的无偏估计之间关系,要用到下面三个结论

  1. 设n维随机向量Y,有E(Y)=a,cov(Y,Y)=σ2In,AE(Y)=a,cov(Y,Y)=\sigma^2I_n,A为n阶对称常数阵,有
    E(YAY)=aAa+σ2tr(A)E(Y'AY)=a'Aa+\sigma^2tr(A)
  2. 设A,B是两个使乘积AB,BA都为方阵的矩阵,则
    tr(AB)=tr(BA)tr(AB)=tr(BA)
  3. tr(A+B)=tr(A)+tr(B)tr(A+B)=tr(A)+tr(B)

性质六
σ^2=Qen(k+1)=Qent\hat\sigma^2=\frac{Q_e}{n-(k+1)}=\frac{Q_e}{n-t}
称为残差方差(Residual Variance),有
E(σ^2)=σ2E(\hat\sigma^2)=\sigma^2
proof.proof.
由(4)式及上面三个结论,有:
E(Qe)=E(YAY)=βX[InX(XX)1X]Xβ+σ2tr[InX(XX)1X]=σ2tr(InX(XX)1X)=σ2[ntr[X(XX)1X]]=σ2[ntr(Ik+1)]=σ2(nk1)\begin{aligned} E(Q_e)=& E(Y'AY)\\ =& \beta'X'[I_n-X(X'X)^{-1}X']X\beta+\sigma^2tr[I_n-X'(X'X)^{-1}X']\\ =& \sigma^2tr(I_n-X(X'X)^{-1}X')\\ =& \sigma^2[n-tr[X(X'X)^{-1}X']]\\ =& \sigma^2[n-tr(I_{k+1})]\\ =& \sigma^2(n-k-1) \end{aligned}
由此即得
E(σ^2)=E(Qenk1)=σ2E(\hat\sigma^2)=E(\frac{Q_e}{n-k-1})=\sigma^2

性质七 在正态性条件下

  1. β^,ε^\hat\beta,\hat\varepsilon相互独立,且β^N(β,σ2N1),ε^N(0,σ2A)\hat\beta\sim N(\beta,\sigma^2N^{-1}),\hat\varepsilon\sim N(\mathbf 0,\sigma^2 A)
  2. β^,Qe\hat\beta,Q_e相互独立
  3. Qeσ2χ2(nk1)\frac{Q_e}{\sigma^2}\sim\chi^2(n-k-1)

由性质一,二,五知1.2.显然成立下证3.
proof.proof.
Qe=YAY=(Xβ+ε)A(Xβ+ε)=βXAXβ+βXAε+εAXβ+εAε\begin{aligned} Q_e=&Y'AY=(X\beta+\varepsilon)'A(X\beta+\varepsilon)\\ =& \beta'X'AX\beta+\beta'X'A\varepsilon+\varepsilon'AX\beta+\varepsilon'A\varepsilon \end{aligned}
注意到A=InX(XX)1XA=I_n-X(X'X)^{-1}X',容易验证
βXAXβ=0\beta'X'AX\beta=0 βXA=AXβ=0\beta'X'A=AX\beta=0
因此
Qe=εAεQ_e=\varepsilon'A\varepsilon
可见,残差平方和是随机误差ε\varepsilon的二次型。
因为矩阵A是对称幂等阵,因此一定存在一个n阶正交阵Γ\Gamma,使A=ΓTΛΓ,A=\Gamma^T\Lambda\Gamma,其中Λ=diag(λ1,,λn),λ1,,λn\Lambda=diag(\lambda_1,\cdots,\lambda_n),\lambda_1,\cdots,\lambda_n是A的特征根,且λi\lambda_i非0即1,非零个数为rank(A)=tr(A)=nk1,rank(A)=tr(A)=n-k-1,不妨设为λ1==λnk1=1\lambda_1=\cdots=\lambda_{n-k-1}=1,记e=Γε/σe=\Gamma\varepsilon/\sigma,由εN(0,σ2In)\varepsilon\sim N(0,\sigma^2I_n)Γ\Gamma的正交性可知eN(0,In)e\sim N(0,I_n),即e=(e1,,en)e=(e_1,\cdots,e_n)'中每个分量eie_i独立,都服从N(0,1)N\sim(0,1),故
Qeσ2=(εσ)ΓΓAΓΓ(εσ)=eΛe=i=1nk1ei2χ2(nk1)\frac{Q_e}{\sigma^2}=(\frac{\varepsilon}{\sigma})'\Gamma'\Gamma A\Gamma'\Gamma(\frac{\varepsilon}{\sigma})=e'\Lambda e=\sum_{i=1}^{n-k-1}e_i^2\sim\chi^2(n-k-1)

性质八εN(0,σ2In)\varepsilon\sim N(0,\sigma^2I_n),则β\beta的最小二乘估计β^\hat\beta也是β\beta的极大似然估计,σ2\sigma^2的极大似然估计为Qe/nQ_e/n
proof.proof.
εN(0,σ2In),YN(Xβ,σ2In)\varepsilon\sim N(0,\sigma^2I_n),Y\sim N(X\beta,\sigma^2I_n),有定义,Y有密度函数
f(Y;β,σ2)=(2πσ2)n/2exp[12σ2(YXβ)(YXβ)]f(Y;\beta,\sigma^2)=(2\pi\sigma^2)^{-n/2}\exp[-\frac{1}{2\sigma^2}(Y-X\beta)'(Y-X\beta)]似然函数
lnL(β,σ2)=n2ln2πn2lnσ212σ2(YXβ)(YXβ)\ln L(\beta,\sigma^2)=-\frac{n}{2}\ln 2\pi-\frac{n}{2}\ln\sigma^2-\frac{1}{2\sigma^2}(Y-X\beta)'(Y-X\beta)
因此
{lnLβ=12σ2(2XY+2XXβ)lnLσ2=n2σ2+12σ4(YXβ)(YXβ)\begin{cases} \frac{\partial \ln L}{\partial \beta}=-\frac{1}{2\sigma^2}(-2X'Y+2X'X\beta)\\ \frac{\partial \ln L}{\partial\sigma^2}=-\frac{n}{2\sigma^2}+\frac{1}{2\sigma^4}(Y-X\beta)'(Y-X\beta) \end{cases}
令它们为0,解得它们的极大似然估计为:
β^L=(XX)1XY\hat\beta_L=(X'X)^{-1}X'Y σ^L2=1n(YXβ^L)(YXβ^L)=Qen\hat\sigma^2_L=\frac{1}{n}(Y-X\hat\beta_L)'(Y-X\hat\beta_L)=\frac{Q_e}{n}

广义最小二乘估计

在线性回归模型(Y,Xβ,σ2In)(Y,X\beta,\sigma^2I_n)中,我们假定各次观测独立进行,即
cov(Y,Y)=σ2Incov(Y,Y)=\sigma^2I_n
我们考虑更一般的情况
cov(Y,Y)=σ2Qcov(Y,Y)=\sigma^2 Q
其中Q是已知的对称阵,且Q0|Q|\neq 0,将相应的线性回归模型记为(Y,Xβ,σ2Q)(Y,X\beta,\sigma^2Q)

定义:设A是一个n阶复矩阵,如果存在一个n阶复矩阵B使A=B2A=B^2,则称B是A的平方根矩阵,记为B=AB=\sqrt A

为了求未知参数β,σ2\beta,\sigma^2的最小二乘估计,作变换Z=Q1/2Y,U=Q1/2XZ=Q^{-1/2}Y,U=Q^{-1/2}X,那么
E(Z)=E(Q1/2Y)=Q1/2E(Y)=Q1/2Xβ=UβE(Z)=E(Q^{-1/2}Y)=Q^{-1/2}E(Y)=Q^{-1/2}X\beta=U\beta cov(Z,Z)=cov(Q1/2Y,Q1/2Y)=σ2Q1/2QQ1/2=σ2Incov(Z,Z)=cov(Q^{-1/2}Y,Q^{-1/2}Y)=\sigma^2Q^{-1/2}\cdot Q\cdot Q^{-1/2}=\sigma^2I_n
这样,线性回归模型(Y,Xβ,σ2Q)(Y,X\beta,\sigma^2 Q)便化为(Z,Uβ,σ2In)(Z,U\beta,\sigma^2I_n),这是前面已经讨论的情况,可得正规方程组
UUβ=UZU'U\beta=U'Z XQ1Xβ=XQ1YX'Q^{-1}X\beta=X'Q^{-1}Y
P=Q1P=Q^{-1}
那么β^=(XPX)1XPY\hat\beta=(X'PX)^{-1}X'PY
这就是β\beta的最小二乘估计,称为广义最小二乘估计。
那么
cov(β^,β^)=σ2(UU)1=σ2(XPX)1cov(\hat\beta,\hat\beta)=\sigma^2(U'U)^{-1}=\sigma^2(X'PX)^{-1}
残差平方和
Qe=ZUβ^2=Q1/2YQ1/2X(XPX)1XPY2=YPYYPX(XPX)1XPY=YPYYPXβ^\begin{aligned} Q_e=||Z-U\hat\beta||^2=&||Q^{-1/2}Y-Q^{-1/2}X(X'PX)^{-1}X'PY||^2\\ =& Y'PY-Y'PX(X'PX)^{-1}X'PY\\ =& Y'PY-Y'PX\hat\beta \end{aligned}

平方和分解公式

记:
Syy=(yiy¯)2=Y1y¯2S_{yy}=\sum(y_i-\bar y)^2=||Y-\mathbf 1\bar y||^2称为总变差平方和(总体平方和,Total Sum of Squares,TSS),
Qe=YY^2Q_e=||Y-\hat Y||^2称为残差平方和(Residual Sum of Squares,RSS),U=Y^1y¯2=(y^iy¯)2U=||\hat Y-\mathbf 1\bar y||^2=\sum (\hat y_i-\bar y)^2称为表示回归值y^i\hat y_i的波动,称为回归平方和(Explained Sum of Squares,ESS,Sum of Squares of Regression)

Syy=Qe+US_{yy}=Q_e+U
称为平方和分解公式
事实上:
Syy=Y1y¯2=(YY^)+(Y^1y¯)2=YY^2+Y^1y¯2+2(YY^)(Y^1y¯)\begin{aligned} S_{yy}=& ||Y-1\bar y||^2=||(Y-\hat Y)+(\hat Y-1\bar y)||^2\\ =&||Y-\hat Y||^2+||\hat Y-1\bar y||^2+2(Y-\hat Y)'(\hat Y-1\bar y) \end{aligned}
引用之前的A,P(投影阵)记号,显然AP=0AP=0,又(YY^)1=0(Y-\hat Y)'\mathbf{1}=0
(YY^)(Y^1y¯)=(YY^)Y^(YY^)1y¯={[InX(XX)1X]Y}[X(XX)1XY]=YAPY=0\begin{aligned} (Y-\hat Y)'(\hat Y-1\bar y)=&(Y-\hat Y)'\hat Y-(Y-\hat Y)'1\bar y\\ =& \{[I_n-X(X'X)^{-1}X']Y \}'[X(X'X)^{-1}X'Y]\\ =& Y'APY=0 \end{aligned}
得证。

参考资料

  1. 应用数理统计(第二版) 关静,张玉环,史道济 主编