第5章〓最小方差无偏估计


对于未知常数的估计，第4章介绍的CRLB的计算，有可能得到有效估计量，同时也是方差最小的无偏估计量，但有效估计量的获得是有条件的，若有效估计量不存在，则还可以求得最小方差无偏(Minimum Variance Unbiased，MVU)估计。要获得好的估计，对信息的利用一定是充分的，即一般的MVU估计一定是根据充分统计量来求得的。本章讨论如何根据充分统计量求取最小方差无偏估计。






5.1最小方差无偏估计的定义

估计的均方误差可以表示为


Mse(θ^)=E{［θ^-θ］2}
=E［(θ^-E(θ^))2］+［E(θ^)-θ］2
=Var(θ^)+［E(θ^)-θ］2(5.1.1)


其中，Var(θ^)=E［(θ^-E(θ^))2］表示估计量的方差，b=E(θ^)-θ表示估计的偏差项，即估计的均方误差等于估计的方差加上偏差项的平方。很显然，对于无偏估计而言，偏差项b=0，估计的均方误差就等于估计的方差。但是，对于有偏估计，估计的均方误差和估计的方差并不相等。对于未知常数的估计，一般不宜直接采用最小均方估计。例如，对于例4.2所述的高斯白噪声中未知常数的估计问题，若假定估计为θ^=αz-，其中z-=1N∑N-1i=0zi为样本均值，由于E(θ^)=αθ，Var(θ^)=α2Var(z-)=α2σ2/N，b=E(θ^)-θ=(α-1)θ，所以，


Mse(θ^)=α2σ2N+(α-1)2θ2
Mse(θ^)α=2ασ2N+2(α-1)θ2


令Mse(θ^)α=0，可解得最佳系数为


αopt=θ2θ2+σ2/N


最佳系数αopt与被估计量θ有关，而θ是未知的，所以，得到的最佳估计是不可用的。

要得到可用的最佳的估计，可约束偏差项b=0，即约束估计为无偏估计的条件下，使均方误差


Mse(θ^)=E［(θ^-θ)2］=∫+∞-∞(θ^-θ)2p(z)dz(5.1.2)


达到最小，这样的估计称为最小方差无偏(MVU)估计。

5.2RBLS定理

第4章介绍的有效估计量，它是一个无偏估计量，而且它的方差达到CRLB，它的方差是最小的，所以，有效估计量就是一个MVU估计。若有效估计量不存在，MVU估计需要寻找另外的方法求解，RBLS（RaoBlackwellLehmannScheffe）定理就是根据充分统计量求解MVU估计的一种方法。

定理： 若θ是θ的一个无偏估计量，T(z)是θ的充分统计量，则θ^=E(θ|T(z))是： 

(1) θ的一个可用的估计； 

(2) 无偏估计； 

(3) 对所有的θ，方差小于或等于θ的方差； 

若充分统计量T(z)是完备的，则θ^=E(θ|T(z))是MVU估计量。

定理的证明从略。定理中完备的含义是指只存在唯一的T(z)的函数，使θ^=E(θ|T(z))是无偏的。

【例5.1】高斯白噪声中未知常数的估计，假定观测为


zi=θ+wi,i=0,1,…,N-1


其中，{wi}是零均值高斯白噪声序列，方差为σ2，θ是未知常数，求θ的MVU估计。

解：  首先找一个无偏估计量，很显然，θ=z0是无偏的。其次，求θ的充分统计量，由例4.12可知，T(z)=∑N-1i=0zi是θ的充分统计量。接着求条件均值θ^=E(θ|T(z))。由高斯随机变量的理论，假定X,Y是联合高斯随机变量，有


E(X|Y)=E(X)+Cov(X,Y)［Var(Y)］-1［Y-E(Y)］(5.1.3)



很显然，θ和T(z)是联合高斯随机变量，由于T(z)~N(Nθ,Nσ2)，


Cov［θ,T(z)］=E(z0-θ)∑N-1i=0zi-Nθ=Ew0∑N-1i=0wi=σ2
θ^=E［θ|T(z)］=θ+σ2(Nσ2)-1∑N-1i=0zi-Nθ=1N∑N-1i=0zi


