4. Optimální odhad stavů systému
algoritmy + příklady
Hana Fitzová
Brno, 2008
Obsah
Kalmanův filtr
Smoothing
Příklady ­ odhad polohy lodi, mezery výstupu HDP,
reálné úrokové míry a reálného kurzu
Rozšířený Kalmanův filtr
Odhad nelineárních systémů metodou Monte Carlo
Vážený Bootstrap algoritmus
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.1/3
Kalmanův filtr
Kalmanův filtr je analogií výše uvedené nelineární
rekurze pro speciální případ, kdy jsou stavy x i
výstupy y rozloženy normálně.
Pak je apriorní i aposteriorní hustota opět normální.
Normálně rozloženou náhodnou veličinu plně určuje
její střední hodnota a rozptyl.
V průběhu výpočtu tedy stačí sledovat pouze tyto
dvě číselné charakteristiky, čímž se celý výpočet
zjednodušuje.
Mají-li být stavy rozloženy normálně, musí být
příslušný dynamický systém lineární s gaussovskými
šumy, což je největší slabina tohoto algoritmu.
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.2/3
Věta (Kalmanův filtr I)
Mějme diskrétní stochastický systém tvaru:
xt+1 = Axt + But + vt(1)
yt = Cxt + Dut + wt(2)
splňující: x0  N(0, 0), vt  N(0, v), wt  N(0, w),
EvtwT
t = 0, Ex0vT
t = 0 t, kde vektor 0 a matice
0, v, w jsou známé.
Nechť jsou známa data Dt. Označme xt|k = xt|Dk.
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.3/3
Věta (Kalmanův filtr II)
Pak apriorní estimátor xt|t-1 a aposteriorní estimátor xt|t
jsou normální pro všechna t:
xt|t-1  N(t|t-1, t|t-1)
xt|t  N(t|t, t|t)
kde střední hodnoty t|t-1, t|t a varianční matice
t|t-1, t|t se vypočtou dle následujících rekurentních
vztahů:
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.4/3
Věta (Kalmanův filtr III)
t|t = t|t-1 + Kt(yt - Ct|t-1 - Dut)(3)
t|t = t|t-1 - KtCt|t-1(4)
t+1|t = At|t + But(5)
t+1|t = At|tAT
+ v(6)
Kt = t|t-1CT
(Ct|t-1CT
+ w)-1
(7)
První dva vztahy tvoří tzv. datový krok algoritmu, druhé
dva vztahy krok predikční, poslední vztah je tzv.
Kalmanovo zesílení (Kalmanův zisk).
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.5/3
Smoothing v KF
Zpětný běh filtru na základě informací z času t = 1, . . . , N
t|N = t|t + Ft(t+1|N - t+1|t)(8)
t|N = t|t - Ft(t+1|t - t+1|N )Ft
T
(9)
Ft = t|tAT
-1
t+1|t
(10)
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.6/3
Příklad ­ loď (1)
Loď se pohybuje po rovníku východním směrem
rychlostí 10 námořních mil za hodinu. Okamžitou
rychlost lodi ovlivňují náhodné poryvy větru a nárazy
vln.
Navigátor lodi odhaduje každou hodinu zeměpisnou
délku lodi l a rychlost lodi s = dl/dt v mph. V čase
t = 0 odhadl navigátor polohu l0 = 0 a rychlost s = 10.
Dále pak zaznamenal údaje:
čas 1 2 3 4 5 6
poloha 9" 19.5" 29" 38.4" 50" 59.5"
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.7/3
Příklad ­ loď (2)
Označíme-li lt, st polohu a rychlost lodi v čase t, pak
je úlohou navigátora optimální odhad veličin lt, st.
Počáteční odhady můžeme modelovat jako nezávislé
náhodné veličiny s normálním rozložením. Rozptyly
odhadů jsou sledovány a jejich odhady jsou Dl0 = 2,
Ds0 = 3.
Nejprve popíšeme chování systému. Během hodiny t
se loď pohybuje rychlostí st, takže její poloha se
změní na:
lt+1 = lt + st
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.8/3
Příklad ­ loď (3)
Rychlost kolísá díky náhodným vlivům, což
popíšeme pomocí náhodné veličiny et  N(0, 1):
st+1 = st + et
Rozptyl měrení sextantu udávaný výrobcem je
w = 2.
Stavový vektor definujeme jako xt =
lt
st
S využitím vztahů (9)-(13) odhadněte optimlání stavy
tohoto systému pomocí Kalmanova filtru.
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.9/3
Návod (loď)
Načtěte matice systému a data (A, B, C, D, v, w,
x0, x0 , U, y, . . . ) 0|-1 = x0, 0|-1 = x0
V cyklu spočtěte t|t, t|t, t+1|t, t+1|t, t|N , t|N
Výsledky vykreslete do obrázků: predikce + filtrace +
smoothing pro polohu, totéž pro rychlost; dále vývoj
rozptylu (f+p+s) obou veličin.
Vytvořte funkci, která bude provádět jednotlivé kroky
Kalmanova filtru včetně smoothingu a příklad upravte
tak, aby tuto funkci využíval. Vstupními parametry
funkce bude: y, U, A, B, C, D, v, w, x0 a x0 .
Výstupy budou: t|t, t|t, t+1|t, t+1|t, t|N , t|N .
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.10/3
Příklad ­ mezera výstupu (1)
Příklad na odhad české mezery výstupu (odhad
hospodářského cyklu)
Předpokládáme, že mezera výstupu je bíly šum
Předpokládáme, že potenciální produkt je náhodná
procházka s tempem růstu
Budeme pracovat s logaritmy (rozdíl logaritmů je
zhruba tempo růstu)
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.11/3
Příklad ­ mezera výstupu (2)
Y je produkt
Y EQ je potenciální produkt
Y GAP je mezera výstupu
y = log(Y )
v1(t), v2(t) jsou šumy
gr, grss jsou aktuální a steady-state tempo růstu
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.12/3
Základní model (ZM)
yEQ
(t) = yEQ
(t - 1) + gr + v1(t)
yGAP
(t) = v2(t)
y(t) = yEQ
(t) + yGAP
(t)
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.13/3
Počáteční nastavení (ZM)
Potenciál je mnohem rigidnější než mezera výstupu,
proto zvolíme v = [0.01 0; 0 1]
Výstupní rovnice je deficitní, tj. w = 0
Vyjdeme z toho, že počátkem 90. let jsme byli na
potenciálu, tj. x0 = [y(1); 0]
Nejistota o volbě počáteční podmínky je veliká:
x0 = 10  v
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.14/3
Rozšířený model (RM)
yEQ
(t) = yEQ
(t - 1) + gr(t) + v1(t)
yGAP
(t) = v2(t)
gr(t) = 0.6gr(t - 1) + 0.4grss
y(t) = yEQ
(t) + yGAP
(t)
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.15/3
Počáteční nastavení (RM)
Další stavovou proměnnou je tempo růstu, je
persistentní (autoregrese), blíží se steady-state
růstu, který je aproximován průměrným tempem
růstu HDP.
Obdobně jako výše zvolíme:
v = [0.0001 0 0; 0 1 0; 0 0 0.002]
w = 0
x0 = [y(1); 0; 2]
x0 = 100  v
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.16/3
Návod
Převedťe systém do stavového tvaru a matice
systému dopište do předchystaného souboru
mezera.m
Data k příkladu naleznete v souboru czech_data.m
Provedťe odhad stavů systému pomocí Kalmanova
filtru, získané výsledky (odhadnutý vývoj
potenciálního produktu a mezery výstupu) vykreslete
do obrázků a okomentujte je
Totéž provedťe pro rozšířený model a sledujte vliv
změn nastavení rovnovážného tempa růstu a
počátečních podmínek na dosažené výsledky
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.17/3
Příklad ­ transmisní kanály
Příklad na odhad rovnovážné úrovně české a
německé reálné úrokové míry a rovnovážné úrovně
reálného kurzu CZ/GER
Data k příkladu naleznete v souboru cz_data.m
Odhad neznámých stavů se nalézá v souboru
kanaly.m
Jedná se o interdependentní systém se 4 stavovými
rovnicemi, je nutno jej nejprve převést
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.18/3
Označení
Jedná se o čtvrtletní údaje za období Q1-1995 až
Q4-2004
Předpona gr v názvu označuje německá data
Přípona eq v názvu označuje rovnovážnou úroveň
veličiny
rs4, rr4 je nominální resp. reálná úroková míra
CPI, pie je index CPI, resp. inflace CPI
prem je riziková prémie
ls, lz je logaritmus nominálního resp. reálného
směnného kurzu CZ/GER
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.19/3
Rozšířený Kalmanův filtr I
Princip rozšířeného Kalmanova filtru spočívá v tom,
že nelineární dynamický systém v každém kroku
približně nahradíme lineárním dynamickým
systémem a na tento pak aplikujeme obyčejný
Kalmanův filtr.
Vztahy odvozené výše nelze užít přímo, poněvadž
linearizace vede na jiný typ systému (viz dále).
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.20/3
Rozšířený Kalmanův filtr II
Mějme nelineární diskrétní stochastický systém:
xt+1 = f(xt, ut) + vt(11)
yt = g(xt, ut) + wt(12)
s počáteční podmínkou x0  N(0, 0).
Rozšířený Kalmanův filtr definujeme jako algoritmus
výpočtů následujících estimátorů:
xt|t-1  N(t|t-1, t|t-1)(13)
xt|t  N(t|t, t|t)(14)
xt|N  N(t|N , t|N )(15)
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.21/3
Rozšířený Kalmanův filtr III
t|t = t|t-1 + Kt(yt - g(t|t-1, ut))(16)
t|t = t|t-1 - KtCtt|t-1(17)
t+1|t = f(t|t, ut)(18)
t+1|t = Att|tAt
T
+ v(19)
Kt = t|t-1Ct
T
(Ctt|t-1Ct
T
+ w)-1
(20)
t|N = t|t + Ft(t+1|N - t+1|t)(21)
t|N = t|t - Ft(t+1|t - t+1|N )Ft
T
(22)
Ft = t|tAt
T
-1
t+1|t
(23)
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.22/3
Rozšířený Kalmanův filtr IV
kde matice At, Ct jsou následující Jakobiány:
At =
f
x
(t|t) Ct =
g
x
(t|t-1)
Vztahy 19-20 představují filtrační (datový) krok.
Vztahy 21-22 představují predikční krok.
Vztahy 24-25 představují zpětný běh filtru
(smoothing).
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.23/3
Princip metody Monte Carlo I
Metoda Monte Carlo je založena na tom, že
informaci o rozložení náhodné veličiny nese náhodný
výběr z tohoto rozložení.
Čím větší je rozsah výběru, tím je informace
přesnější.
Pro náhodný vektor x  Lnx
máme náhodný výběr
(x1, . . . , xn) rozsahu n, který obsahuje vzorky
xi  Rnx
.
Empirická hustota pravděpodobnosti získaná
z náhodného výběru je aproximací skutečné hustoty
pravděpodobnosti náhodného vektoru x.
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.24/3
Princip metody Monte Carlo II
Empirickou hustotu pravděpodobnosti lze psát jako:
pn(x) =
1
n
n
i=1
(x - xi)
kde (x) : Rnx
 R je tzv. Diracova funkce.
