Transcript Chapter 14: Statistical Least Square Estimation

Chapter 27: Linear Filtering -
Part I: Kalman Filter
Standard Kalman filtering – Linear
dynamics
Kalman Filter

Model dynamics – discrete time





xk+1 = M(xk) + Wk+1 or Mxk + wk+1
xk: true state
wk: model error
xk ∊ Rn, M: Rn  Rn, M ∊ Rn x n
Observation



zk = h(xk) + vk or Hxk + vk
h: Rn  Rm, H ∊ Rm x n , vk ∊ Rn
E(vk) = 0, COV(vk) = Rk
Filtering, Smoothing, Prediction
FN = {zi | 1 ≤ I ≤ N }
[Wiener, 1942]
[Kolmogorov, 1942]
k<N
k=N
k>N
Smoothing
filtering
Prediction
Go to page 464 – classification
Statement of Problem – Linear case





x0 ~ N(m0, P0)
xk+1 = Mkxk + wk+1, wk ~ N(0, Qk)
zk = Hkxk + vk, vk ~ N(0, Rk)
Given Fk = { zj | 1 ≤ j ≤ k }, find the best estimate
of
xk that minimizes the mean squared error
E[(xk - )T(xk - )] = tr[ E(xk - ) (xk - )T] = tr[ ]
If
is also unbiased => it is min. variance!
Model Forecast Step

At time k = 0, F0 – initial information is given



Given
, the predictable part of x1 is


Error in prediction

0
1
Forecast covariance

Covariance


Predicted Observation


E[z1/ x1= x1f] = E[H1x1+v1|x1=x1f] = H1x1f
COV(z1|x1f) = E{[z1 – E(z1|x1f)][z1 – E(z1|x1f)]T}
= E(v1v1T) = R1
Basic idea

At time k = 1
0

x1f
z1
P1f
R1
1
Fast Forward

Time k-1 to k
k-1
k
xkf
zk
Pkf
Rk
Forecast from k-1 to k


∴
Observations at time k

Model predicted observations


Data Assimilation Step

prior



Kalman Gain
Innovation
Posterior estimate – also known as analysis

Substitute and simplify



Covaiance of analysis

=>
where
Conditions on the Kalman gain –
minimization of total variance

Kk = PkfHkTDk-1 = PkfHkT[HkPkfHkT + Rk]-1
1.

2.

Comments on the Kalman gain

3. An Interpretation of Kk: n = m, Hk = I



Let Pkf = Diag [ P11f P22f … Pnnf]
Rk = Diag [ R11 R22 … Rnn]
Kk = PkfHkT[HkPkfHkT + Rk]-1
= Pkf[Pkf + Rk]-1
= Diag[ P11f/(P11f+R11) , P22f/(P22f+R22), ….., Pnnf/(Pnnf+Rnn)]
= xkf + Kk[z-Hxkf] = xkf + Kk[z-xkf] = (I-Kk)xkf + Kkz


∴ If Piif is large, zi,k has a larger weight
Comments – special cases


4.
is independent of observations
5. No observations




xkf =
Pkf =
for all k ≥ 0
xkf = Mk-1xk-1f = Mk-1Mk-2Mk-3…..M1M0x0f
Pkf = Mk-1Pk-1fMk-1T + Qk
= M(k-1:0)P0MT(k-1:0) + ∑ M(k-1:j+1)Qj+1MT(k-1:j+1)
where M(i,j) = MiMi-1Mi-2….Mj, Qj ≡ 0
Pk+1f = M(k-1:0)P0MT(k-1:0)
Special cases - continued

6. No Dynamics

Mk =I, Wk ≡ 0, Qk ≡ 0
xk+1 = xk = x
zk = Hkx + vk
xkf =
with
= E(x0)
Pkf =
with
= P0
=>

Same as (17.2.11) – (17.2.12) Static case





Special cases

7. When observations are perfect


Rk ≡ 0
=> Kk = PkfHkT[HkPkfHkT + Rk]-1
= PkfHkT[HkPkfHkT]-1


Hk: m x n, Pkf: n x n, HkT: n x m
=> [HkPkfHkT]-1: m x m

Recall:

From (27.2.19):
Special cases






∴ ( I- KkHk ) = ( I – KkHk )2, idempotent
Fact: Idempotent matrices are singular
=> Rank of ( I- KkHk ) ≤ n – 1
∴ Rank(
) ≤ min {Rank( I- KkHk ), Rank( Pkf )} ≤ n – 1
∴ Rank of
≤n–1
∴ When Rk is small, this will cause computational
instability
Special cases

