最优估计准则与方法(5)加权最小二乘估计(WLS)_学习笔记
前言
最优估计理论中研究的最小二乘估计(LS)为线性最小二乘估计(LLS),包括古典最小二乘估计(CLS)[1]、加权最小二乘估计(WLS)和递推最小二乘估计(RLS)。本文将详细介绍加权最小二乘估计(WLS)。
线性参数估计问题描述
这里重复文章[1]的相关描述。设XXX为nnn维未知参数向量,ZZZ为kkk维观测向量,表示经过kkk组实验观测得到的观测值向量,其中元素ziz_{i}zi表示第i次观测实验得到的观测值,显然其是1维观测标量,VVV为kkk维观测噪声向量,其中元素viv_{i}vi表示第i次观测实验的观测噪声,显然其是1维噪声标量。一般情况下k>nk > nk>n且希望kkk比nnn大得多。单次观测值为多维的情况将在其他篇幅讨论。观测实验依据的自变量为θ\thetaθ,则将观测量ziz_{i}zi表示为关于θ\thetaθ的未知函数f(θ,X)f(\theta,X)f(θ,X):
zi=f(θ,X)=∑j=1n[xjhi,j(θ)]+vi=x1hi,1(θ)+x2hi,2(θ)+⋯+xnhi,n(θ)+vi\begin{align*} z_{i} = f(\theta,X) = \sum_{j=1}^{n} \left [ x_{j}h_{i,j}(\theta) \right ]+ v_{i} = x_{1}h_{i,1}(\theta)+ x_{2}h_{i,2}(\theta) + \cdots + x_{n}h_{i,n}(\theta) + v_{i} \tag{1} \\ \end{align*} zi=f(θ,X)=j=1∑n[xjhi,j(θ)]+vi=x1hi,1(θ)+x2hi,2(θ)+⋯+xnhi,n(θ)+vi(1)
其中
X=[x1x2⋮xn]Z=[z1z2⋮zk]V=[v1v2⋮vk]\begin{align*} X = \begin{bmatrix} x_{1} \\ x_{2} \\ \vdots \\ x_{n} \end{bmatrix} Z = \begin{bmatrix} z_{1} \\ z_{2} \\ \vdots \\ z_{k} \end{bmatrix} V = \begin{bmatrix} v_{1} \\ v_{2} \\ \vdots \\ v_{k} \end{bmatrix} \end{align*} X=x1x2⋮xnZ=z1z2⋮zkV=v1v2⋮vk
式(1)中hi,j(θ)h_{i,j}(\theta)hi,j(θ)表示第iii次观测第jjj个基函数,常用为多项式、三角函数或自然指数函数形式:
hi,j(θ)=θj−1hi,j(θ)=sin(jθ)hi,j(θ)=exp(λjθ)\begin{align*} h_{i,j}(\theta) &= \theta ^{j-1} \\ h_{i,j}(\theta) &= sin(j\theta) \\ h_{i,j}(\theta) &= exp(\lambda_{j} \theta) \\ \end{align*} hi,j(θ)hi,j(θ)hi,j(θ)=θj−1=sin(jθ)=exp(λjθ)
其中,λj\lambda_{j}λj为自然数指数参数。
当观测实验进行,上述基函数均可根据θ\thetaθ求得。令hi=[hi,1(θ)hi,2(θ)⋯hi,n(θ)]h_{i} = \begin{bmatrix} h_{i,1}(\theta) & h_{i,2}(\theta) & \cdots & h_{i,n}(\theta) \\ \end{bmatrix}hi=[hi,1(θ)hi,2(θ)⋯hi,n(θ)]且为已知,其为nnn维常向量,将式(1)改写为:
Z=HX+V\begin{align*} Z= HX+ V \tag{2} \\ \end{align*} Z=HX+V(2)
其中,HHH为参数向量XXX到观测向量ZZZ的k×nk \times nk×n维转移矩阵:
H=[h1h2⋮hk]=[h1,1(θ)h1,2(θ)⋯h1,n(θ)h2,1(θ)h2,2(θ)⋯h2,n(θ)⋮⋮⋱⋮hk,1(θ)hk,2(θ)⋯hk,n(θ)]\begin{align*} H = \begin{bmatrix} h_{1} \\ h_{2} \\ \vdots \\ h_{k} \end{bmatrix} = \begin{bmatrix} h_{1,1}(\theta) & h_{1,2}(\theta) & \cdots & h_{1,n}(\theta) \\ h_{2,1}(\theta) & h_{2,2}(\theta) & \cdots & h_{2,n}(\theta) \\ \vdots & \vdots & \ddots & \vdots\\ h_{k,1}(\theta) & h_{k,2}(\theta) & \cdots & h_{k,n}(\theta) \end{bmatrix} \\ \end{align*} H=h1h2⋮hk=h1,1(θ)h2,1(θ)⋮hk,1(θ)h1,2(θ)h2,2(θ)⋮hk,2(θ)⋯⋯⋱⋯h1,n(θ)h2,n(θ)⋮hk,n(θ)
显然,观测向量ZZZ与被估参数向量XXX存在线性关系,依据最优准则求对XXX的估计值X^\hat{X}X^是一个线性参数估计问题,自然对应线性最小二乘估计(LLS)。
这里讨论下超定方程组的矛盾:当k=nk = nk=n时,线性方程组有唯一精确解,但当k>nk > nk>n,线性方程数大于未知被估参数向量的维度,线性方程组变成线性超定方程组,其解不唯一。最小二乘法的思想是需求统计意义上的近似解,使线性超定方程组中各方程能得到近似相等。
加权最小二乘估计(Weighted Least Squares Estimation, WLSE)
最小二乘估计(LS) 假设每次观测量对于估计结果的影响程度相同,但实际上观测数据的权重与该次观测的残差平方呈反比更为合理,因此引出加权最小二乘估计(WLS)。
加权最小二乘估计(WLS) 估计准则为:加权残差平方和最小。
根据式(3)代价函数改写如下:
J=E^TWE^=(Z−HX^)TW(Z−HX^)=∑i=1kwie^i2=∑i=1kwi(zi−hiX^)2=min\begin{align*} J = \hat{E}^{T}W\hat{E} &= (Z-H\hat{X})^{T}W(Z-H\hat{X}) = \sum_{i=1}^{k} w_{i}\hat{e}_{i}^{2} = \sum_{i=1}^{k}w_{i}(z_{i}-h_{i}\hat{X})^{2}=min \tag{3} \\ \end{align*} J=E^TWE^=(Z−HX^)TW(Z−HX^)=i=1∑kwie^i2=i=1∑kwi(zi−hiX^)2=min(3)
其中,e^i\hat{e}_{i}e^i为第iii次观测的残差(Residual Error),E^\hat{E}E^为kkk维残差向量有:
e^i=zi−hiX^E^=Z−HX^=[e^1e^2⋮e^k]\begin{align*} \hat{e}_{i} &= z_{i}-h_{i}\hat{X} \\ \hat{E} &= Z-H\hat{X} = \begin{bmatrix} \hat{e}_{1} \\ \hat{e}_{2} \\ \vdots \\ \hat{e}_{k} \end{bmatrix} \\ \end{align*} e^iE^=zi−hiX^=Z−HX^=e^1e^2⋮e^k
WWW为可根据实际情况适当选取的k×kk\times kk×k阶对称正定加权矩阵,但当W=IW=IW=I时,加权最小二乘估计退化为最小二乘估计。
W=[w10⋯00w2⋯0⋮⋮⋱⋮00⋯wk]\begin{align*} W &= \begin{bmatrix} w_{1} & 0& \cdots& 0\\ 0& w_{2} & \cdots& 0 \\ \vdots& \vdots& \ddots & \vdots\\ 0& 0& \cdots& w_{k} \end{bmatrix} \\ \end{align*} W=w10⋮00w2⋮0⋯⋯⋱⋯00⋮wk
加权最小二乘估计(WLS)方法
根据式(3)进行对如下代价函数进行最小化:
J=(Z−HX^)TW(Z−HX^)\begin{align*} J &= (Z-H\hat{X})^{T}W(Z-H\hat{X}) \tag{4} \\ \end{align*} J=(Z−HX^)TW(Z−HX^)(4)
令JJJ对X^\hat{X}X^求偏导,并令其为0,有:
∂J∂X^=0∂J∂(Z−HX^)∂(Z−HX^)∂X^=0−2HTW(Z−HX^)=0X^=(HTWH)−1HTWZ\begin{align*} \frac{\partial J}{\partial \hat{X}} &= 0 \\ \frac{\partial J}{\partial (Z-H\hat{X})}\frac{\partial (Z-H\hat{X})}{\partial \hat{X}} &=0 \\ -2H^{T}W(Z-H\hat{X}) &= 0 \\ \hat{X} &= (H^{T}WH)^{-1}H^{T}WZ \tag{5} \end{align*} ∂X^∂J∂(Z−HX^)∂J∂X^∂(Z−HX^)−2HTW(Z−HX^)X^=0=0=0=(HTWH)−1HTWZ(5)
再由∂2J∂X^2=2HTWH>0\frac{\partial^{2} J}{\partial \hat{X}^{2}}=2H^{T}WH > 0∂X^2∂2J=2HTWH>0,为X^\hat{X}X^为被估参数向量XXX的加权最小二乘估计,显然其是观测向量ZZZ的线性估计。
Jmin=(Z−HX^)TW(Z−HX^)=ZT(I−H(HTWH)−1HTW)T(I−H(HTWH)−1HTW)Z\begin{align*} J_{min} &= (Z-H\hat{X})^{T}W(Z-H\hat{X}) \\ &= Z^{T}(I-H(H^{T}WH)^{-1}H^{T}W)^{T}(I-H(H^{T}WH)^{-1}H^{T}W)Z \tag{6} \\ \end{align*} Jmin=(Z−HX^)TW(Z−HX^)=ZT(I−H(HTWH)−1HTW)T(I−H(HTWH)−1HTW)Z(6)
加权最小二乘估计(WLS)无偏性
令估计误差为X~\tilde{X}X~,定义被估参数向量XXX与估计值向量X^\hat{X}X^的偏差,有:
X~=X−X^=(HTWH)−1HTWHX−(HTWH)−1HTWZ=(HTWH)−1HTW(HX−Z)=−(HTWH)−1HTWV\begin{align*} \tilde{X} &= X - \hat{X} \tag{7} \\ &= (H^{T}WH)^{-1}H^{T}WHX - (H^{T}WH)^{-1}H^{T}WZ \\ &= (H^{T}WH)^{-1}H^{T}W(HX - Z) \\ &= -(H^{T}WH)^{-1}H^{T}WV \tag{8} \\ \end{align*} X~=X−X^=(HTWH)−1HTWHX−(HTWH)−1HTWZ=(HTWH)−1HTW(HX−Z)=−(HTWH)−1HTWV(7)(8)
估计误差X~\tilde{X}X~的数学期望为:
E[X~]=E[X−X^]=E[−(HTWH)−1HTWV]=−(HTWH)−1HTWE[V]\begin{align*} E[\tilde{X}] &= E[X - \hat{X}] \tag{9} \\ &= E[-(H^{T}WH)^{-1}H^{T}WV] \\ &= -(H^{T}WH)^{-1}H^{T}WE[V] \tag{10} \\ \end{align*} E[X~]=E[X−X^]=E[−(HTWH)−1HTWV]=−(HTWH)−1HTWE[V](9)(10)
由式(10)可知,如果观测噪声VVV为白噪声,即E[V]=0E[V]=0E[V]=0,则加权最小二乘估计X^\hat{X}X^为无偏线性估计。在该无偏估计情况下,估计误差X^\hat{X}X^的方差矩阵与估计量X^\hat{X}X^的均方误差矩阵相等,推导见[1],即:
Var(X~)=MSE(X^)=E[X~X~T]=E[(−(HTWH)−1HTWV)(−(HTWH)−1HTWV)T]=(HTWH)−1HTWE[VVT]WH(HTWH)−1=(HTWH)−1HTWRWH(HTWH)−1\begin{align*} Var(\tilde{X}) &= MSE(\hat{X}) \tag{11} \\ &= E[\tilde{X}\tilde{X}^{T}] \\ &= E[(-(H^{T}WH)^{-1}H^{T}WV)(-(H^{T}WH)^{-1}H^{T}WV)^{T}] \\ &= (H^{T}WH)^{-1}H^{T}WE[VV^{T}]WH(H^{T}WH)^{-1} \\ &= (H^{T}WH)^{-1}H^{T}WRWH(H^{T}WH)^{-1} \tag{12} \\ \end{align*} Var(X~)=MSE(X^)=E[X~X~T]=E[(−(HTWH)−1HTWV)(−(HTWH)−1HTWV)T]=(HTWH)−1HTWE[VVT]WH(HTWH)−1=(HTWH)−1HTWRWH(HTWH)−1(11)(12)
其中,RRR为观测噪声向量VVV的方差矩阵:
R=[σ120⋯00σ22⋯0⋮⋮⋱⋮00⋯σk2]\begin{align*} R &= \begin{bmatrix} \sigma_{1}^{2} & 0& \cdots& 0\\ 0& \sigma_{2}^{2} & \cdots& 0 \\ \vdots& \vdots& \ddots & \vdots\\ 0& 0& \cdots& \sigma_{k}^{2} \end{bmatrix} \\ \end{align*} R=σ120⋮00σ22⋮0⋯⋯⋱⋯00⋮σk2
由式(8)和(14)可知,即使在无偏估计前提下,二者并不一定相等。因此,加权最小二乘无偏估计只能保证加权残差平方和最小,但不能保证估计误差方差最小。
最优加权最小二乘估计
由于WWW为可设定的对称正定加权矩阵,在无偏估计前提下,WWW取某个值可使估计误差方差矩阵式(12)最小,令R=CTCR=C^{T}CR=CTC,则:
Var(X~)=(HTWH)−1HTWRWH(HTWH)−1=(CWH(HTWH)−1)TCWH(HTWH)−1\begin{align*} Var(\tilde{X}) &= (H^{T}WH)^{-1}H^{T}WRWH(H^{T}WH)^{-1} \\ &= (CWH(H^{T}WH)^{-1})^{T} CWH(H^{T}WH)^{-1} \tag{13} \\ \end{align*} Var(X~)=(HTWH)−1HTWRWH(HTWH)−1=(CWH(HTWH)−1)TCWH(HTWH)−1(13)
令A=CWH(HTWH)−1A=CWH(H^{T}WH)^{-1}A=CWH(HTWH)−1,B=C−1HB=C^{-1}HB=C−1H,根据施瓦次(Schwarz)不等式:
Var(X~)=ATA≥(ATB)T(BTB)−1(BTA)=(HTR−1H)−1\begin{align*} Var(\tilde{X}) &= A^{T} A \geq (A^{T}B)^{T}(B^{T}B)^{-1} (B^{T}A) = (H^{T}R^{-1}H)^{-1} \tag{14} \\ \end{align*} Var(X~)=ATA≥(ATB)T(BTB)−1(BTA)=(HTR−1H)−1(14)
若式(14)取最小值,W=R−1W=R^{-1}W=R−1,此时有
X^=(HTR−1H)−1HTR−1ZVar(X~)=(HTR−1H)−1\begin{align*} \hat{X} &= (H^{T}R^{-1}H)^{-1}H^{T}R^{-1}Z \tag{15} \\ Var(\tilde{X}) &= (H^{T}R^{-1}H)^{-1} \tag{16} \\ \end{align*} X^Var(X~)=(HTR−1H)−1HTR−1Z=(HTR−1H)−1(15)(16)
当噪声向量VVV的统计均值为E[V]=0E[V]=0E[V]=0,且加权残差平方和中的最优加权矩阵W=R−1W=R^{-1}W=R−1时,最优加权最小二乘估计是缺少初值条件下的线性无偏最小方差估计,又称为马尔可夫(Markov)估计。
综上,根据加权最小二乘估计原理,做如下总结:
- 求加权最小二乘估计量X^\hat{X}X^不需要任何观测噪声向量VVV的任何统计信息;
- 加权最小二乘估计的无偏性取决于噪声向量VVV的数学期望,如VVV为白噪声,即为无偏估计;
- 无论是否具备无偏性,最小二乘估计只能保证加权残差平方和最小而不是估计误差方差最小;
- 当噪声向量VVV的均值为0,且已知其方差矩阵RRR,最优加权矩阵W=R−1W=R^{-1}W=R−1,此时为最优加权最小二乘估计,即马尔可夫估计。
参考文献
[1] 最优估计准则与方法(4)最小二乘估计(LS)_学习笔记
https://blog.csdn.net/jimmychao1982/article/details/149656745
[2] 《最优估计理论》,周凤歧,2009,高等教育出版社。
[3] 《最优估计理论》,刘胜,张红梅著,2011,科学出版社。