最小二乘法的无偏估计

最小二乘法的无偏估计

目前我已知的最小二乘估计法有两种,第一种是基本的最小二乘法,对于白噪声,这类方法可以得到无偏估计。但对于有色噪声,这类方法只能得到有偏估计。为了解决这个问题,就导致了第二类最小二乘法的产生。这类改进算法可以在有色噪声下也能得到无偏估计。

  1. 噪声视为白噪声的最小二乘法
    • 一般最小二乘法
    • 加权最小二乘法
    • 递推最小二乘法(RLS)
    • 渐消记忆RLS法
  2. 噪声视为有色噪声的最小二乘法
    • 广义最小二乘法(GLS)
    • 增广最小二乘法(RELS)
    • 多级最小二乘法(MSLS)

这篇主要介绍改进后的最小二乘法。 在进入正题之前,先简要阐述白噪声和有色噪声的区别:

  • 白噪声:不同时刻的噪声是不相关的,自相关函数为脉冲函数;
  • 有色噪声:不同时刻的噪声之间存在相关性。在工程实践中,往往是这类噪声。

1.广义最小二乘法(GLS)

1.1.系统模型

系统模型如下所示:
A ( q − 1 ) y ( k ) = B ( q − 1 ) u ( k ) + ε ( k ) , ε ( k ) = ξ ( k ) C ( q − 1 ) (1) A(q^{-1})y(k)=B(q^{-1})u(k)+\varepsilon(k), \varepsilon(k)=\frac{\xi(k)}{C(q^{-1})} \tag 1 A(q1)y(k)=B(q1)u(k)+ε(k),ε(k)=C(q1)ξ(k)(1)
上式中, ξ ( k ) \xi(k) ξ(k)为白噪声, ε ( k ) \varepsilon(k) ε(k)为有色噪声序列。三个多项式形式如下所示:
A ( q − 1 ) = 1 + a 1 q − 1 + . . . + a n a q − n a B ( q − 1 ) = b 0 + b 1 q − 1 + . . . + b n b q − n b C ( q − 1 ) = 1 + c 1 q − 1 + . . . + c n c q − n c (2) \begin{aligned} A(q^{-1})=1+a_1q^{-1}+...+a_{n_a}q^{-{n_a}} \tag 2\\ B(q^{-1})=b_0+b_1q^{-1}+...+b_{n_b}q^{-{n_b}} \\ C(q^{-1})=1+c_1q^{-1}+...+c_{n_c}q^{-{n_c}} \end{aligned} A(q1)=1+a1q1+...+anaqnaB(q1)=b0+b1q1+...+bnbqnbC(q1)=1+c1q1+...+cncqnc(2)
设系统的输入为白噪声 ξ ( k ) \xi(k) ξ(k),输出为有色噪声 ε ( k ) \varepsilon(k) ε(k),这种线性系统通常称为形成滤波器,如式(1)。

把式(1)进行等效变换后,可以得到式(3):
A ( q − 1 ) C ( q − 1 ) y ( k ) = B ( q − 1 ) C ( q − 1 ) u ( k ) + ξ ( k ) (3) A(q^{-1})C(q^{-1})y(k)=B(q^{-1})C(q^{-1})u(k)+\xi(k) \tag 3 A(q1)C(q1)y(k)=B(q1)C(q1)u(k)+ξ(k)(3)
我们再把上式简化一下,令 { y ˉ ( k ) = C ( q − 1 ) y ( k ) u ˉ ( k ) = C ( q − 1 ) u ( k ) \begin{cases} \bar{y}(k)=C(q^{-1})y(k) \\ \bar{u}(k)=C(q^{-1})u(k) \end{cases} {yˉ(k)=C(q1)y(k)uˉ(k)=C(q1)u(k),就可以得到式(4):
A ( q − 1 ) y ˉ ( k ) = B ( q − 1 ) u ˉ ( k ) + ξ ( k ) (4) A(q^{-1})\bar{y}(k)=B(q^{-1})\bar{u}(k)+\xi(k) \tag 4 A(q1)yˉ(k)=B(q1)uˉ(k)+ξ(k)(4)

此时的噪声 ξ ( k ) \xi(k) ξ(k)已经是白噪声,可以直接使用LS法对系统参数进行无偏估计。

1.2.噪声模型

式(4)就是GLS的关键公式。但为了获得这个公式,我们必然先得到 C ( q − 1 ) C(q^{-1}) C(q1),这就需要我们对噪声模型进行估计。

根据系统模型表达式(1),可以得到关于噪声的差分方程描述:
ε ( k ) = 1 C ( q − 1 ) ξ ( k ) ⟹ C ( q − 1 ) ε ( k ) = ξ ( k ) ⟹ ε ( k ) = − c 1 ε ( k − 1 ) − . . . c n c ε ( k − n c ) + ξ ( k ) \begin{aligned} \varepsilon(k)=\frac{1}{C(q^{-1})}\xi(k) &\Longrightarrow C(q^{-1})\varepsilon(k)=\xi(k) \\ &\Longrightarrow \varepsilon(k)=-c_1\varepsilon(k-1)-...c_{n_c}\varepsilon(k-n_c)+\xi(k) \end{aligned} ε(k)=C(q1)1ξ(k)C(q1)ε(k)=ξ(k)ε(k)=c1ε(k1)...cncε(knc)+ξ(k)
由上式得到 N N N个方程,将其写成矩阵形式如下:
ε = Ω f + ξ \varepsilon=\Omega f + \xi ε=Ωf+ξ
根据最小二乘法,我们可以得到 C ( q − 1 ) C(q^{-1}) C(q1)的估计值 f f f
f = ( Ω T Ω ) − 1 Ω T ε f=(\Omega^T \Omega)^{-1}\Omega^T\varepsilon f=(ΩTΩ)1ΩTε
但这里存在一个问题:我们无法对噪声进行直接测量。也就是说,我们无法得到 Ω , ε \Omega,\varepsilon Ω,ε

一种方法是对误差进行估计,用残差 e ( k ) e(k) e(k)代替噪声 ε ( k ) \varepsilon(k) ε(k)
e ( k ) = A ^ ( q − 1 ) y ( k ) − B ^ ( q − 1 ) u ( k ) (5) e(k)=\hat{A}(q^{-1})y(k)-\hat{B}(q^{-1})u(k) \tag 5 e(k)=A^(q1)y(k)B^(q1)u(k)(5)
以上公式(3)、(4)、(5)就是推导GLS的主要公式了。

1.3.主要步骤

广义最小二乘法的基本思想:

  • 先不考虑有色噪声,用一般LS法估计系统参数 a i , b i a_i,b_i ai,bi,得到有偏估计 a ^ i ( 1 ) , b ^ i ( 1 ) \hat{a}^{(1)}_i,\hat{b}^{(1)}_i a^i(1),b^i(1)
  • 再用有偏估计 a ^ i ( 1 ) , b ^ i ( 1 ) \hat{a}^{(1)}_i,\hat{b}^{(1)}_i a^i(1),b^i(1)去估计噪声模型 c i c_i ci
  • 已知噪声模型后,根据等式(4)再去估计系统参数 a ^ i ( 2 ) , b ^ i ( 2 ) \hat{a}^{(2)}_i,\hat{b}^{(2)}_i a^i(2),b^i(2)

以上反复迭代后即可得到系统参数的无偏估计。

(1)用一般LS法得到系统参数的有偏估计

在GLS一开始,我们先不考虑噪声是白噪声还是有色噪声,直接当作白噪声用LS法进行估计,当然这次的估计值不会很准确,因为是一个有偏估计。但作为参数估计的初始值已经可以了,后面会再进行修正。
A ( q − 1 ) y ( k ) = B ( q − 1 ) u ( k ) + ε ( k ) A(q^{-1})y(k)=B(q^{-1})u(k)+\varepsilon(k) A(q1)y(k)=B(q1)u(k)+ε(k)
第一次参数估计的结果为
θ ^ ( 1 ) = ( Φ T Φ ) − 1 Φ T y \hat{\theta}^{(1)}=(\Phi^T\Phi)^{-1}\Phi^Ty θ^(1)=(ΦTΦ)1ΦTy

(2)用估计值 θ ^ \hat{\theta} θ^估计噪声模型参数 f f f

在已有系统参数估计值的情况下,通过式(5),就可以得到噪声 ε ( k ) \varepsilon(k) ε(k)的估计值 e ( k ) e(k) e(k)。根据 N N N个等式(5),我们可以得到

  • Ω \Omega Ω N × n c N \times n_c N×nc
  • e e e N × 1 N \times 1 N×1

然后使用LS法对噪声参数进行估计,结果为
f ^ ( 1 ) = ( Ω T Ω ) − 1 Ω T e \hat{f}^{(1)}=(\Omega^T \Omega)^{-1}\Omega^Te f^(1)=(ΩTΩ)1ΩTe

(3)用噪声估计参数 f ^ \hat{f} f^对系统估计参数 θ ^ \hat{\theta} θ^进行修正

得到噪声模型后,我们就可以得到 y ˉ ( k ) , u ˉ ( k ) \bar{y}(k),\bar{u}(k) yˉ(k),uˉ(k),然后应用式(4),通过最小二乘法对系统参数进行修正:
θ ^ ( 2 ) = ( Φ ˉ T Φ ˉ ) − 1 Φ ˉ T y ˉ ( k ) \hat{\theta}^{(2)}=(\bar{\Phi}^T\bar{\Phi})^{-1}\bar{\Phi}^T\bar{y}(k) θ^(2)=(ΦˉTΦˉ)1ΦˉTyˉ(k)
最后,重复步骤(2)~(3),直到估计值 θ ^ ( i ) \hat{\theta}^{(i)} θ^(i)收敛。

2.增广矩阵法(RELS)

2.1.系统模型与公式推导

考虑这样的一个系统(注意GLS和RELS的系统模型的差别):
A ( q − 1 ) y ( k ) = B ( q − 1 ) u ( k ) + C ( q − 1 ) ξ ( k ) (1) A(q^{-1})y(k)=B(q^{-1})u(k)+C(q^{-1})\xi(k) \tag 1 A(q1)y(k)=B(q1)u(k)+C(q1)ξ(k)(1)
上式中, ξ ( k ) \xi(k) ξ(k)是新息序列,具有白噪声特征。三个多项式形式如下所示:
A ( q − 1 ) = 1 + a 1 q − 1 + . . . + a n a q − n a B ( q − 1 ) = b 0 + b 1 q − 1 + . . . + b n b q − n b C ( q − 1 ) = 1 + c 1 q − 1 + . . . + c n c q − n c (2) \begin{aligned} A(q^{-1})=1+a_1q^{-1}+...+a_{n_a}q^{-{n_a}} \tag 2\\ B(q^{-1})=b_0+b_1q^{-1}+...+b_{n_b}q^{-{n_b}} \\ C(q^{-1})=1+c_1q^{-1}+...+c_{n_c}q^{-{n_c}} \end{aligned} A(q1)=1+a1q1+...+anaqnaB(q1)=b0+b1q1+...+bnbqnbC(q1)=1+c1q1+...+cncqnc(2)
增广矩阵法的特点是把噪声模型参数也加入到了估计参数中,相当于扩充被估计参数的维数,然后再用最小二乘法同时估计系统参数和噪声参数。
θ = [ a 1 . . . a n a b 0 . . . b n b c 1 . . . c n c ] T \theta=\begin{bmatrix} a_1 & ... & a_{n_a} & b_0 & ... & b_{n_b} & c_1 & ... & c_{n_c} \end{bmatrix}^T θ=[a1...anab0...bnbc1...cnc]T
根据式(1),我们可以得到以下形式:
y ( k ) = Ψ k T θ + ξ ( k ) ,    k = n + 1 , . . . , n + N (3) y(k)= \Psi_k^T \theta + \xi(k), \; k=n+1,...,n+N \tag 3 y(k)=ΨkTθ+ξ(k),k=n+1,...,n+N(3)

Ψ k = [ − y ( k − 1 ) . . . − y ( k − n a ) u ( k ) . . . u ( k − n b ) ξ ( k − 1 ) . . . ξ ( k − n c ) ] T \Psi_k = \begin{bmatrix} -y(k-1) & ... & -y(k-n_a) & u(k) & ... & u(k-n_b) & \xi(k-1) & ... & \xi(k-n_c) \end{bmatrix}^T Ψk=[y(k1)...y(kna)u(k)...u(knb)ξ(k1)...ξ(knc)]T

把上式扩展为 N N N个方程,就是一个标准的最小二乘法形式。但与广义最小二乘法遇到的问题一样,我们无法直接测量噪声数据,因此需要对 ξ ( k ) \xi(k) ξ(k)进行估计。 方法也与GLS一样,用残差 e ( k ) e(k) e(k)来估计噪声 ξ ( k ) \xi(k) ξ(k)
Ψ k = [ − y ( k − 1 ) . . . − y ( k − n a ) u ( k ) . . . u ( k − n b ) ξ ^ ( k − 1 ) . . . ξ ^ ( k − n c ) ] T \Psi_k = \begin{bmatrix} -y(k-1) & ... & -y(k-n_a) & u(k) & ... & u(k-n_b) & \hat{\xi}(k-1) & ... & \hat{\xi}(k-n_c) \end{bmatrix}^T Ψk=[y(k1)...y(kna)u(k)...u(knb)ξ^(k1)...ξ^(knc)]T
根据式(3),可以得到噪声的估计 ξ ^ ( k ) \hat{\xi}(k) ξ^(k)表示为:
ξ ^ ( k ) = y ( k ) − Ψ k T θ ^ (4) \hat{\xi}(k)=y(k)-\Psi_k^T \hat{\theta} \tag 4 ξ^(k)=y(k)ΨkTθ^(4)
根据递推最小二乘法的方式,我们可以得到增广矩阵法的递推方程:
K N + 1 = P N Ψ ^ N + 1 ( 1 + Ψ ^ N + 1 T P N Ψ ^ N + 1 ) − 1 P N + 1 = P N − K N + 1 Ψ ^ N + 1 T P N θ ^ N + 1 = θ ^ N + K N + 1 ( y N + 1 − Ψ ^ N + 1 T θ ^ N ) (5) \begin{aligned} K_{N+1}&=P_N \hat{\Psi}_{N+1}(1+\hat{\Psi}_{N+1}^TP_N\hat{\Psi}_{N+1})^{-1} \tag 5 \\ P_{N+1}&=P_N-K_{N+1}\hat{\Psi}_{N+1}^TP_N \\ \hat{\theta}_{N+1}&=\hat{\theta}_{N}+K_{N+1}(y_{N+1}-\hat{\Psi}_{N+1}^T \hat{\theta}_N) \end{aligned} KN+1PN+1θ^N+1=PNΨ^N+1(1+Ψ^N+1TPNΨ^N+1)1=PNKN+1Ψ^N+1TPN=θ^N+KN+1(yN+1Ψ^N+1Tθ^N)(5)

3.GLS vs RELS

由上面可知,GLS和RELS对于有色噪声的处理方式是不同的,下式(1)与(2)分别是GLS和RELS对有色噪声的处理方法, ξ ( k ) \xi(k) ξ(k)为白噪声。
A ( q − 1 ) y ( k ) = B ( q − 1 ) u ( k ) + ε ( k ) , ε ( k ) = 1 D ( q − 1 ) ξ ( k ) (1) A(q^{-1})y(k)=B(q^{-1})u(k)+\varepsilon(k),\varepsilon(k)=\frac{1}{D(q^{-1})}\xi(k) \tag 1 A(q1)y(k)=B(q1)u(k)+ε(k),ε(k)=D(q1)1ξ(k)(1)

A ( q − 1 ) y ( k ) = B ( q − 1 ) u ( k ) + ε ( k ) , ε ( k ) = C ( q − 1 ) ξ ( k ) (2) A(q^{-1})y(k)=B(q^{-1})u(k)+\varepsilon(k),\varepsilon(k)=C(q^{-1})\xi(k) \tag 2 A(q1)y(k)=B(q1)u(k)+ε(k),ε(k)=C(q1)ξ(k)(2)

虽然两者对噪声的表示形式不同,但这并不会影响最后的参数估计结果。 在进行参数辨识时,系统参数的估计值应该是一致的,区别在于对噪声模型的估计会不一样。

下面对一组数据分别用GLS、RELS进行参数辨识。

系统模型如下所示,参数真值为 a 1 = − 1.5 , a 2 = 0.7 , b 0 = 0 , b 1 = 1 , b 2 = 0.5 a_1=-1.5,a_2=0.7,b_0=0,b_1=1,b_2=0.5 a1=1.5,a2=0.7,b0=0,b1=1,b2=0.5
y ( k ) + a 1 y ( k − 1 ) + a 2 y ( k − 2 ) = b 0 u ( k ) + b 1 u ( k − 1 ) + b 2 u ( k − 2 ) + ε ( k ) y(k)+a_1y(k-1)+a_2y(k-2)=b_0u(k)+b_1u(k-1)+b_2u(k-2)+\varepsilon(k) y(k)+a1y(k1)+a2y(k2)=b0u(k)+b1u(k1)+b2u(k2)+ε(k)

ε ( k ) = ξ ( k ) + c 1 ξ ( k − 1 ) + c 2 ξ ( k − 2 ) \varepsilon(k)=\xi(k)+c_1\xi(k-1)+c_2\xi(k-2) ε(k)=ξ(k)+c1ξ(k1)+c2ξ(k2)

上式的噪声表示形式与RELS的表示形式相同,可以直接用RELS进行参数辨识。RELS的辨识结果为
θ ^ = [ − 1.5071 0.7006 − 0.0372 1.0410 0.4801 ] T , J = 40.6337 \hat{\theta}=\begin{bmatrix} -1.5071 & 0.7006 & -0.0372 & 1.0410 & 0.4801 \end{bmatrix}^T,J=40.6337 θ^=[1.50710.70060.03721.04100.4801]T,J=40.6337
在这里插入图片描述

对于GLS,噪声 ε ( k ) \varepsilon(k) ε(k)可以设为另外一种形式,即
ε ( k ) + f 1 ε ( k − 1 ) + f 2 ε ( k − 2 ) = ξ ( k ) ⟹ ε ( k ) = 1 1 + f 1 q − 1 + f 2 q − 2 ξ ( k ) \varepsilon(k)+f_1\varepsilon(k-1)+f_2\varepsilon(k-2)=\xi(k) \Longrightarrow \varepsilon(k)=\frac{1}{1+f_1q^{-1}+f_2q^{-2}}\xi(k) ε(k)+f1ε(k1)+f2ε(k2)=ξ(k)ε(k)=1+f1q1+f2q21ξ(k)
GLS的辨识结果为
θ ^ = [ − 1.5170 0.7115 − 0.0220 1.0441 0.4572 ] T , J = 40.7475 \hat{\theta}=\begin{bmatrix} -1.5170 & 0.7115 & -0.0220 & 1.0441 & 0.4572 \end{bmatrix}^T,J=40.7475 θ^=[1.51700.71150.02201.04410.4572]T,J=40.7475
在这里插入图片描述

由上可见,虽然两种方式对噪声的表示方式不同,但它们最终得到的系统参数的估计结果是相近,最终都收敛到了真实值。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值