很容易验证，T(z)=∑N-1i=0zi是完备的，所以，θ^=1N∑N-1i=0zi是MVU估计。

由于完备的充分统计量只存在一个唯一的函数使其无偏，所以最小方差无估计量也可以通过下面的方法求解。

假定T(z)是完备的充分统计量，则θ^=g［T(z)］，只要找到一个函数g(·)，使之满足


E{g［T(z)］}=θ(5.1.4)


在例5.1中，T(z)=∑N-1i=0zi，E［T(z)］=E∑N-1i=0zi=Nθ，所以，很容易找到这个函数为g(x)=xN，θ^=g［T(z)］=1N∑N-1i=0zi。

【例5.2】假定观测为


zi=wi，i=0,1,…,N-1


其中{wi}是独立同分布噪声，且wi在(0,2θ)区间服从均匀分布，求θ的MVU估计。

解：    首先要找一个无偏估计量，很显然，θ=1N∑N-1i=0zi是无偏的，其次，要求一个充分统计量，由习题4.13可知，充分统计量为T(z)=max{z0,z1，…,zN-1}，接下来求一个充分统计量的函数，使其无偏。

T(z)的分布函数为


FT(t)=P{T(z)≤t}=P{z0≤t,…,zN-1≤t}
=∏N-1i=0P{zi≤t}=(P{zi≤t})N


而


P{zi≤t}=
0，t≤0
t2θ，0<t≤2θ
1，t>2θ


T(z)的概率密度为


pT(t)=FT(t)t=N(P{zi≤t})N-1P{zi≤t}t=
0，t≤0
Nt2θN-112θ，0<t≤2θ
0，t>2θ
E(T)=∫2θ0tpT(t)dt=∫2θ0tNt2θN-112θdt=2NN+1θ


因此，


θ^=N+12Nmax{z0,z1,…,zN-1}


统计量T(z)的完备性参见习题5.3。






5.3线性最小方差无偏估计

RBLS定理给出了求解MVU估计的一种方法，但求解的过程是比较复杂的，并且不能保证每次都能得到MVU估计，这时可以考虑采用线性最小方差无偏估计(Best Linear Unbiased Estimator，BLUE)。BLUE是一种线性估计，在约束估计为无偏的条件下，使估计的方差最小。

假定估计为

θ^=∑N-1i=0aizi=aTz(5.3.1)


其中a=［a0a1…aN-1］T，z=［z0z1…zN-1］T，估计的方差为


Var(θ^)=E{［aTz-aTE(z)］2}=E{［aT(z-E(z))］2}

=E［aT(z-E(z))(z-E(z))Ta］=aTCza(5.3.2)


其中Cz为观测矢量z的方差阵。

注意，根据无偏的约束条件，E(θ^)=∑N-1i=0aiE(zi)=θ，要求E(zi)是θ的线性函数，即


E(zi)=siθ，i=0,1,…,N-1(5.3.3)


其中，si(i=0,1,…,N-1)是已知量，否则有可能不满足约束条件。例如，如果E(zi)=cosθ，显然不存在满足约束条件∑N-1i=0aicosθ=θ的系数。

由式(5.3.3)，无偏的约束条件可表示为


E(θ^)=∑N-1i=0aiE(zi)=∑N-1i=0aisiθ=θ


即


∑N-1i=0aisi=1(5.3.4)


用矢量表示为


aTs=1(5.3.5)


其中s=［s0s1…sN-1］T。

 构造目标函数，


J=aTCza+λ(aTs-1)(5.3.6)


对系数a求导，并令导数等于零，即


Ja=2Cza+λs=0(5.3.7)


由此可得最佳系数为


a=-λ2C-1zs(5.3.8)


将式(5.3.8)代入式(5.3.5)，可得-12λsTC-1zs=1，或者表示为-12λ=1sTC-1zs，于是，最佳系数可表示为

a=C-1zssTC-1zs(5.3.9)


线性最小方差无偏估计为


θ^=sTC-1zzsTC-1zs(5.3.10)


将式(5.3.9)代入式(5.3.2)，可得估计的方差为


Var(θ^)=sTC-1zCzC-1zs(sTC-1zs)2=sTC-1zs(sTC-1zs)2=1sTC-1zs(5.3.11)