Diracova funkce
(x) = lim
h0
h(x)
h(x) =
1
h
pro 0 < x < h, h(x) = 0 jinak
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.25/3
Princip metody Monte Carlo III
Dále přiřadíme jednotlivým vzorkům xi váhy wi  0
tak, aby suma všech vah byla rovna jedné.
Empirická hustota pravděpodobnosti je pak tvaru:
pn(x) =
n
i=1
wi (x - xi)
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.26/3
Vážený Bootstrap algoritmus I
Mějme dynamický systém tvaru:
xt+1 = f(xt, ut, vt)
yt = g(xt, ut, wt)
s počátečním stavem x0 a s empirickou hustotou
pn(x) =
n
i=1
wi(0| - 1)(x - xi(0| - 1))
Mějme data Dt. Pak odhady xt|t-1, xt|t, xt|N s empirickými
hustotami
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.27/3
Vážený Bootstrap algoritmus II
pn(xt|t-1) =
n
i=1
wi(t|t - 1)(x - xi(t|t - 1))
pn(xt|t) =
n
i=1
wi(t|t)(x - xi(t|t))
pn(xt|N ) =
n
i=1
wi(t|N)(x - xi(t|N))
jsou spočteny váženým bootstrap algoritmem právě když
platí:
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.28/3
Vážený Bootstrap algoritmus III
Filtrační (datový) krok: hustoty a váhy se
přepočítávají; vzorky zůstávají stejné.
xi(t|t) = xi(t|t - 1)
wi(t|t) = p(y(t)|xi(t|t - 1), u(t)) wi(t|t - 1)
wi(t|t) =
wi(t|t)
n
j=1 wj(t|t)
Predikční krok: hustoty a vzorky se přepočítávají
(vzhledem k jejich významnosti); váhy zůstávají
stejné.
xi(t + 1|t) = f(xi(t|t), u(t)) + vi(t)
wi(t + 1|t) = wi(t|t)
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.29/3
Vážený Bootstrap algoritmus IV
Zpětný běh (smoothing): hustoty a váhy se
přepočítávají; vzorky zůstávají stejné.
xi(t|N) = xi(t|t)
wi(t|N) = wi(t|t)
n
j=1
wj(t + 1|N) p(xj(t + 1|N)|xi(t|N), u(t
wi(t|N) =
wi(t|t)
n
j=1 wj(t|t)
4. Optimální odhad stavů systémualgoritmy + příklady ­ p.30/3