In [2]:

    

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
from pandas import DataFrame
from numpy.random import multivariate_normal
from numpy.linalg  import inv
from pykalman import KalmanFilter
%matplotlib inline


    

In [3]:

    

def Relu(z):
    return np.max((z, np.zeros(z.shape)), axis=0)

def dRelu(z):
    return (z > 0) * 1

p = 3
q = 10
T = 500

data = loadmat('Tut11_file1.mat')
locals().update(data)
data.keys()


    

    
Out[3]:





dict_keys(['__header__', '__version__', '__globals__', 'A', 'B', 'Gamma', 'Sigma', 'W', 'mu0'])

In [4]:

    

z = np.zeros((p, T + 1))
x = np.zeros((q, T + 1))
np.random.seed(10)
z[:,[0]] = mu0 #+ multivariate_normal(np.zeros(p), Sigma).reshape(mu0.shape)

def transition_f(z):
    return A.dot(z) + W.dot(np.maximum(z, 0.))

def obs_g(x):
    return B.dot(x)

for t in range(T):
    x[:, t + 1] = obs_g(z[:, t]) #+ multivariate_normal(np.zeros(q), Gamma)
    z[:, t + 1] = transition_f(z[:, t]) #+ multivariate_normal(np.zeros(p), Sigma)
x[:, 0] = np.ma.masked


    

In [5]:

    

# Filter

mu = np.zeros(z.shape)
K = np.zeros((p, q, T))
V = np.zeros((p, p, T))
L = np.zeros((p, p, T))

L[..., 0] = Sigma
# mu[..., [0]] = mu0
K[...,0] = L[..., 0].dot(B.T.dot(inv(B.dot(L[..., 0].dot(B.T)) + Gamma)))
mu[..., [0]] = A.dot(mu0) + K[..., 0].dot(x[:, [0]] - B.dot(A.dot(mu0)))
V[..., 0] = (np.eye(p) - K[..., 0].dot(B)).dot(Sigma)
#L[..., 0] = A.dot(V[..., 0].dot(A.T)) + Sigma

for t in range(1, T):
    L[..., t] = A.dot(V[..., t - 1].dot(A.T)) + Sigma

    K[...,t] = L[..., t].dot(B.T.dot(inv(B.dot(L[..., t].dot(B.T)) + Gamma)))
    mu[..., [t]] = A.dot(mu[..., [t-1]]) + K[..., t].dot(x[:, [t]] - B.dot(A.dot(mu[..., [t-1]])))
    V[..., t] = (np.eye(p) - K[..., t].dot(B)).dot(L[..., t])
    #L[..., t] = A.dot(V[..., t].dot(A.T)) + Sigma


    

In [6]:

    

mu[:, 0:2]


    

    
Out[6]:





array([[-0.49309647, -0.42085202],
       [ 0.2719109 , -0.26398786],
       [ 0.61657095, -0.41404933]])

In [7]:

    

print ('Non smoothed result:', np.sum((mu[:, :-1] - z[:, 1:]).T ** 2))
print('Smoothed result:', np.sum((mu_tilde - z).T ** 2))
plt.plot(mu_tilde.T)


    

    




Non smoothed result: 2.06663480058






    




---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-63d22bea7055> in <module>()
      1 print ('Non smoothed result:', np.sum((mu[:, :-1] - z[:, 1:]).T ** 2))
----> 2 print('Smoothed result:', np.sum((mu_tilde - z).T ** 2))
      3 plt.plot(mu_tilde.T)

NameError: name 'mu_tilde' is not defined

In [134]:

    

mu[:, :-1].shape


    

    
Out[134]:





(3, 500)

In [146]:

    

# Filter

mu = np.zeros(z.shape)
K = np.zeros((p, q, T))
V = np.zeros((p, p, T))
L = np.zeros((p, p, T))

K[...,0] = Sigma.dot(B.T.dot(inv(B.dot(Sigma.dot(B.T)) + Gamma)))
mu[..., [0]] = A.dot(mu0) + W.dot(Relu(mu0)) + \
               K[..., 0].dot(x[:, [0]]  - B.dot(A.dot(mu0) + W.dot(Relu(mu0))))
V[..., 0] = (np.eye(p) - K[..., 0].dot(B)).dot(Sigma)
L[..., 0] = (A + W.dot(dRelu(mu0))).dot(V[..., 0].dot((A + W.dot(dRelu(mu0))).T)) + Sigma

for t in range(1, T):
    K[...,t] = L[..., t - 1].dot(B.T.dot(inv(B.dot(L[..., t - 1].dot(B.T)) + Gamma)))
    mu[..., [t]] = A.dot(mu[..., [t-1]]) + W.dot(Relu(mu[..., [t-1]])) + \
                   K[..., t].dot(x[:, [t]] - B.dot(A.dot(mu[..., [t-1]]) + W.dot(Relu(mu[..., [t-1]]))))
    V[..., t] = (np.eye(p) - K[..., t].dot(B)).dot(L[..., t-1])
    L[..., t] = (A + W.dot(dRelu(mu[..., [t-1]]))).dot(V[..., t].dot(A + W.dot(dRelu(mu[..., [t-1]])))) + Sigma