【例5.3】高斯白噪声中未知常数的估计，假定观测为


zi=θ+wi,i=0,1,…,N-1


其中，θ为未知常数，{wi}为非平稳的零均值高斯噪声，且协方差为E(wiwj)=σ2iδij，δij=1，i=j
0，i≠j，求θ的线性最小方差无偏估计。

解：   根据无偏约束要求，


E(zi)=siθ，i=0,1,…,N-1


所以，si=1，i=0,1,…,N-1，s=［11…1］T=1N×1，测量噪声的方差阵为Cz=σ2diag{σ20,σ21,…,σ2N-1}，由式(5.3.10)得


θ^=sTC-1zzsTC-1zs=1TC-1zz1TC-1z1=∑N-1i=0ziσ2i∑N-1i=01σ2i


由式(5.3.11)得


Var(θ^)=1sTC-1zs=11TC-1z1=1∑N-1i=01σ2i







5.4信号处理实例——系统辨识

传递函数是描述惯性系统的重要工具，利用传递函数可进行系统预报、动态特性分析、分类等用途。通常，由系统的输入和输出来确定系统传递函数的过程也称为系统辨识。系统辨识在控制、通信、地质勘探、目标识别等领域有重要应用。系统辨识的典型方法包括阶跃响应法、脉冲响应法、相关分析法和最小二乘法等。本信号处理实例主要采用节拍延迟线模型阐述系统辨识的基本方法，尤其是线性最小方差无偏估计的应用。

节拍延迟线模型本质上可以看作一个FIR滤波器。输入信号un经过多步的延迟以后与抽头系数(加权系数)对应相乘并相加，即可得到系统的输出，如图5.1(a)所示。从系统模型的角度看，图5.1(a)可以简洁地表达为图5.1(b)，其中H(z)=∑p-1k=0hkz-k为系统函数，wn表示系统观测噪声。系统辨识就是根据输入序列un和输出序列zn估计出节拍延迟线模型的权系数hk。



图5.1系统辨识的基本模型




为简单起见，可将观测模型表达如下： 


zn=∑p-1k=0hkun-k+wn，n=0,1,…,N-1(5.4.1)


令z=［z0z1…zN-1］T，θ=［h0h1…hp-1］T，w=［w0w1…wN-1］T，H=u00…0
u1u0…0

uN-1uN-2…uN-p,于是式(5.4.1)可写为


z=Hθ+w(5.4.2)


其中，E(w)=0,Cov(w,w)=C。根据5.3节可知参数θ的线性最小方差无偏估计及其方差阵分别为


θ^=(HTC-1H)-1HTC-1z(5.4.3)
Cθ^=(HTC-1H)-1(5.4.4)


特别地，当C=σ2I时，θ^=(HTH)-1HTz，Cθ^=σ2(HTH)-1。可见，参数估计精度与矩阵H或输入信号{un}有关，提高参数估计精度可设计选取不同的输入信号。因此，在回归分析中矩阵H也称为设计矩阵。需要说明的是，如果输入信号给定了，最小方差无偏估计及其性能也就确定了。本实例讨论在输入信号可变的情况下如何降低估计量的方差，即形成最佳的设计矩阵。由于估计量θ^第i个分量θ^i的方差为


Var(θ^i)=［Cθ^］ii=σ2［(HTH)-1］ii(5.4.5)


直接寻求估计量的最小方差有一定困难。事实上，HTH是实对称矩阵，故存在正交阵V=［v0,…,vp-1］T使得HTH=VΛVT，其中Λ=diag(λ0,λ1,…,λp-1)。于是［Cθ^］ii=σ2［(HTH)-1］ii=σ2vTiΛ-1vi。最小化估计方差可通过约束vTivi=1的条件下，求出使σ2vTiΛ-1vi最小的vi。

利用拉格朗日乘因子法求解该优化问题。构造目标函数


J(vi,γ)=σ2vTiΛ-1vi+γ(1-vTivi)(5.4.6)


令Jvi=0,Jγ=0，可得(σ2Λ-1-γI)vi=0，vTivi=1，i=0,1,…,p-1。结合Λ的对角特性可知，γi=σ2λ-1i，vi=ei=［0…010…0］T。显然，HTH=Λ，minH［Cθ^］ii=σ2λ-1i。

进一步，由于［H］ij=ui-j,i≥j
0,i<j，则［HTH］ij=∑N-1n=min(i,j)un-iun-j≈∑N-1-|i-j|n=0unun+|i-j|。
令r(uu)k=1N∑N-1-kn=0unun+k表示序列{un}的自相关函数，则


HTH=Nr(uu)0r(uu)1…r(uu)p-1
r(uu)1r(uu)0…r(uu)p-2
……
r(uu)p-1r(uu)p-2…r(uu)0(5.4.7)


由式(5.4.5)、式(5.4.6)可知，最小化估计量的方差Var(θ^i)=［Cθ^］ii，就是要使得矩阵HTH对角化，即r(uu)k=r(uu)0δk，其中δk表示离散δ函数。通常，可选择输入信号un为伪随机白噪声序列。此时，HTH=Nr(uu)0I，θ^=(HTH)-1HTz,即节拍延迟线模型的抽头系数估计量可写为


h^i=1Nr(uu)0∑N-1n=iun-izn=1N∑N-1-in=0unzn+iruu0=r(uz)ir(uu)0(5.4.8)


其中，r(uz)i=1N∑N-1-in=0unzn+i为互相关函数，r(uu)0=σ2为噪声方差。可见，系统参数的最小方差无偏估计量由输入输出序列互相关和输入序列自相关共同决定。

事实上，式(5.4.8)也可从频谱的角度理解。在无噪声情况下，当N足够大时，由式(5.4.1)可知

Puz(f)=H(f)Puu(f)(5.4.9)


其中Puz(f)=∑+∞k=-∞r-(uz)ke-j2πfk，r-(uz)k=E(unzn+k)，H(f)=∑P-1k=0hke-j2πfk，Puu(f)=∑+∞k=-∞r(uz)ke-j2πfk。可见，式(5.4.8)的结果与式(5.4.9)也是相符的。二者的差异在于一个使用期望，一个使用样本自相关函数。总的来说，从时域和频域看二者是等效的。要想获得更好的参数估计，需要对输入的信号有一定的要求，即频谱要宽且没有零点。因此，伪随机白噪声序列是一个比较好的选择，其功率谱是常数。

下面通过仿真实验进行验证。观测模型如式(5.4.1)所示，模型参数设置如下： h0=1.0，h1=-0.2，p=2； 观测样本数N=100,300； 信噪比SNR在-10~20dB均匀变化，步长为3dB。通过蒙特卡洛仿真后，统计参数估计的均方根误差和CRLB进行比较，仿真结果如图5.2所示。




图5.2模型参数的估计性能曲线



由图5.2可以看出： ①参数估计的克拉美罗下限随着样本数的增加而降低，也随着信噪比的增大而减低； ②随着样本数的增加，系数估计误差根方差更加接近克拉美罗界。


习题

5.1什么是完备的充分统计量？如何利用完备的充分统计量估计未知参数？由充分统计量是否能够获得MVUE？

5.2有效估计量、充分统计量、MVUE之间是否有联系？可举例验证你的观点。

5.3假定观测为


zi=wi，i=0,1,…,N-1


其中{wi}是独立同分布噪声，且wi在(0,2θ)区间服从均匀分布。试证明T(z)=max{z0,z1，…,zN-1}为参数θ的充分完备统计量。

5.4设zn，n=0,1,…,N-1表示N次独立贝努利试验(投掷硬币)的观测，且P(zn=1)=θ，P(zn=0)=1-θ。试利用RBLS定理求θ的MVU估计量(RBLS定理对概率分布函数也成立，可以假定充分统计量是完备的)。

5.5考虑观测zn=Asn+B+wn，n=0,1,…,N-1，其中sn为已知信号，{wn}为零均值噪声序列，协方差阵为C。试求待估参数的BLUE及方差。特别地，若B已知，情况如何？

5.6观测模型为

zn=Acos2πf0n+wn，n=0,1,…,N-1

其中wn是零均值高斯白噪声，方差为σ2，f0、σ2均已知。

(1) 求参数A的线性最小方差无偏估计量； 

(2) 参数A的最小方差无偏估计量是否存在？若存在，试求之。