Skip to content
GitLab
Commit e3693e15 authored by Kamil Dedecius's avatar Kamil Dedecius
Browse files

typo

parent 7dab1f81
Hide whitespace changes
Inline Side-by-side
cviceni/8_lab/Strukturni_model.ipynb
View file @ e3693e15
%% Cell type:code id: tags:
``` python
import numpy as np
from scipy.stats import norm
import matplotlib.pylab as plt
import sys
np.set_printoptions(precision=2)
```
%% Cell type:code id: tags:
``` python
class KF():
def __init__(self, A, B, H, R, Q):
self.A = A
self.B = B
self.H = H
self.R = R
self.Q = Q
self.P = np.eye(self.Q.shape[0]) * 10000
self.x = np.ones(self.Q.shape[0])
self.log_x = []
def predict(self, u=None):
xminus = self.A.dot(self.x)
Pminus = self.A.dot(self.P).dot(self.A.T) + self.Q
self.x = xminus
self.P = Pminus
def correct(self, yt):
HPHT = self.H.dot(self.P).dot(self.H.T)
zavorka_inv = np.linalg.inv(HPHT + self.R)
K = self.P.dot(self.H.T).dot(zavorka_inv)
innovation = yt - self.H.dot(self.x)
xplus = self.x + K.dot(innovation)
zavorka = np.eye(K.shape[0]) - K.dot(self.H)
KRKT = K.dot(self.R).dot(K.T)
Pplus = zavorka.dot(self.P).dot(zavorka.T) + KRKT
self.x = xplus
self.P = Pplus
def log(self):
self.log_x.append(self.x)
```
%% Cell type:markdown id: tags:
# Strukturní modely
Už známe klasické modely časových řad typu ARIMA. Existuje ale ještě jedna možnost, jak časové řady modelovat, a sice prostřednictvím tzv. strukturních modelů. Ty jsou modely stavovými, ovšem namísto koeficientů AR či MA částí modelujme přímo vývoj střední hodnoty, sezónních složek atd.
V tomto cvičení (které je současně testem) si zkusíme jeden z nejjednodušších a přitom ohromně mocných modelů: lokální lineární trend. Opět budeme mít pozorovanou veličinu $y_t$, která má tentokrát střední hodnotu - level - $l_t$ a je zatížena normálním šumem $\varepsilon_t$. Střední hodnotu budeme rovněž modelovat, a sice jako veličinu, která je ovlivněna svou přechozí hodnotou $l_{t-1}$ a šumem $\xi_t$, ale vstupuje do ní i drift $b_t$. Jedná se tedy o náhodnou procházku s driftem. A abychom dokázali pochytat i dynamiku na střední hodnotě, budeme i samotný drift modelovat jako náhodnou procházku. V řeči statového modelu to bude vypadat následovně (dobře rozmyslete!):
%% Cell type:markdown id: tags:
$$
\begin{alignat}{2}
y_t &= \mu_t + \varepsilon_t, \qquad &\varepsilon \sim \mathcal{N}(0, \sigma_{\varepsilon}^2), \\
\mu_t &= \mu_{t-1} + \nu_{t-1} + \xi_t, \qquad &\xi_t \sim \mathcal{N}(0, \sigma_{\xi}^2), \\
\nu_t &= \nu_{t-1} + \zeta_t, \qquad &\zeta_t \sim \mathcal{N}(0, \sigma_{\zeta}^2).
\end{alignat}
$$
%% Cell type:markdown id: tags:
Je důležité myslet na to, že pozorovaná veličina je $y_t$, zatímco nepozorovatelné $l_t$ a $b_t$ budeme uvažovat jako stavy.
Zkusíme si systém nasimulovat se známými hodnotami variancí (v kódu níže zadávám směrodatné odchylky kvůli scipy!):
%% Cell type:code id: tags:
``` python
ndat = 300
npred = 50
sigma_epsilon = 10
sigma_xi = .9
sigma_zeta = .3
mu = np.zeros(ndat)
nu = np.zeros(ndat)
l = np.zeros(ndat)
b = np.zeros(ndat)
y = np.zeros(ndat)
for t in range(1, ndat):
nu[t] = nu[t-1] + norm.rvs(scale=sigma_zeta)
mu[t] = mu[t-1] + nu[t-1] + norm.rvs(scale=sigma_xi)
y[t] = mu[t] + norm.rvs(scale=sigma_epsilon)
b[t] = b[t-1] + norm.rvs(scale=sigma_zeta)
l[t] = l[t-1] + b[t-1] + norm.rvs(scale=sigma_xi)
y[t] = l[t] + norm.rvs(scale=sigma_epsilon)
plt.plot(y, label=r'$y_t$')
plt.legend()
plt.show()
```
%% Cell type:markdown id: tags:
Zkuste si pustit simulaci víckrát, ať se seznámíte s flexibilitou modelu. Rovněž zkuste měnit variance jednotlivých šumů a pozorujte, jak se změny projevují (chce to pokaždé více běhů).
%% Cell type:markdown id: tags:
## Úkol 1:
Využijte Kalmanův filtr pro modelování výše uvedeného systému. Tj. sestavte matice (uveďte je zde), pusťte odhad a udělejte predikce i na např. 50 kroků za horizont dat.
%% Cell type:code id: tags:
``` python
A = ***
H = ***
Q = ***
R = *** # pozor, musí být 2D
```
%% Cell type:code id: tags:
``` python
kf = KF(A=A, B=None, H=H, R=R, Q=Q)
for t in range(ndat+npred):
kf.predict()
if t < ndat:
kf.correct(np.atleast_1d(y[t]))
kf.log()
log_x = np.array(kf.log_x).T
```
%% Cell type:code id: tags:
``` python
plt.plot(y, label=r'$y_t$')
plt.plot(log_x[0].T, label=r'$\hat{\mu}_t$')
plt.legend()
```
%% Cell type:markdown id: tags:
## Úkol 2:
V souboru '20180224_Gp_part_5m.txt' naleznete měření z přístroje zjišťujícího tok elektronů a protonů v slunečním větru. Tento instrument je umístěn v družici [GOES-15](https://en.wikipedia.org/wiki/GOES_15). Vyberte sloupec s měřeními "Electrons at >2.0 Mev" a zkuste jej namodelovat pomocí modelu lokálního lineárního trendu. To znamená, že budete muset vyladit i (neznámé) hodnoty variancí šumu. Nedá to moc práce, ale chce to použít mozek :)
%% Cell type:code id: tags:
``` python
fn = './20180224_Gp_part_5m.txt'
data = np.genfromtxt(fn, skip_header=26)
E2 = data[:,-2]
plt.plot(E2)
```
%% Cell type:code id: tags:
``` python
sigma_epsilon = ***
sigma_xi = ***
sigma_zeta = ***
A = np.array([[1, 1],
[0, 1]])
H = np.array([[1, 0]])
R = np.atleast_2d(sigma_epsilon**2)
Q = np.diag([sigma_xi**2, sigma_zeta**2])
```
%% Cell type:code id: tags:
``` python
kf_E2 = KF(A=A, B=None, H=H, R=R, Q=Q)
kf_E2.x = [E2[0], 0]
kf_E2.P = np.eye(2) * 1000
print(kf_E2.x)
for t, E2t in enumerate(E2):
kf_E2.predict()
kf_E2.correct(np.atleast_1d(E2t))
kf_E2.log()
E2_log_x = np.array(kf_E2.log_x).T
print(E2_log_x[:,:5])
```
%% Cell type:code id: tags:
``` python
plt.plot(E2[:200], label=r'E>2')
plt.plot(E2_log_x[0,:200].T, label=r'$\hat{\mu}_t$')
plt.legend()
```
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment