最优估计 – kalman and lsm
kalman Filter 和 least square 目的均为最优化某一指标,指标是优化的关键:
常用的估计准则有:
-
无偏估计:估计值的期望等于被估计参数的真实值。
-
线性最小方差估计:将估计量限制为观测值的线性函数,已知观测量Z和和被估计量X一二阶矩(EX,Var{X},EZ,Var{Z},Cov{X,Z}),使估计误差的方差最小,即最小化 t r { E [ X ~ − E X ~ ] [ X ~ − E X ~ ] T } tr\{E[\tilde{X}-E\tilde{X}][\tilde{X}-E\tilde{X}]^{T}\} tr{E[X~−EX~][X~−EX~]T} , X ~ \tilde{X} X~为估计误差(等价于最小化均方误差阵,若为无偏估计)可得其无偏估计值为 X ~ L M V ( Z ) = E X + c o v ( X , Z ) ( v a r ( Z ) ) − 1 [ Z − E Z ] \tilde{X}_{LMV}(Z)=EX+cov(X,Z)(var(Z))^{-1}[Z-EZ] X~LMV(Z)=EX+cov(X,Z)(var(Z))−1[Z−EZ]对于观测模型Z=HX+V,上述条件若已知
{ E X = μ x , V a r ( X ) = P x , E V = 0 , V a r ( V ) = R , E ( X V T ) = 0 } \{EX=\mu_x,Var(X)=P_x,EV=0,Var(V)=R,E(XV^T)=0\} {EX=μx,Var(X)=Px,EV=0,Var(V)=R,E(XVT)=0} 即可得到。
-
最小二乘估计:对数据(X、Z)的统计特性一无所知,但仍需对X进行估计,目标是最小化残差1平方和。
满足最小方差必满足残差平方和最小,反之则不成立。
经典最小二乘
针对隐状态X,若其无法直接观测,但间接获取其观测值
Z
=
[
z
1
,
z
2
,
…
,
z
n
]
T
Z=[z_1,z_2,\dots,z_n]^T
Z=[z1,z2,…,zn]T ,若其观测值为状态值的线性函数:
Z
i
=
H
i
X
+
V
i
,
i
=
1
,
…
,
n
Z_i=H_iX+V_i,i=1,\dots,n
Zi=HiX+Vi,i=1,…,n
z
i
z_i
zi为第i次测量的观测值,
H
i
H_i
Hi为第i次测量的观测模型(设计矩阵,实验的观测值),
V
i
V_i
Vi为第i次测量的噪声(误差)。
则第i次测量的估计误差:
e
i
^
=
z
i
−
H
i
X
^
\hat{e_i}=z_i-H_i\hat{X}
ei^=zi−HiX^
则n次测量的误差(残差)平方和为优化指标:
J
(
X
^
)
=
∑
i
=
1
n
(
z
i
−
H
i
X
^
)
2
=
(
Z
−
H
X
^
)
T
(
Z
−
H
X
^
)
=
t
r
[
(
Z
−
H
X
^
)
(
Z
−
H
X
^
)
T
]
J(\hat{X})=\sum_{i=1}^{n}{(z_i-H_i\hat{X})^2}=(Z-H\hat{X})^T(Z-H\hat{X}) \\ =tr[(Z-H\hat{X})(Z-H\hat{X})^T]
J(X^)=i=1∑n(zi−HiX^)2=(Z−HX^)T(Z−HX^)=tr[(Z−HX^)(Z−HX^)T]
令
∂
J
∂
X
^
=
0
\frac{\partial{J}}{\partial{\hat{X}}}=0
∂X^∂J=0 ,可得最小二乘估计值:
X
^
L
S
=
(
H
T
H
)
−
1
H
T
Z
\hat{X}_{LS}=(H^TH)^{-1}H^TZ
X^LS=(HTH)−1HTZ
将
Z
=
H
X
+
V
Z=HX+V
Z=HX+V此时状态的估计误差:
X
~
L
S
=
X
−
X
^
L
S
=
−
(
H
T
H
)
−
1
H
T
V
\tilde{X}_{LS}=X-\hat{X}_{LS}=-(H^TH)^{-1}H^TV
X~LS=X−X^LS=−(HTH)−1HTV
若测量噪声均值为0,则
E
(
X
~
L
S
)
=
0
E(\tilde{X}_{LS})=0
E(X~LS)=0,此时最小二乘估计为**无偏估计**,状态估计误差的(协)方差[^ 2]
V
a
r
(
X
~
L
S
)
=
E
[
(
X
~
−
E
X
~
)
(
X
~
−
E
X
~
)
T
]
Var(\tilde{X}_{LS})=E[(\tilde{X}-E\tilde{X})(\tilde{X}-E\tilde{X})^T]
Var(X~LS)=E[(X~−EX~)(X~−EX~)T]与估计量的均方误差矩阵
E
[
X
−
X
^
]
[
X
−
X
^
]
T
E[X-\hat{X}][X-\hat{X}]^T
E[X−X^][X−X^]T相等。可见标准最小二乘不需要噪声V的任何统计信息。
由(5)式可得:
V
a
r
(
X
~
L
S
)
=
E
[
X
−
X
^
]
[
X
−
X
^
]
T
=
(
H
T
H
)
−
1
H
T
E
(
V
V
T
)
H
(
H
T
H
)
−
1
=
(
H
T
H
)
−
1
H
T
R
H
(
H
T
H
)
−
1
\begin{aligned} Var(\tilde{X}_{LS})=E[X-\hat{X}][X-\hat{X}]^T & = (H^TH)^{-1}H^TE(VV^T)H(H^TH)^{-1}\\ &=(H^TH)^{-1}H^TRH(H^TH)^{-1} \end{aligned}
Var(X~LS)=E[X−X^][X−X^]T=(HTH)−1HTE(VVT)H(HTH)−1=(HTH)−1HTRH(HTH)−1
其中
R
=
E
(
V
V
T
)
R=E(VV^T)
R=E(VVT)为测量误差(噪声)的(协)方差阵。
加权最小二乘(weighted least square)
在经典最小二乘中,假定每一次测量的权重相同,但是一般来说近期数据比远期数据影响更大,因此引入加权最小二乘,其指标形式:
J
W
(
X
^
)
=
∑
i
=
1
n
(
z
i
−
H
i
X
^
)
2
=
(
Z
−
H
X
^
)
T
W
(
Z
−
H
X
^
)
J_W(\hat{X})=\sum_{i=1}^{n}{(z_i-H_i\hat{X})^2}=(Z-H\hat{X})^TW(Z-H\hat{X})
JW(X^)=i=1∑n(zi−HiX^)2=(Z−HX^)TW(Z−HX^)
同样使其偏导数为0,可得
X
^
L
S
W
=
(
H
T
W
H
)
−
1
H
T
W
Z
\hat{X}_{LSW}=(H^TWH)^{-1}H^TWZ
X^LSW=(HTWH)−1HTWZ
由附录2,若噪声不满足同方差,则普通最小二乘(4)并不是BLUE,此时噪声的协方差阵
E [ V V T ] = σ 2 R , R ≠ I E[VV^T]=\sigma^2R,R\neq{I} E[VVT]=σ2R,R=I , R = [ r 1 ⋱ r n ] R=\begin{bmatrix}r_1\\&\ddots\\&& r_n\end{bmatrix} R=⎣⎡r1⋱rn⎦⎤,即原模型存在异方差性。
设
R
=
D
D
T
,
D
=
[
r
1
⋱
r
n
]
R=DD^T,D=\begin{bmatrix}\sqrt{r_1}\\&\ddots\\&& \sqrt{r_n}\end{bmatrix}
R=DDT,D=⎣⎡r1⋱rn⎦⎤ ,用
D
−
1
D^{-1}
D−1同时左乘
Z
=
H
X
+
V
Z=HX+V
Z=HX+V两端得到新的模型:
D
−
1
Z
=
D
−
1
H
X
+
D
−
1
V
Z
⋆
=
H
⋆
X
+
V
⋆
\begin{aligned} D^{-1}Z&=D^{-1}HX+D^{-1}V \\ Z^{\star}&=H^{\star}X+V^{\star} \end{aligned}
D−1ZZ⋆=D−1HX+D−1V=H⋆X+V⋆
此时,原模型的加权最小二乘估计量为无偏的。
E
[
V
⋆
V
⋆
T
]
=
E
[
D
−
1
V
V
T
D
−
1
T
]
=
D
−
1
E
[
V
V
T
]
D
−
1
T
=
σ
2
D
−
1
R
D
−
1
T
=
σ
2
I
\begin{aligned} E[V^{\star}V^{\star T}]&=E[D^{-1}VV^TD^{-1\ T}]\\ &=D^{-1}E[VV^T]D^{-1\ T}\\ &=\sigma^2D^{-1}RD^{-1\ T}\\ &=\sigma^2I \end{aligned}
E[V⋆V⋆T]=E[D−1VVTD−1 T]=D−1E[VVT]D−1 T=σ2D−1RD−1 T=σ2I
此时得到的参数估计为:
X
^
L
S
W
=
(
H
⋆
T
H
⋆
)
−
1
H
⋆
T
Z
⋆
=
(
H
T
R
−
1
H
)
−
1
H
T
R
−
1
Z
\begin{aligned} \hat{X}_{LSW}&=(H^{\star T}H^{\star})^{-1}H^{\star T}Z^{\star}\\ &=(H^TR^{-1}H)^{-1}H^TR^{-1}Z \end{aligned}
X^LSW=(H⋆TH⋆)−1H⋆TZ⋆=(HTR−1H)−1HTR−1Z
可以证明(见附录),当
W
=
R
−
1
W=R^{-1}
W=R−1时,最小二乘估计时缺少初值条件下的线性无偏最小方差估计(BLUE,Best Linear Unbiased Estimation)——即能够使估计误差的方差阵最小,又称马尔可夫估计,其中
R
=
E
[
V
V
T
]
R=E[VV^T]
R=E[VVT]
为随机噪声的(协)方差阵(对称正定阵)。
递推最小二乘(Recursive Least Square,RLS)
上述方法进行一次估计需要所有历史数据,不利于在线估计,考虑前n次测量:
Z
n
=
H
n
X
+
V
n
Z_n=H_nX+V_n
Zn=HnX+Vn
则加权的最小二乘估计为:
X
^
L
S
W
(
n
)
=
(
H
n
T
R
n
−
1
H
n
)
−
1
H
n
T
R
n
−
1
Z
n
\hat{X}_{LSW}(n)=(H_{n}^TR_{n}^{-1}H_n)^{-1}H_{n}^TR_{n}^{-1}Z_n
X^LSW(n)=(HnTRn−1Hn)−1HnTRn−1Zn
估计误差的(协)方差矩阵为:
P
n
=
E
[
X
~
L
S
W
(
n
)
X
~
L
S
W
T
(
n
)
]
=
E
[
−
(
H
T
R
−
1
H
)
−
1
]
H
T
R
−
1
V
V
T
R
−
1
H
(
H
T
R
−
1
H
)
−
1
=
(
H
T
R
−
1
H
)
−
1
H
T
R
−
1
H
(
H
T
R
−
1
H
)
−
1
=
(
H
T
R
−
1
H
)
−
1
\begin{aligned} P_n&=E[\tilde{X}_{LSW}(n)\tilde{X}_{LSW}^T(n)]\\ &=E[-(H^TR^{-1}H)^{-1}]H^TR^{-1}VV^TR^{-1}H(H^TR^{-1}H)^{-1}\\ &=(H^TR^{-1}H)^{-1}H^TR^{-1}H(H^TR^{-1}H)^{-1}\\ &=(H^TR^{-1}H)^{-1} \end{aligned}
Pn=E[X~LSW(n)X~LSWT(n)]=E[−(HTR−1H)−1]HTR−1VVTR−1H(HTR−1H)−1=(HTR−1H)−1HTR−1H(HTR−1H)−1=(HTR−1H)−1
结合上述两式,可得:
X
^
L
S
W
(
n
)
=
P
n
H
n
T
R
n
−
1
Z
n
\hat{X}_{LSW}(n)=P_nH_{n}^TR_{n}^{-1}Z_n
X^LSW(n)=PnHnTRn−1Zn
现得到一个新的测量值:
z
n
+
1
=
H
n
+
1
X
+
v
n
+
1
z_{n+1}=H_{n+1}X+v_{n+1}
zn+1=Hn+1X+vn+1
添加到矩阵中:
X
^
L
S
W
(
n
+
1
)
=
(
H
n
+
1
T
R
n
+
1
−
1
H
n
+
1
)
−
1
H
n
+
1
T
R
n
+
1
−
1
Z
n
+
1
\hat{X}_{LSW}(n+1)=(H_{n+1}^TR_{n+1}^{-1}H_{n+1})^{-1}H_{n+1}^TR_{n+1}^{-1}Z_{n+1}
X^LSW(n+1)=(Hn+1TRn+1−1Hn+1)−1Hn+1TRn+1−1Zn+1
将新的测量噪声加入到原本的测量噪声矩阵中:R阵应为对角阵:
R
k
+
1
−
1
=
[
R
n
−
1
0
0
r
n
+
1
−
1
]
R_{k+1}^{-1}= \begin{bmatrix} R_n^{-1} & 0 \\0&r^{-1}_{n+1} \end{bmatrix}
Rk+1−1=[Rn−100rn+1−1]
将式子展开:
P
n
+
1
−
1
=
H
n
+
1
T
R
n
+
1
−
1
H
n
+
1
=
[
H
n
T
,
h
n
+
1
T
]
[
R
n
−
1
0
0
r
n
+
1
−
1
]
[
H
n
h
n
+
1
]
=
H
n
T
R
n
−
1
H
n
+
h
n
+
1
T
r
n
+
1
−
1
h
n
+
1
P_{n+1}^{-1}=H_{n+1}^TR_{n+1}^{-1}H_{n+1}=[H_n^T,h_{n+1}^T] \begin{bmatrix} R_n^{-1} & 0 \\0&r^{-1}_{n+1} \end{bmatrix} \begin{bmatrix} H_n\\h_{n+1} \end{bmatrix} =H_n^TR_n^{-1}H_n+h_{n+1}^Tr_{n+1}^{-1}h_{n+1}
Pn+1−1=Hn+1TRn+1−1Hn+1=[HnT,hn+1T][Rn−100rn+1−1][Hnhn+1]=HnTRn−1Hn+hn+1Trn+1−1hn+1
即:
P
n
+
1
−
1
=
P
n
−
1
+
h
n
+
1
T
r
n
+
1
−
1
h
n
+
1
P_{n+1}^{-1}=P_n^{-1}+h_{n+1}^Tr_{n+1}^{-1}h_{n+1}
Pn+1−1=Pn−1+hn+1Trn+1−1hn+1
综上,可以推得:
P
n
+
1
=
P
n
−
P
n
h
n
+
1
T
[
h
n
+
1
P
n
h
n
+
1
T
+
r
n
+
1
]
−
1
h
n
+
1
P
n
K
n
+
1
=
P
n
+
1
h
n
+
1
T
r
n
+
1
−
1
X
^
L
S
W
(
n
+
1
)
=
X
^
L
S
W
(
n
)
+
K
n
+
1
[
z
n
+
1
−
h
n
+
1
X
^
L
S
W
(
n
)
]
\begin{aligned} P_{n+1}&=P_n-P_nh_{n+1}^T[h_{n+1}P_nh_{n+1}^T+r_{n+1}]^{-1}h_{n+1}P_n\\ K_{n+1} &= P_{n+1}h_{n+1}^Tr_{n+1}^{-1}\\ \hat{X}_{LSW}(n+1)&=\hat{X}_{LSW}(n)+K_{n+1}[z_{n+1}-h_{n+1}\hat{X}_{LSW}(n)] \end{aligned}
Pn+1Kn+1X^LSW(n+1)=Pn−Pnhn+1T[hn+1Pnhn+1T+rn+1]−1hn+1Pn=Pn+1hn+1Trn+1−1=X^LSW(n)+Kn+1[zn+1−hn+1X^LSW(n)]
其中
K
n
+
1
K_{n+1}
Kn+1可将(31)代入展开为:
K
n
+
1
=
P
n
h
n
+
1
T
[
h
n
+
1
P
n
h
n
+
1
T
+
r
n
+
1
]
−
1
K_{n+1} = P_nh_{n+1}^T[h_{n+1}P_nh_{n+1}^T+r_{n+1}]^{-1}
Kn+1=Pnhn+1T[hn+1Pnhn+1T+rn+1]−1
因此
P
n
+
1
P_{n+1}
Pn+1亦可表示为:
P
n
+
1
=
P
n
−
K
n
+
1
h
n
+
1
P
n
P_{n+1}=P_n-K_{n+1}h_{n+1}P_n
Pn+1=Pn−Kn+1hn+1Pn
卡尔曼滤波
若被估计量X不随时间变化,或随时间缓慢变化则为“静态估计”,而被估计量随时间变化为“动态估计”。
(待填)
参考:
https://blog.csdn.net/qinruiyan/article/details/50793114
《最优估计理论》刘胜,张红梅,科学出版社
最佳线性无偏估计(GM假设)
假设多元线性回归模型: Z = H X + V Z=HX+V Z=HX+V
Z = ( z 1 , … , z n ) T H = [ h i j ] n × p X = ( x o , … , x p ) V = ( v 0 , … , v n ) \begin{aligned}Z&=(z_1,\dots,z_n)^T\\H&=\begin{bmatrix}h_{ij}\end{bmatrix}_{n\times{p}}\\X&=(x_o,\dots,x_p)\\V&=(v_0,\dots,v_n)\end{aligned} ZHXV=(z1,…,zn)T=[hij]n×p=(xo,…,xp)=(v0,…,vn)
则GM假设:
E ( V ∣ H ) = 0 , ∀ H ( 零 均 值 ) V a r ( V ∣ H ) = E ( V V T ∣ H ) = σ 2 I n ( 同 方 差 且 不 相 关 ) \begin{aligned}E(V|H)&=0,\forall H\ (零均值)\\Var(V|H)&=E(VV^T|H)=\sigma^2I_n\ (同方差且不相关)\end{aligned} E(V∣H)Var(V∣H)=0,∀H (零均值)=E(VVT∣H)=σ2In (同方差且不相关)
则此时对参数X的最佳线性无偏估计为:
X ^ = ( H T H ) − 1 H T Z \hat{X}=(H^TH)^{-1}H^TZ X^=(HTH)−1HTZ最小二乘估计与最小方差估计等价条件证明:
各种估计方法的比较:
[外链图片转存中…(img-PDFlGI9x-1576551483788)]
残差在数理统计中是指实际观察值和估计值之间的差。若设线性回归模型为 Z = H X + V Z=HX+V Z=HX+V ,其中Z为n维输出向量,H是 n × ( p + 1 ) n\times(p+1) n×(p+1) 阶设计矩阵,X是p+1维向量,V为n维随机变量(扰动)。则回归系数的估计值 X ^ = ( H T H ) − 1 H T Z \hat{X}=(H^TH)^{-1}H^TZ X^=(HTH)−1HTZ ,拟合值 Z ^ = H X ^ = H ( H T H ) − 1 H T Z \hat{Z} = H\hat{X}=H(H^TH)^{-1}H^TZ Z^=HX^=H(HTH)−1HTZ,残差为 ϵ ^ = z i − z i ^ = z i − H i X ^ \hat{\epsilon}=z_i-\hat{z_i}=z_i-H_i\hat{X} ϵ^=zi−zi^=zi−HiX^ ,其由观测真值和H阵给出,不考虑噪声V。
[^ 2]: https://zh.wikipedia.org/wiki/%E5%8D%8F%E6%96%B9%E5%B7%AE%E7%9F%A9%E9%98%B5 ↩︎在线性回归模型中,如果随机噪声(误差)满足零均值、同方差且互不相关,则回归系数的最优线性无偏估计(BLUE,Best Linear unbiased estimator)就是普通最小二乘估计。 ↩︎