8. Residual Checking





rk = zk – Hkxkf = innovation
= xkf + Kkrk
rk = zk –Hkxkf
= Hkxk + vk – Hkxkf
= Hk(xk – xkf) + vk
= Hkekf + vk
∴ COV(rk) = HkPkfHkT + Rk
∴ By computing rk and its covariance, we can check if the filter is
working O.K.

10. Computational Cost
Example 27.2.1 Scalar Dynamics with No
Observation


a > 0, wk ~ N( 0, q), x0 ~ N(m0, P0)
xk = axk-1 + wk




E(xk) = akE(x0) = akm0
Pk = Var(xk) = Var(axk-1 +wk) = a2Pk-1 + q
∴ Pk = a2kP0 + q[(a2k -1)/(a2 – 1)]
Scalar dynamics


Note: For a given m0, P0, q, the behavior of the moments
depends on a
1. 0 < a < 1


2. 1 < a < ∞


lim E(xk)  0, lim Pk = q/(1-a2)
lim E(xk) = 0, lim Pk = ∞
3. a = 1



xk = x0 + ∑wk
E(xk) = m0
Pk = P0 + kq
Example 27.2.2 Kalman Filtering







xk+1 = axk + wk+1, wk+1~(0,q)
zk = hxk + vk, vk~N(0,r)
xk+1f = a
Pk+1f = a2 + q
= xkf + Kk[zk – hxkf]
Kk = Pkfh[h2Pkf + r]-1 = hr-1
= Pkf – (Pkf)2h2[h2Pkf + r]-1 = Pkfr[h2Pkf+r]-1
Recurrences: Analysis of Stability


HW. 1
xk+1f = a(1-Kkh)xkf + aKkzk


ek+1f = a(1-Kkh)ekf + aKkvk + wk+1



Pk+1f = a2 +q
= Pkfr(h2Pkf+r)-1
Example continued



Pk+1f = a2Pkfr / (h2Pkf+r) + q
Pk+1f / r = a2Pkf / (h2Pkf + r) + q/r
= a2(Pkf/r) / [h2(Pkf/r) + 1] + q/r
Pk+1 = a2Pk/(h2Pk+1) + α



α = q/r (ratio)
Pk = Pkf/r
Riccati equation ( First-order, scalar, nonlinear)
Asymptotic Properties




Let h = 1
Pk+1 = a2Pk/(Pk+1) + α
Let δk = Pk+1 – Pk
= a2Pk/(Pk+1) – Pk + α
= [-Pk2 + Pk(a2- 1 + α) + α]/(Pk+1)
∴ δk = g(Pk)/(Pk+1)

where g(Pk) = -Pk2 + Pk(a2 + α - 1) + α
Example continued




When Pk+1 = Pk, => δk = 0 => equilibrium
=> δk = 0 if g(Pk) = 0
β
-Pk2 + Pk(a2 + α – 1) + α = 0
-Pk2 + β Pk + α = 0




Evaluate the derivative of g(.) at P* and P*
g’(Pk) = -2Pk + β
Example continued



∴ P* is an attractor – stable
P* is a repellor – unstable
Thus,
=> P* = a2P*/(P*+1) + α
Rate of Convergence



Let yk = pk – p* , pk+1 = a2pk/(pk+1) + α
yk+1 = pk+1 – p*
= [a2pk/(pk+1) + α] - [a2p*/(p*+1) + α]
= a2pk/(pk+1) - a2p*/(p*+1)
= a2(pk – p*)/[(1+pk)(1+p*)]
= a2yk/[(1+pk)(1+p*)]
∴ 1/yk+1 = [(1+pk)(1+p*)]/ a2yk
= [(1+yk+p*)(1+p*)] / a2yk
= [(1+p*)/a]2/yk + (1+p*)/a2
Rate convergence - continued

zk = 1/yk
=> zk+1 = czk + b
where c = [(1+p*)/a]2 and b = (1+p*)/a2
Iterating:

∴

when c> 1 (ie) c = [(1+p*)/a]2 >1
When this is true, yk  0 at exp. rate


Rate of convergence - continued

From (27.2.38) h = 1

Analysis Covariance Converges
Stability of the Filter

h=1
ek+1f = a(1-Kkh)ekf + aKkvk + wk+1
-- homo part
Kk = pkf/(pkf+r) = pk/(pk+1)
1 – Kk = 1/(pk+1)

∴

∴



