text:
```python
%matplotlib inline
```

=====================================================================
Spectro-temporal receptive field (STRF) estimation on continuous data
=====================================================================

This demonstrates how an encoding model can be fit with multiple continuous
inputs. In this case, we simulate the model behind a spectro-temporal receptive
field (or STRF). First, we create a linear filter that maps patterns in
spectro-temporal space onto an output, representing neural activity. We fit
a receptive field model that attempts to recover the original linear filter
that was used to create this data.

References
----------
Estimation of spectro-temporal and spatio-temporal receptive fields using
modeling with continuous inputs is described in:

.. [1] Theunissen, F. E. et al. Estimating spatio-temporal receptive
       fields of auditory and visual neurons from their responses to
       natural stimuli. Network 12, 289-316 (2001).

.. [2] Willmore, B. & Smyth, D. Methods for first-order kernel
       estimation: simple-cell receptive fields from responses to
       natural scenes. Network 14, 553-77 (2003).

.. [3] Crosse, M. J., Di Liberto, G. M., Bednar, A. & Lalor, E. C. (2016).
       The Multivariate Temporal Response Function (mTRF) Toolbox:
       A MATLAB Toolbox for Relating Neural Signals to Continuous Stimuli.
       Frontiers in Human Neuroscience 10, 604.
       doi:10.3389/fnhum.2016.00604

.. [4] Holdgraf, C. R. et al. Rapid tuning shifts in human auditory cortex
       enhance speech intelligibility. Nature Communications, 7, 13654 (2016).
       doi:10.1038/ncomms13654

```python
# Authors: Chris Holdgraf <choldgraf@gmail.com>
#          Eric Larson <larson.eric.d@gmail.com>
#
# License: BSD (3-clause)

# sphinx_gallery_thumbnail_number = 7

import numpy as np
import matplotlib.pyplot as plt

import mne
from mne.decoding import ReceptiveField, TimeDelayingRidge

from scipy.stats import multivariate_normal
from scipy.io import loadmat
from sklearn.preprocessing import scale
rng = np.random.RandomState(1337)  # To make this example reproducible
```

Load audio data
---------------

We'll read in the audio data from [3]_ in order to simulate a response.

In addition, we'll downsample the data along the time dimension in order to
speed up computation. Note that depending on the input values, this may
not be desired. For example if your input stimulus varies more quickly than
1/2 the sampling rate to which we are downsampling.

```python
# Read in audio that's been recorded in epochs.
path_audio = mne.datasets.mtrf.data_path()
data = loadmat(path_audio + '/speech_data.mat')
audio = data['spectrogram'].T
sfreq = float(data['Fs'][0, 0])
n_decim = 2
audio = mne.filter.resample(audio, down=n_decim, npad='auto')
sfreq /= n_decim
```

Create a receptive field
------------------------

We'll simulate a linear receptive field for a theoretical neural signal. This
defines how the signal will respond to power in this receptive field space.

```python
n_freqs = 20
tmin, tmax = -0.1, 0.4

# To simulate the data we'll create explicit delays here
delays_samp = np.arange(np.round(tmin * sfreq),
                        np.round(tmax * sfreq) + 1).astype(int)
delays_sec = delays_samp / sfreq
freqs = np.linspace(50, 5000, n_freqs)
grid = np.array(np.meshgrid(delays_sec, freqs))

# We need data to be shaped as n_epochs, n_features, n_times, so swap axes here
grid = grid.swapaxes(0, -1).swapaxes(0, 1)

# Simulate a temporal receptive field with a Gabor filter
means_high = [.1, 500]
means_low = [.2, 2500]
cov = [[.001, 0], [0, 500000]]
gauss_high = multivariate_normal.pdf(grid, means_high, cov)
gauss_low = -1 * multivariate_normal.pdf(grid, means_low, cov)
weights = gauss_high + gauss_low  # Combine to create the "true" STRF
kwargs = dict(vmax=np.abs(weights).max(), vmin=-np.abs(weights).max(),
              cmap='RdBu_r', shading='gouraud')

fig, ax = plt.subplots()
ax.pcolormesh(delays_sec, freqs, weights, **kwargs)
ax.set(title='Simulated STRF', xlabel='Time Lags (s)', ylabel='Frequency (Hz)')
plt.setp(ax.get_xticklabels(), rotation=45)
plt.autoscale(tight=True)
mne.viz.tight_layout()
```

Simulate a neural response
--------------------------

Using this receptive field, we'll create an artificial neural response to
a stimulus.

To do this, we'll create a time-delayed version of the receptive field, and
then calculate the dot product between this and the stimulus. Note that this
is effectively doing a convolution between the stimulus and the receptive
field. See `here <https://en.wikipedia.org/wiki/Convolution>`_ for more
information.

```python
# Reshape audio to split into epochs, then make epochs the first dimension.
n_epochs, n_seconds = 16, 5
audio = audio[:, :int(n_seconds * sfreq * n_epochs)]
X = audio.reshape([n_freqs, n_epochs, -1]).swapaxes(0, 1)
n_times = X.shape[-1]

# Delay the spectrogram according to delays so it can be combined w/ the STRF
# Lags will now be in axis 1, then we reshape to vectorize
delays = np.arange(np.round(tmin * sfreq),
                   np.round(tmax * sfreq) + 1).astype(int)

# Iterate through indices and append
X_del = np.zeros((len(delays),) + X.shape)
for ii, ix_delay in enumerate(delays):
    # These arrays will take/put particular indices in the data
    take = [slice(None)] * X.ndim
    put = [slice(None)] * X.ndim
    if ix_delay > 0:
        take[-1] = slice(None, -ix_delay)
        put[-1] = slice(ix_delay, None)
    elif ix_delay < 0:
        take[-1] = slice(-ix_delay, None)
        put[-1] = slice(None, ix_delay)
    X_del[ii][tuple(put)] = X[tuple(take)]

# Now set the delayed axis to the 2nd dimension
X_del = np.rollaxis(X_del, 0, 3)
X_del = X_del.reshape([n_epochs, -1, n_times])
n_features = X_del.shape[1]
weights_sim = weights.ravel()

# Simulate a neural response to the sound, given this STRF
y = np.zeros((n_epochs, n_times))
for ii, iep in enumerate(X_del):
    # Simulate this epoch and add random noise
    noise_amp = .002
    y[ii] = np.dot(weights_sim, iep) + noise_amp * rng.randn(n_times)

# Plot the first 2 trials of audio and the simulated electrode activity
X_plt = scale(np.hstack(X[:2]).T).T
y_plt = scale(np.hstack(y[:2]))
time = np.arange(X_plt.shape[-1]) / sfreq
_, (ax1, ax2) = plt.subplots(2, 1, figsize=(6, 6), sharex=True)
ax1.pcolormesh(time, freqs, X_plt, vmin=0, vmax=4, cmap='Reds')
ax1.set_title('Input auditory features')
ax1.set(ylim=[freqs.min(), freqs.max()], ylabel='Frequency (Hz)')
ax2.plot(time, y_plt)
ax2.set(xlim=[time.min(), time.max()], title='Simulated response',
        xlabel='Time (s)', ylabel='Activity (a.u.)')
mne.viz.tight_layout()
```

Fit a model to recover this receptive field
-------------------------------------------

Finally, we'll use the :class:`mne.decoding.ReceptiveField` class to recover
the linear receptive field of this signal. Note that properties of the
receptive field (e.g. smoothness) will depend on the autocorrelation in the
inputs and outputs.

```python
# Create training and testing data
train, test = np.arange(n_epochs - 1), n_epochs - 1
X_train, X_test, y_train, y_test = X[train], X[test], y[train], y[test]
X_train, X_test, y_train, y_test = [np.rollaxis(ii, -1, 0) for ii in
                                    (X_train, X_test, y_train, y_test)]
# Model the simulated data as a function of the spectrogram input
alphas = np.logspace(-3, 3, 7)
scores = np.zeros_like(alphas)
models = []
for ii, alpha in enumerate(alphas):
    rf = ReceptiveField(tmin, tmax, sfreq, freqs, estimator=alpha)
    rf.fit(X_train, y_train)

    # Now make predictions about the model output, given input stimuli.
    scores[ii] = rf.score(X_test, y_test)
    models.append(rf)

times = rf.delays_ / float(rf.sfreq)

# Choose the model that performed best on the held out data
ix_best_alpha = np.argmax(scores)
best_mod = models[ix_best_alpha]
coefs = best_mod.coef_[0]
best_pred = best_mod.predict(X_test)[:, 0]

# Plot the original STRF, and the one that we recovered with modeling.
_, (ax1, ax2) = plt.subplots(1, 2, figsize=(6, 3), sharey=True, sharex=True)
ax1.pcolormesh(delays_sec, freqs, weights, **kwargs)
ax2.pcolormesh(times, rf.feature_names, coefs, **kwargs)
ax1.set_title('Original STRF')
ax2.set_title('Best Reconstructed STRF')
plt.setp([iax.get_xticklabels() for iax in [ax1, ax2]], rotation=45)
plt.autoscale(tight=True)
mne.viz.tight_layout()

# Plot the actual response and the predicted response on a held out stimulus
time_pred = np.arange(best_pred.shape[0]) / sfreq
fig, ax = plt.subplots()
ax.plot(time_pred, y_test, color='k', alpha=.2, lw=4)
ax.plot(time_pred, best_pred, color='r', lw=1)
ax.set(title='Original and predicted activity', xlabel='Time (s)')
ax.legend(['Original', 'Predicted'])
plt.autoscale(tight=True)
mne.viz.tight_layout()
```

Visualize the effects of regularization
---------------------------------------

Above we fit a :class:`mne.decoding.ReceptiveField` model for one of many
values for the ridge regularization parameter. Here we will plot the model
score as well as the model coefficients for each value, in order to
visualize how coefficients change with different levels of regularization.
These issues as well as the STRF pipeline are described in detail
in [1]_, [2]_, and [4]_.

```python
# Plot model score for each ridge parameter
fig = plt.figure(figsize=(10, 4))
ax = plt.subplot2grid([2, len(alphas)], [1, 0], 1, len(alphas))
ax.plot(np.arange(len(alphas)), scores, marker='o', color='r')
ax.annotate('Best parameter', (ix_best_alpha, scores[ix_best_alpha]),
            (ix_best_alpha, scores[ix_best_alpha] - .1),
            arrowprops={'arrowstyle': '->'})
plt.xticks(np.arange(len(alphas)), ["%.0e" % ii for ii in alphas])
ax.set(xlabel="Ridge regularization value", ylabel="Score ($R^2$)",
       xlim=[-.4, len(alphas) - .6])
mne.viz.tight_layout()

# Plot the STRF of each ridge parameter
for ii, (rf, i_alpha) in enumerate(zip(models, alphas)):
    ax = plt.subplot2grid([2, len(alphas)], [0, ii], 1, 1)
    ax.pcolormesh(times, rf.feature_names, rf.coef_[0], **kwargs)
    plt.xticks([], [])
    plt.yticks([], [])
    plt.autoscale(tight=True)
fig.suptitle('Model coefficients / scores for many ridge parameters', y=1)
mne.viz.tight_layout()
```

Using different regularization types
------------------------------------
In addition to the standard ridge regularization, the
:class:`mne.decoding.TimeDelayingRidge` class also exposes
`Laplacian <https://en.wikipedia.org/wiki/Laplacian_matrix>`_ regularization
term as:

\begin{align}\left[\begin{matrix}
        1 & -1 &   &   & & \\
       -1 &  2 & -1 &   & & \\
          & -1 & 2 & -1 & & \\
          & & \ddots & \ddots & \ddots & \\
          & & & -1 & 2 & -1 \\
          & & &    & -1 & 1\end{matrix}\right]\end{align}

This imposes a smoothness constraint of nearby time samples and/or features.
Quoting [3]_:

   Tikhonov [identity] regularization (Equation 5) reduces overfitting by
   smoothing the TRF estimate in a way that is insensitive to
   the amplitude of the signal of interest. However, the Laplacian
   approach (Equation 6) reduces off-sample error whilst preserving
   signal amplitude (Lalor et al., 2006). As a result, this approach
   usually leads to an improved estimate of the system’s response (as
   indexed by MSE) compared to Tikhonov regularization.

```python
scores_lap = np.zeros_like(alphas)
models_lap = []
for ii, alpha in enumerate(alphas):
    estimator = TimeDelayingRidge(tmin, tmax, sfreq, reg_type='laplacian',
                                  alpha=alpha)
    rf = ReceptiveField(tmin, tmax, sfreq, freqs, estimator=estimator)
    rf.fit(X_train, y_train)

    # Now make predictions about the model output, given input stimuli.
    scores_lap[ii] = rf.score(X_test, y_test)
    models_lap.append(rf)

ix_best_alpha_lap = np.argmax(scores_lap)
```

Compare model performance
-------------------------
Below we visualize the model performance of each regularization method
(ridge vs. Laplacian) for different levels of alpha. As you can see, the
Laplacian method performs better in general, because it imposes a smoothness
constraint along the time and feature dimensions of the coefficients.
This matches the "true" receptive field structure and results in a better
model fit.

```python
fig = plt.figure(figsize=(10, 6))
ax = plt.subplot2grid([3, len(alphas)], [2, 0], 1, len(alphas))
ax.plot(np.arange(len(alphas)), scores_lap, marker='o', color='r')
ax.plot(np.arange(len(alphas)), scores, marker='o', color='0.5', ls=':')
ax.annotate('Best Laplacian', (ix_best_alpha_lap,
                               scores_lap[ix_best_alpha_lap]),
            (ix_best_alpha_lap, scores_lap[ix_best_alpha_lap] - .1),
            arrowprops={'arrowstyle': '->'})
ax.annotate('Best Ridge', (ix_best_alpha, scores[ix_best_alpha]),
            (ix_best_alpha, scores[ix_best_alpha] - .1),
            arrowprops={'arrowstyle': '->'})
plt.xticks(np.arange(len(alphas)), ["%.0e" % ii for ii in alphas])
ax.set(xlabel="Laplacian regularization value", ylabel="Score ($R^2$)",
       xlim=[-.4, len(alphas) - .6])
mne.viz.tight_layout()

# Plot the STRF of each ridge parameter
xlim = times[[0, -1]]
for ii, (rf_lap, rf, i_alpha) in enumerate(zip(models_lap, models, alphas)):
    ax = plt.subplot2grid([3, len(alphas)], [0, ii], 1, 1)
    ax.pcolormesh(times, rf_lap.feature_names, rf_lap.coef_[0], **kwargs)
    ax.set(xticks=[], yticks=[], xlim=xlim)
    if ii == 0:
        ax.set(ylabel='Laplacian')
    ax = plt.subplot2grid([3, len(alphas)], [1, ii], 1, 1)
    ax.pcolormesh(times, rf.feature_names, rf.coef_[0], **kwargs)
    ax.set(xticks=[], yticks=[], xlim=xlim)
    if ii == 0:
        ax.set(ylabel='Ridge')
fig.suptitle('Model coefficients / scores for laplacian regularization', y=1)
mne.viz.tight_layout()
```

Plot the original STRF, and the one that we recovered with modeling.

```python
rf = models[ix_best_alpha]
rf_lap = models_lap[ix_best_alpha_lap]
_, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(9, 3),
                                  sharey=True, sharex=True)
ax1.pcolormesh(delays_sec, freqs, weights, **kwargs)
ax2.pcolormesh(times, rf.feature_names, rf.coef_[0], **kwargs)
ax3.pcolormesh(times, rf_lap.feature_names, rf_lap.coef_[0], **kwargs)
ax1.set_title('Original STRF')
ax2.set_title('Best Ridge STRF')
ax3.set_title('Best Laplacian STRF')
plt.setp([iax.get_xticklabels() for iax in [ax1, ax2, ax3]], rotation=45)
plt.autoscale(tight=True)
mne.viz.tight_layout()
```


text:
# TSSL Lab 2 - Structural model, Kalman filtering and EM
We will continue to work with the Global Mean Sea Level (GMSL) data that we got acquainted with in lab 1. The data is taken from https://climate.nasa.gov/vital-signs/sea-level/ and is available on LISAM in the file `sealevel.csv`.

In this lab we will analyse this data using a structural time series model. We will first set up a model and implement a Kalman filter to infer the latet states of the model, as well doing long-term prediction. We will then implement a disturbance smoother and an expectation maximization algorithm to tune the parameters of the model.

We load a few packages that are useful for solving this lab assignment.

```python
import pandas  # Loading data / handling data frames
import numpy as np
import scipy.linalg
import math
import matplotlib.pyplot as plt
plt.rcParams["figure.figsize"] = (12,8)  # Increase default size of plots
```

## 2.1 Setting up a structural state space model

We start by loading and plotting data to remind ourselves what it looks like.

```python
data=pandas.read_csv('sealevel.csv',header=0)
```

```python
y = data['GMSL'].values
u = data['Year'].values
ndata = len(y)
plt.plot(u,y)
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
plt.show()
```

In this lab we will use a structural time series model to analys this data set. Specifically, we assume that the data $\{y_t\}_{t\geq 1}$ is generated by

\begin{align}
    y_t = \mu_t + \gamma_t + \varepsilon_t
\end{align}

where $\mu_t$ is a trend component, $\gamma_t$ is a seasonal component, and $\varepsilon_t$ is an observation noise. The model is expressed using a state space representation,

\begin{align}
    \alpha_{t+1} &= T \alpha_t + R\eta_t, & \eta_t&\sim N(0,Q), \\
    y_t &= Z \alpha_t + \varepsilon_t, & \varepsilon_t&\sim N(0,\sigma_\varepsilon^2).
\end{align}

**Q0:** Let $d = \dim(\alpha_t)$ denote the _state dimension_ and $d_\eta = \dim(\eta_t)$ denote the dimension of the state noise. Then, what are the dimenisons of the matrices $T$, $R$, and $Z$ of the state space model?

**A:** The matrices of the state space model can be expressed in terms of the dimensions of the state and state noise as follows:

$T$ has dimensions $d$ x $d$
$R$ has dimensions $d$ x $d_\eta$
$Z$ has dimensions $1$ x $d$

**Q1:** Create the state space matrices $T_{[\mu]}$, $R_{[\mu]}$, and $Z_{[\mu]}$ corresponding to the trend component $\mu_t$. We assume a local linear trend (that is, of order $k=2$).

_Hint:_ Use **2-dimensional** `numpy.ndarray`s of the correct sizes to represent all the matrices.

```python
# State space matrices corresponding to the trend component
k = 2
T_mu = np.zeros(shape=(k, k))
R_mu = np.zeros(shape=(k, 1))
Z_mu = np.zeros(shape=(1, k))

# Function to compute factorial
def nCr(n,r):
    f = math.factorial
    return f(n) / f(r) / f(n-r)

# Fill in the state space matrices, T, R and Z
T_mu[0, ] = [(-1)**(i+2) * nCr(k, i+1) for i in range(k)]
for i in range(1, k):
    T_mu[i, i-1] = 1

R_mu[0, 0] = 1
Z_mu[0, 0] = 1
```

```python
print("Matrix T for the trend component : \n", T_mu)
print("Matrix R for the trend component : \n", R_mu)
print("Matrix Z for the trend component : \n", Z_mu)
```

    Matrix T for the trend component :
     [[ 2. -1.]
     [ 1.  0.]]
    Matrix R for the trend component :
     [[1.]
     [0.]]
    Matrix Z for the trend component :
     [[1. 0.]]

**Q2:** There is a yearly seasonal pattern present in the data. What should we set the periodicity $s$ of the seasonal component to, to capture this pattern?

_Hint:_ Count the average number of observations per (whole) year and round to the closest integer.

```python
# Calculate average counts per year in the observations
years = data['Year'].astype(str).str[0:4]
yearCounts = years.value_counts()

print(yearCounts)
print("Average number of observations per year : ", yearCounts.mean())
```

    2017    37
    2019    37
    2009    37
    2007    37
    2013    37
    2012    37
    1994    37
    2000    37
    1993    37
    2003    37
    2005    37
    1997    37
    1999    37
    2001    37
    2015    37
    2011    37
    1995    37
    2006    37
    2016    37
    2018    37
    1998    37
    2010    37
    2004    37
    1996    36
    2002    36
    2014    36
    2008    36
    2020     2
    Name: Year, dtype: int64
    Average number of observations per year :  35.607142857142854

```python
# Calculate average counts per year in the observations excluding the two data points at the end limit
yearCounts = years[:-2].value_counts()
print("Average number of observations per year (excluding the two observations made in 2020): ", yearCounts.mean())
```

    Average number of observations per year (excluding the two observations made in 2020):  36.851851851851855

The periodicity s should be set to 37 in order to capture the yearly seasonal pattern present in the data.

**Q3:** What is the _state dimension_ of a seasonal component with periodicity $s$? That is, how many states are needed in the corresponding state space representation?

**A:** The state dimension of the seasonal component is (s-1)*(s-1) = 36 x 36 i.e. 36 states are needed in the state space representation of a seasonal component model with periodicity s=37.

**Q4:** Create the state space matrices $T_{[\gamma]}$, $R_{[\gamma]}$, and $Z_{[\gamma]}$ corresponding to the seasonal component $\gamma_t$.

_Hint:_ Use **2-dimensional** `numpy.ndarray`s of the correct sizes to represent all the matrices.

```python
# State space matrices corresponding to the seasonal component
s = 37
T_gamma = np.zeros(shape=(s-1, s-1))
R_gamma = np.zeros(shape=(s-1, 1))
Z_gamma = np.zeros(shape=(1, s-1))

# Fill in the state space matrices, T, R and Z
T_gamma[0, ] = -1
for i in range(1, s-1):
    T_gamma[i, i-1] = 1

R_gamma[0, 0] = 1
Z_gamma[0, 0] = 1
```

```python
print("Matrix T for the seasonal component : \n", T_gamma)
print("Matrix R for the seasonal component : \n", R_gamma)
print("Matrix Z for the seasonal component : \n", Z_gamma)
```

    Matrix T for the seasonal component :
     [[-1. -1. -1. ... -1. -1. -1.]
     [ 1.  0.  0. ...  0.  0.  0.]
     [ 0.  1.  0. ...  0.  0.  0.]
     ...
     [ 0.  0.  0. ...  0.  0.  0.]
     [ 0.  0.  0. ...  1.  0.  0.]
     [ 0.  0.  0. ...  0.  1.  0.]]
    Matrix R for the seasonal component :
     [[1.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]
     [0.]]
    Matrix Z for the seasonal component :
     [[1. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
      0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

**Q5:** Using the matrices that you have constructed above, create the state space matrices for the complete structural time series model. Print out the shapes of the resulting system matrices and check that they correspond to what you expect (cf **Q0**).

_Hint:_ Use `scipy.linalg.block_diag` and `numpy.concatenate`.

```python
# State space matrices for the complete structural time series model
T = scipy.linalg.block_diag(T_mu, T_gamma)
R = scipy.linalg.block_diag(R_mu, R_gamma)
Z = np.concatenate((Z_mu, Z_gamma), axis=1)

# Print the dimensions
print("Dimensions of T : ", T.shape)
print("Dimensions of R : ", R.shape)
print("Dimensions of Z : ", Z.shape)
```

    Dimensions of T :  (38, 38)
    Dimensions of R :  (38, 2)
    Dimensions of Z :  (1, 38)

The dimensions of the state space matrices for the complete structural time series model checks out with our expectations in question 0. Here $d = 38$ and $d_\eta = 2$.

We also need to specify the variances of the process noise $\eta_t$ and measurement noise $\varepsilon_t$. Below, we will estimate (two of) these variances from data, but for now we set them arbitrarily to get an initial model to work with.

```python
# Some arbitrary noise values for now
sigma_trend = 0.01
sigma_seas = 1
sigma_eps = 1
Q = np.array([[sigma_trend**2, 0.], [0., sigma_seas**2]])  # Process noise covariance matrix
```

Finally, to complete the model we need to specify the distribution of the initial state. This encodes our _a priori_ belief about the actual values of the trend and seasonality, i.e., before observing any data.

**Q6:** Set up the mean vector of the initial state $a_1 = \mathbb{E}[\alpha_1]$ such that:
* The trend component starts at the first observation, $\mathbb{E}[\mu_1] = y_1$,
* The slope of the trend is _a priori_ zero in expectation, $\mathbb{E}[\mu_1 - \mu_0] = 0$,
* The initial mean of all states related to the seasonal component are zero.

Also, create an initial state covariance matrix $P_1 = \text{Cov}(\alpha_1)$ as an identity matrix of the correct dimension, multiplied with a large value (say, 100) to represent our uncertainty about the initial state.

```python
# Initial state vector
d = k+s-1
a1 = np.zeros(shape=(d, 1))
a1[0, 0] = y[0]

# Initial state covariance matrix
P1 = np.identity(d) * 100
```

We have now defined all the matrices etc. that make up the structural state space model. For convenience, we can create an object of the class `LGSS` available in the module `tssltools_lab2` as a container for these quantities.

```python
from tssltools_lab2 import LGSS  # Module available in LISAM
model = LGSS(T, R, Q, Z, sigma_eps**2, a1, P1)

help(model.get_params)
```

    Help on method get_params in module tssltools_lab2:

    get_params() method of tssltools_lab2.LGSS instance
        Return all model parameters.

        T, R, Q, Z, H, a1, P1 = model.get_params()

## 2.2 Kalman filtering for the structural model
Now we have the data and a model available. Next, we will turn our attention to the inference problem, which is a central task when analysing time series data using the state space framework.

State inference is the problem of estimating the unknown (latent) state variables given the data. For the time being we assume that the _model parameters_ are completely specified, according to above, and only consider how to estimate the states using the Kalman filter.

In the questions below we will treat the first $n=800$ time steps as training data and the remaining $m=197$ observations as validation data.

```python
n = 800
m = ndata-n
```

**Q7:** Complete the Kalman filter implementation below. The function should be able to handle missing observations, which are encoded as "not a number", i.e. `y[t] = np.nan` for certain time steps `t`.

_Hint:_ The Kalman filter involves a lot of matrix-matrix and matrix-vector multiplications. It turns out to be convient to store sequences of vectors (such as the predicted and filtered state estimates) as `(d,1,n)` arrays, instead of `(d,n)` or `(n,d)` arrays. In this way the matrix multiplications will result in 2d-arrays of the correct shapes without having to use a lot of explicit `reshape`. However, clearly, this is just a matter of coding style preferences!

```python
from tssltools_lab2 import kfs_res  # Module available in LISAM. kfs_res is a container class for storing the result.

def kalman_filter(y, model: LGSS):
    """Kalman filter for LGSS model with one-dimensional observation.

    :param y: (n,) array of observations. May contain nan, which encodes missing observations.
    :param model: LGSS object with the model specification.

    :return kfs_res: Container class with member variables,
        alpha_pred: (d,1,n) array of predicted state means.
        P_pred: (d,d,n) array of predicted state covariances.
        alpha_filt: (d,1,n) array of filtered state means.
        P_filt: (d,d,n) array of filtered state covariances.
        y_pred: (n,) array of means of p(y_t | y_{1:t-1})
        F_pred: (n,) array of variances of p(y_t | y_{1:t-1})
    """

    n = len(y)
    d = model.d  # State dimension
    alpha_pred = np.zeros((d, 1, n))
    P_pred = np.zeros((d, d, n))
    alpha_filt = np.zeros((d, 1, n))
    P_filt = np.zeros((d, d, n))
    y_pred = np.zeros(n)
    F_pred = np.zeros(n)

    T, R, Q, Z, H, a1, P1 = model.get_params()  # Get all model parameters (for brevity)

    for t in range(n):
        # Time update (predict)
        if t == 0:  # Initialize predictions at first time step
            alpha_pred[:, :, t] = a1
            P_pred[:, :, t]     = P1

        else:  # All consecutive time steps
            alpha_pred[:, :, t] = T.dot(alpha_filt[:, :, t-1])
            P_pred[:, :, t]     = T.dot(P_filt[:, :, t-1].dot(T.T)) + R.dot(Q.dot(R.T))

        # Compute prediction of current output
        y_pred[t] = Z.dot(alpha_pred[:, :, t])
        F_pred[t] = Z.dot(P_pred[:, :, t].dot(Z.T)) + H

        # Measurement update (correct)
        if np.isnan(y[t]):  # Handle missing data
            alpha_filt[:, :, t] = alpha_pred[:, :, t]
            P_filt[:, :, t]     = P_pred[:, :, t]

        else:
            K = P_pred[:, :, t].dot(Z.T) / F_pred[t] # Kalman gain
            alpha_filt[:, :, t] = alpha_pred[:, :, t] + K.dot(y[t]-y_pred[t])
            P_filt[:, :, t]     = (np.identity(d) - K.dot(Z)).dot(P_pred[:, :, t])

    # Container for storing all results
    kf = kfs_res(alpha_pred, P_pred, alpha_filt, P_filt, y_pred, F_pred)
    return kf
```

**Q8:** Use the Kalman filter to infer the states of the structural time series applied to the sealevel data. Run the filter on the training data (i.e., first $n=800$ time steps), followed by a long-range prediction of $y_t$ for the remaining time points.

Generate a plot which shows:
1. The data $y_{1:n+m}$,
2. The one-step predictions $\hat y_{t|t-1} \pm 1$ standard deviation for the training data, i.e., $t = 1,...,n$,
3. The long-range predictions $\hat y_{t|n} \pm 1$ standard deviation for the validation data, i.e., $t= n+1,...,n+m$,
4. A vertical line indicating the switch between training and validation data, using, e.g., `plt.axvline(x=u[n])`.

_Hint:_ It is enough to call the `kalman_filter` function once. Make use of the missing data functionality!

```python
# Fill NaNs for m time steps to predict using the kalman filter
y_new = np.concatenate((y[:n], np.repeat(np.nan, m)))
kf = kalman_filter(y=y_new, model=model)

# Plot the data and predictions
plt.plot(u, y, label='Data')
plt.plot(u, kf.y_pred, label='Predictions')
plt.plot(u, kf.y_pred + np.sqrt(kf.F_pred), label='Upper CI', ls='-.')
plt.plot(u, kf.y_pred - np.sqrt(kf.F_pred), label='Lower CI', ls='-.')
plt.axvline(x=u[n], color='grey', ls='--', label='Switch between training and validation data')
plt.title('Predictions')
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
```

**Q9:** Based on the output of the Kalman filter, compute the training data log-likelihood $\log p(y_{1:n})$.

```python
# Compute log-likelihood
def log_likelihood(y, y_hat, F):
    llik = -n/2*np.log(2*np.pi) - 1/2*np.sum(np.log(abs(F))+(y-y_hat).T/F*(y-y_hat))
    return(llik)

# Log-likelihood of training data based on the output of Kalman filter
llik = log_likelihood(y=y[:n], y_hat=kf.y_pred[:n], F=kf.F_pred[:n])
print("The log-likelihood of the training data : ", llik)
```

    The log-likelihood of the training data :  -2848.96385321561

## 2.3 Identifying the noise variances using the EM algorithm
So far we have used fixed model parameters when running the filter. In this section we will see how the model parameters can be learnt from data using the EM algorithm. Specifically, we will try to learn the variance of the state noise affecting the seasonal component as well as the variance of the observation noise,

\begin{align}
    \theta = (\sigma_\gamma^2, \sigma_\varepsilon^2).
\end{align}

For brevity, the variance of the trend component $\sigma_\mu^2$ is fixed to the value $\sigma_\mu^2 = 0.01^2$ as above. (See Appendix A below for an explanation.)

Recall that we consider $y_{1:n}$ as the training data, i.e., we will estimate $\theta$ using only the first $n=800$ observations.

**Q10:** Which optimization problem is it that the EM algorithm is designed to solve? Complete the line below!

**A:**

$\hat \theta = \arg\max_{\theta} \mathbb{E} \left[ log p_{\theta} (\alpha_{1:n}, y_{1:n}) | y_{1:n}, \tilde \theta \right] = const - \frac {1}{2} \sum \limits_{t=1}^{n} \left[ log |\sigma_{\epsilon}^2| + log |Q| + \left\{ \hat \epsilon_{t|n}^2 + Var[\epsilon_t | y_{1:n}] \right\} \sigma_{\epsilon}^{-2} + tr \left[ \left\{ \hat \eta_{t|n} \hat \eta_{t|n}^T + Var[\eta | y_{1:n}] \right\} Q^{-1} \right] \right] $

**Q11:** Write down the updating equations on closed form for the M-step in the EM algorithm.

_Hint: Look at Exercise Session 2_

**A:**

$\hat \sigma_{\epsilon}^2 = \frac{1}{n} \sum \limits_{t=1}^{n} \left\{ \hat \epsilon_{t|n}^2 + Var[\epsilon_t | y_{1:n}] \right\}$

$\hat \sigma_{\gamma}^2 = \frac{1}{n} \sum \limits_{t=1}^{n} \left\{ \hat \eta_{t|n} \hat \eta_{t|n}^T + Var[\eta | y_{1:n}] \right\} [d \eta, d \eta]$

To implement the EM algorithm we need to solve a _smoothing problem_. The Kalman filter that we implemented above is based only on a forward propagation of information. The _smoother_ complements the forward filter with a backward pass to compute refined state estimates. Specifically, the smoothed state estimates comprise the mean and covariances of

\begin{align}
    &p(\alpha_t \mid y_{1:n}), & t=&1,\dots,n
\end{align}

Furthermore, when implementing the EM algorithm it is convenient to work with the (closely related) smoothed estimates of the disturbances, i.e., the state and measurement noise,

\begin{align}
    &p(\eta_t \mid y_{1:n}), & t=&1,\dots,n-1 \\
    &p(\varepsilon_t \mid y_{1:n}), & t=&1,\dots,n
\end{align}

An implementation of a state and disturbance smoother is available in the `tssltools_lab2` module. You may use this when implementing the EM algorithm below.

```python
from tssltools_lab2 import kalman_smoother
help(kalman_smoother)
```

    Help on function kalman_smoother in module tssltools_lab2:

    kalman_smoother(y, model: tssltools_lab2.LGSS, kf: tssltools_lab2.kfs_res)
        Kalman (state and disturbance) smoother for LGSS model with one-dimensional observation.

        :param y: (n,) array of observations. May contain nan, which encodes missing observations.
        :param model: LGSS object with the model specification.
        :parma kf: kfs_res object with result from a Kalman filter foward pass.

        :return kfs_res: Container class. The original Kalman filter result is augmented with the following member variables,
            alpha_sm: (d,1,n) array of smoothed state means.
            V: (d,d,n) array of smoothed state covariances.
            eps_hat: (n,) array of smoothed means of observation disturbances.
            eps_var: (n,) array of smoothed variances of observation disturbances.
            eta_hat: (deta,1,n) array of smoothed means of state disturbances.
            eta_cov: (deta,deta,n) array of smoothed covariances of state disturbances.

```python
# Fill NaNs for m time steps to predict using the kalman filter and smoother
y_new = np.concatenate((y[:n], np.repeat(np.nan, m)))

# Run the kalman filter
kf = kalman_filter(y=y_new, model=model)

# Run the kalman smoother
ks = kalman_smoother(y=y_new, model=model, kf=kf)
```

**Q12:** Implement an EM algorithm by completing the code below. Run the algorithm for 100 iterations and plot the traces of the parameter estimates, i.e., the values $\theta_r$, for $r = 0,\dots,100$.

_Note:_ When running the Kalman filter as part of the EM loop you should only filter the _training data_ (i.e. excluding the prediction for validation data).

```python
num_iter = 100

def compute_theta(ks, Q):
    """Implements the M-step of the EM-algorithm, based on the Kalman smoother results (ks)"""

    eps_hat  = ks.eps_hat
    eps_var  = ks.eps_var
    eta_hat  = ks.eta_hat
    eta_cov  = ks.eta_cov

    n = len(eps_hat)
    sigma_eps_hat = 1/n * np.sum(eps_hat**2 + eps_var)

    deta = eta_hat.shape[0]
    sigma_seas_hat = 0
    for t in range(n):
        eta_var_hat    = eta_hat[:, :, t].dot(eta_hat[:, :, t].T) + eta_cov[:, :, t]
        sigma_seas_hat = sigma_seas_hat + eta_var_hat[deta-1, deta-1]
    sigma_seas_hat = sigma_seas_hat/n

    theta_new = [np.sqrt(sigma_eps_hat), np.sqrt(sigma_seas_hat)]
    return theta_new

sigma_seas_hat = 1
sigma_eps_hat  = 1

for r in range(1, num_iter):
    # E-step
    Q_hat = np.array([[sigma_trend**2, 0.], [0., sigma_seas_hat**2]])
    model = LGSS(T, R, Q_hat, Z, sigma_eps_hat**2, a1, P1)

    kf = kalman_filter(y=y[:n], model=model)
    ks = kalman_smoother(y=y[:n], model=model, kf=kf)

    # M-step
    theta_hat = compute_theta(ks, Q_hat)
    sigma_eps_hat  = theta_hat[0]
    sigma_seas_hat = theta_hat[1]

```

```python
print("Estimated sigma_eps_hat : ", sigma_eps_hat,
      "\nEstimated sigma_seas_hat : ", sigma_seas_hat)
```

    Estimated sigma_eps_hat :  2.7389155054413994
    Estimated sigma_seas_hat :  0.19126034276958498

## 2.4 Further analysing the data
We will now fix the model according to the final output from the EM algorithm and further analyse the data using this model.

**Q13:** Rerun the Kalman filter to compute a _long range prediction for the validation data points,_ analogously to **Q8** (you can copy-paste code from that question). That is, generate a plot which shows:
1. The data $y_{1:n+m}$,
2. The one-step predictions $\hat y_{t|t-1} \pm 1$ standard deviation for the training data, i.e., $t = 1,...,n$,
3. The long-range predictions $\hat y_{t|n} \pm 1$ standard deviation for the validation data, i.e., $t= n+1,...,n+m$,
4. A vertical line indicating the switch between training and validation data, using, e.g., `plt.axvline(x=u[n])`.

Furthermore, compute the training data log-likelihood $\log p(y_{1:n})$ using the estimated model (cf. **Q9**).

```python
# Fix the model according to the final output from the EM algorithm
Q_hat = np.array([[sigma_trend**2, 0.], [0., sigma_seas_hat**2]])
model = LGSS(T, R, Q_hat, Z, sigma_eps_hat**2, a1, P1)

# Run the kalman filter
y_new = np.concatenate((y[:n], np.repeat(np.nan, m)))
kf = kalman_filter(y=y_new, model=model)
```

```python
# Plot the data and predictions
plt.plot(u, y, label='Data')
plt.plot(u, kf.y_pred, label='Predictions')
plt.plot(u, kf.y_pred + np.sqrt(kf.F_pred), label='Upper CI', ls='-.')
plt.plot(u, kf.y_pred - np.sqrt(kf.F_pred), label='Lower CI', ls='-.')
plt.axvline(x=u[n], color='grey', ls='--', label='Switch between training and validation data')
plt.title('Predictions')
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
```

```python
# Log-likelihood of training data based on updated filter predictions
llik = log_likelihood(y=y[:n], y_hat=kf.y_pred[:n], F=kf.F_pred[:n])
print("The log-likelihood of the training data : ", llik)

# Log-likelihood of validation data based on updated filter predictions
llik = log_likelihood(y=y[n:], y_hat=kf.y_pred[n:], F=kf.F_pred[n:])
print("The log-likelihood of the validation data : ", llik)
```

    The log-likelihood of the training data :  -2112.0321419185084
    The log-likelihood of the validation data :  -1227.305036587401

Note that we can view the model for the data $y_t$ as being comprised of an underlying "signal", $s_t = \mu_t + \gamma_t$ plus observation noise $\varepsilon_t$

\begin{align}
    y_t = s_t + \varepsilon_t
\end{align}

We can obtain refined, _smoothed,_ estimates of this signal by conditioning on all the training data $y_{1:n}$.

**Q14:** Run a Kalman smoother to compute smoothed estimates of the signal, $\mathbb{E}[s_t | y_{1:n}]$, conditionally on all the _training data_. Then, similarly to above, plot the following:
1. The data $y_{1:n+m}$,
2. The smoothed estimates $\mathbb{E}[s_t | y_{1:n}] \pm 1$ standard deviation for the training data, i.e., $t = 1,...,n$,
3. The predictions $\mathbb{E}[s_t | y_{1:n}] \pm 1$ standard deviation for the validation data, i.e., $t = n+1,...,n+m$,
4. A vertical line indicating the switch between training and validation data, using, e.g., `plt.axvline(x=u[n])`.

_Hint:_ Express $s_t$ in terms of $\alpha_t$. Based on this expression, compute the smoothed mean and variance of $s_t$ based on the smoothed mean and covariance of $\alpha_t$.

```python
# Run the kalman smoother using updated model
y_new = np.concatenate((y[:n], np.repeat(np.nan, m)))
ks = kalman_smoother(y=y_new, model=model, kf=kf)

# Get smoothed estimates of the signal
y_pred = np.ndarray(ndata)
for t in range(ndata):
    y_pred[t] = model.Z.dot(ks.alpha_sm[:, :, t])
```

```python
# Plot the data and smoothed estimates
plt.plot(u, y, label='Data')
plt.plot(u, y_pred, label='Predictions')
plt.plot(u, y_pred + np.sqrt(ks.eps_var), label='Upper CI', ls='-.')
plt.plot(u, y_pred - np.sqrt(ks.eps_var), label='Lower CI', ls='-.')
plt.axvline(x=u[n], color='grey', ls='--', label='Switch between training and validation data')
plt.title('Predictions')
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
```

```python
# Log-likelihood of training data based on kalman smoother predictions using updated model
llik = log_likelihood(y=y[:n], y_hat=y_pred[:n], F=ks.eps_var[:n])
print("The log-likelihood of the training data : ", llik)

# Log-likelihood of validation data based on kalman smoother predictions using updated model
llik = log_likelihood(y=y[n:], y_hat=y_pred[n:], F=ks.eps_var[n:])
print("The log-likelihood of the validation data : ", llik)
```

    The log-likelihood of the training data :  -4488.68214451553
    The log-likelihood of the validation data :  -1184.1845001726308

**Q15:** Explain, using a few sentences, the qualitative differences (or similarities) between the Kalman filter predictions plotted in **Q13** and the smoothed signal estimates plotted in **Q14** for,
1. Training data points, $t \leq n$
2. Validation data points, $t > n$

**A:** The kalman smoother gives better starting predictions than the kalman filter. But we noticed a decrease in the log-likelihood of the smoother when compared to the filter on the training data. This result was a bit surprising. We expected better estimates from smoothing since it uses all the data and adds a backward pass on top of filtering.

However, both the smoothed and filter estimates are quite similar on the validation data, with the log-likelihood of kalman smoother slightly better than that of the kalman filter on validation data. The confidence intervals were also found to be tighter for the smoothed estimates.

We can shed additional light on the properties of the process under study by further decomposing the signal into its trend and seasonal components.

**Q16:** Using the results of the state smoother, compute and plot the _smoothed estimates_ of the two signal components, i.e.:

1. Trend: $\hat \mu_{t|n} = \mathbb{E}[\mu_t | y_{1:n}]$ for $t = 1,\dots,n$
2. Seasonal: $\hat \gamma_{t|n} = \mathbb{E}[\gamma_t | y_{1:n}]$ for $t = 1,\dots,n$

_(You don't have to include confidence intervals here if don't want to, for brevity.)_

```python
# Compute the trend and seasonal components from smoothed mean of states
mu_hat    = np.ndarray(ndata)
gamma_hat = np.ndarray(ndata)

# Compute
for t in range(ndata):
    mu_hat[t]    = Z_mu.dot(ks.alpha_sm[:k, :, t])
    gamma_hat[t] = Z_gamma.dot(ks.alpha_sm[k:, :, t])

# Plot the data and smoothed estimates of trend and seasonal components
plt.plot(u, y, label='Data')
plt.plot(u, mu_hat, label='Trend component estimates')
plt.plot(u, gamma_hat, label='Seasonal component estimates')
plt.axvline(x=u[n], color='grey', ls='--', label='Switch between training and validation data')
plt.title('Predictions')
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
```

## 2.5 Missing data
We conclude this section by illustrating one of the key merits of the state space approach to time series analysis, namely the simplicity of handling missing data. To this end we will assume that a chunk of observations in the middle of the training data is missing.

**Q17:** Let the values $y_{t}$ for $ 300 < t \leq 400$ be missing (set to `np.nan`). Modify the data and rerun the Kalman filter and smoother. Plot,
1. The Kalman filter predictions, analogously to **Q8**
2. The Kalman smoother predictions, analogously to **Q13**

Comment on the qualitative differences between the filter and smoother estimates and explain what you see (in a couple of sentences).

```python
# Set observations from t=300 to t=400 to NaNs and predict using kalman filter and smoother
y_new = np.concatenate((y[:300], np.repeat(np.nan, 100), y[400:n], np.repeat(np.nan, m)))

# Kalman filter predictions
kf = kalman_filter(y=y_new, model=model)

# Kalman smoother predictions
ks = kalman_smoother(y=y_new, model=model, kf=kf)
ks_y_pred = np.ndarray(len(y_new))
for t in range(len(y_new)):
    ks_y_pred[t] = model.Z.dot(ks.alpha_sm[:, :, t])

# Plot the missing data along with the kalman filter predictions and smoother predictions
plt.plot(u[:300], y[:300], label='Data', color='royalblue')
plt.plot(u[300:400], y[300:400], label='Missing data', color='grey')
plt.plot(u[400:n], y[400:n], color='royalblue')
plt.plot(u[n:], y[n:], color='grey')
plt.plot(u, kf.y_pred, label='Kalman filter predictions', color='green')
plt.plot(u, ks_y_pred, label='Kalman smoother predictions', color='red')
plt.axvline(x=u[n], color='grey', ls='--', label='Switch between training and validation data')
plt.title('Predictions')
plt.legend()
plt.xlabel('Year')
plt.ylabel('GMSL deviation (mm)')
```

We see that the filter and smoother predictions are almost exactly the same on the missing validation data (we observed a slight increase in the log-likelihood of the the smoother compared to the filter). Since we are essentially performing online predictions i.e. predicting the current time point based on all the previous time points upto the current one, kalman filtering is sufficient in this scenario. Furthermore, smoothing would achieve the same results as the filtering because the additional backpass through the data in smoothing does not add any new information and hence the quality of predictions is the same.

## Appendix A. Why didn't we learn the trend noise variance as well?
In the assignment above we have fixed $\sigma_\mu$ to a small value. Conceptually it would have been straightforward to learn also this parameter with the EM algorithm. However, unfortunately, the maximum likelihood estimate of $\sigma_\mu$ often ends up being too large to result in accurate _long term predictions_. The reason for this issue is that the structural model

\begin{align}
    y_t = \mu_t + \gamma_t + \varepsilon_t
\end{align}

is not a perfect description of reality. As a consequence, when learning the parameters the mismatch between the model and the data is compensated for by increasing the noise variances. This results in a trend component which does not only capture the long term trends of the data, but also seemingly random variations due to a model misspecification, possibly resulting in poor _long range predictions_.

Kitagawa (Introduction to Time Series Modeling, CRC Press, 2010, Section 12.3) discusses this issue and proposes two solutions. The first is a simple and pragmatic one: simply fix $\sigma_\mu^2$ to a value smaller than the maximum likelihood estimate. This is the approach we have taken in this assignment. The issue is of course that in practice it is hard to know what value to pick, which boild down to manual trial and error (or, if you are lucky, the designer of the lab assignment will tell you which value to use!).

The second, more principled, solution proposed by Kitagawa is to augment the model with a stationary AR component as well. That is, we model

\begin{align}
    y_t = \mu_t + \gamma_t + \nu_t + \varepsilon_t
\end{align}

where $\nu_t \sim$ AR$(p)$. By doing so, the stationary AR component can compensate for the discrepancies between the original structural model and the "true data generating process". It is straightforward to include this new component in the state space representation (how?) and to run the Kalman filter and smoother on the resulting model. Indeed, this is one of the beauties with working with the state space representation of time series data! However, the M-step of the EM algorithm becomes a bit more involved if we want to use the method to estimate also the AR coefficients of the $\nu$-component, which is beyond the scope of this lab assignment.

ext: ipynb

max_stars_repo_path: LGSS and Kalman Filter/Code.ipynb

max_issues_repo_path: LGSS and Kalman Filter/Code.ipynb

max_forks_repo_path: LGSS and Kalman Filter/Code.ipynb


text:
# Robust Kalman Filtering for Path Retrieval

#### Vinícius Miranda

Please consult the accompanying report for a fuller explanation of the goal of this notebook.

_References:_ [Robust Kalman filtering for vehicle tracking](https://www.cvxpy.org/examples/applications/robust_kalman.html)

```python
import numpy as np
import matplotlib.pyplot as plt
import cvxpy as cvx
```

### Simulating the data

Across all experimentts, we simulate 1000 time steps for a total duration of 60 seconds. The damping factor is set to 0.05. We shall vary the parametrizations of the unknown input vector and of the measurement noise.

_Experiment 1_

$$\textbf{u}(k)\sim \mathcal{N}(0, I_2)$$

$$
\begin{equation}
  \boldsymbol{\eta}(k) \sim P(X),\text{ where }
  P(X)=
  \begin{cases}
    \mathcal{N}(0, I_2), & \text{if}\ \rho > 0.2 \\
    \mathcal{N}(0, 15I_2), & \text{otherwise}
  \end{cases}
\end{equation}
$$
$$\rho \sim \mathcal{U}(0, 1)$$

_Experiment 2_

$$\textbf{u}(k)\sim \mathcal{N}(0, I_2)$$

$$
\begin{equation}
  \boldsymbol{\eta}(k) \sim P(X),\text{ where }
  P(X)=
  \begin{cases}
    \mathcal{N}(0, I_2), & \text{if}\ \rho > 0.4 \\
    \mathcal{N}(0, 15I_2), & \text{otherwise}
  \end{cases}
\end{equation}
$$

$$\rho \sim \mathcal{U}(0, 1)$$

_Experiment 3_

$$\textbf{u}(k)\sim \mathcal{N}(0, 3I_2)$$

$$
\begin{equation}
  \boldsymbol{\eta}(k) \sim P(X),\text{ where }
  P(X)=
  \begin{cases}
    \mathcal{N}(0, 3I_2), & \text{if}\ \rho > 0.4 \\
    \mathcal{N}(0, 15I_2), & \text{otherwise}
  \end{cases}
\end{equation}
$$

$$\rho \sim \mathcal{U}(0, 1)$$

```python
n_experiments = 3
n = 1000 # timesteps
T = 60 # total duration
ts, delta = np.linspace(0,T, n, endpoint=True, retstep=True)
gamma = .05 # damping factor

A = np.zeros((4,4))
B = np.zeros((4,2))
C = np.zeros((2,4))

A[0,0] = 1
A[1,1] = 1
A[0,2] = (1-gamma*delta/2)*delta
A[1,3] = (1-gamma*delta/2)*delta
A[2,2] = 1 - gamma*delta
A[3,3] = 1 - gamma*delta

B[0,0] = delta**2/2
B[1,1] = delta**2/2
B[2,0] = delta
B[3,1] = delta

C[0,0] = 1
C[1,1] = 1
```

```python
np.random.seed(164)
outlier_ps      = [0.2, 0.4, 0.4]
input_sigmas    = [1, 1, 3]
noise_sigmas    = [1, 1, 3]
outliers_sigmas = [15, 15, 15]

positions, measurements, inputs = [], [], []
for i, (p, input_s, noise_s, out_s) in enumerate(zip(outlier_ps, input_sigmas, noise_sigmas, outliers_sigmas)):

    positions.append(np.zeros((4,n+1)))
    positions[i][:,0] = [0,0,0,0]
    measurements.append(np.zeros((2,n)))

    # generate random input and noise vectors
    inputs.append(input_s*np.random.randn(2,n))
    noise = noise_s*np.random.randn(2,n)

    # add outliers to noise
    inds = np.random.rand(n) <= p
    noise[:,inds] = out_s*np.random.randn(2,n)[:,inds]

    # simulate the system forward in time
    for t in range(n):
        measurements[i][:,t] = C@positions[i][:,t] + noise[:,t]
        positions[i][:,t+1]  = A@positions[i][:,t] + B@inputs[i][:,t]

```

```python
def experiment_plot_helper(true, observed, axis):

    colors = ['ro', 'bo', 'go']
    n_cols = len(true)
    fig, ax = plt.subplots(2, 3, figsize=(15, 10))

    for i in range(n_cols):
        ax[0,i].plot(true[i][0,:], true[i][1,:], colors[i], alpha=.1)
        ax[1,i].plot(observed[i][0,:], observed[i][1,:], colors[i], alpha=.1)
        ax[0,i].axis(axis[i])
        ax[1,i].axis(axis[i])
        ax[0,i].set_title("Experiment {}: True".format(i+1))
        ax[1,i].set_title("Experiment {}: Measured".format(i+1))

    plt.show()
```

```python
experiment_plot_helper(positions, measurements, [[-1, 16, -60, 5], [-2, 33, -6, 1], [-55, 3, -60, 12]])
```

```python
# We copy the true trajectory since we will now predict it
positions_true = [p.copy() for p in positions]

def simulator(filtering, tau = .1, tau_robust = 2, rho = 2):

    if filtering not in ["Kalman", "Robust Kalman"]:
        raise("Invalid filter type.")

    # Parameters
    tau = tau if filtering == "Kalman" else tau_robust

    positions_pred, inputs_pred = [], []
    for i in range(n_experiments):

        # Optimization variables
        positions = cvx.Variable(shape=(4, n+1))
        inputs    = cvx.Variable(shape=(2, n))
        noise     = cvx.Variable(shape=(2, n))

        # Objective
        obj = cvx.Minimize(cvx.sum_squares(inputs) + tau*cvx.sum_squares(noise)) \
              if filtering == "Kalman" else \
              cvx.Minimize(cvx.sum_squares(inputs) + \
                           cvx.sum([tau*cvx.huber(cvx.norm(noise[:,t]), rho) for t in range(n)]))

        constr = []
        for t in range(n):
            constr += [ positions[:,t+1]       == A*positions[:,t] + B*inputs[:,t] ,
                        measurements[i][:,t]   == C*positions[:,t] + noise[:,t]   ]

        cvx.Problem(obj, constr).solve()

        positions_pred.append(np.array(positions.value))
    return positions_pred
```

```python
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)

pos_results = {}
parametrizations = [("Kalman", 0.01), ("Kalman", 0.1), ("Kalman", 1), ("Kalman", 10),
                    ("Robust Kalman", None, 1, 1), ("Robust Kalman", None, 1, 10),
                    ("Robust Kalman", None, 10, 1), ("Robust Kalman", None, 10, 10)]

for par in parametrizations:
    pos_results[par] = simulator(*par)
```

```python
fig, ax = plt.subplots(6, 4, figsize=(20, 30))
clr = ['r', 'b', 'g']

for k, ax_lim in enumerate([[-1, 16, -60, 5], [-2, 33, -6, 1], [-55, 3, -60, 12]]):
    for i, par in enumerate(parametrizations):

        if i < 4:
            ax[2*k, i].plot(positions_true[k][0,:], positions_true[k][1,:], clr[k]+'-')
            ax[2*k, i].plot(measurements[k][0,:], measurements[k][1,:], clr[k]+'o', alpha=.05)
            ax[2*k, i].plot(pos_results[par][k][0,:], pos_results[par][k][1,:], 'k--', alpha=.8)
            ax[2*k, i].axis(ax_lim)
            ax[2*k, i].set_title(r"Normal Filter. $\tau = {:.2f}$".format(par[1]))

        else:
            ax[2*k+1, i - 4].plot(positions_true[k][0,:], positions_true[k][1,:], clr[k]+'-')
            ax[2*k+1, i - 4].plot(measurements[k][0,:], measurements[k][1,:], clr[k]+'o', alpha=.05)
            ax[2*k+1, i - 4].plot(pos_results[par][k][0,:], pos_results[par][k][1,:], 'k--', alpha=.8)
            ax[2*k+1, i - 4].axis(ax_lim)
            ax[2*k+1, i - 4].set_title(r"Robust Filter. $\tau = {:.2f}$, $\varrho = {:.2f}$".format(par[2], par[3]))

plt.show()
```

max_stars_repo_path: Robust Kalman Filtering for Path Retrieval.ipynb


text:
# Data Analysis and Machine Learning: Elements of machine learning

**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University

Date: **Sep 28, 2018**

Copyright 1999-2018, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license

## Neural networks

Artificial neural networks are computational systems that can learn to
perform tasks by considering examples, generally without being
programmed with any task-specific rules. It is supposed to mimic a
biological system, wherein neurons interact by sending signals in the
form of mathematical functions between layers. All layers can contain
an arbitrary number of neurons, and each connection is represented by
a weight variable.

## Artificial neurons

The field of artificial neural networks has a long history of
development, and is closely connected with the advancement of computer
science and computers in general. A model of artificial neurons was
first developed by McCulloch and Pitts in 1943 to study signal
processing in the brain and has later been refined by others. The
general idea is to mimic neural networks in the human brain, which is
composed of billions of neurons that communicate with each other by
sending electrical signals.  Each neuron accumulates its incoming
signals, which must exceed an activation threshold to yield an
output. If the threshold is not overcome, the neuron remains inactive,
i.e. has zero output.

This behaviour has inspired a simple mathematical model for an artificial neuron.

$$
\begin{equation}
 y = f\left(\sum_{i=1}^n w_ix_i\right) = f(u)
\label{artificialNeuron} \tag{1}
\end{equation}
$$

Here, the output $y$ of the neuron is the value of its activation function, which have as input
a weighted sum of signals $x_i, \dots ,x_n$ received by $n$ other neurons.

Conceptually, it is helpful to divide neural networks into four
categories:
1. general purpose neural networks for supervised learning,

2. neural networks designed specifically for image processing, the most prominent example of this class being Convolutional Neural Networks (CNNs),

3. neural networks for sequential data such as Recurrent Neural Networks (RNNs), and

4. neural networks for unsupervised learning such as Deep Boltzmann Machines.

In natural science, DNNs and CNNs have already found numerous applications. In
statistical physics, they have been applied to detect phase
transitions in 2D Ising and Potts models, lattice gauge theories, and
different phases of polymers, or solving the Navier-Stokes equation in weather forecasting.
Deep learning has also found interesting applications in quantum
physics. Various quantum phase transitions can be detected and studied
using DNNs and CNNs,
topological phases, and even non-equilibrium many-body
localization. Representing quantum states as DNNs quantum state
tomography are among some of the impressive
achievements to reveal the potential of DNNs to facilitate the study
of quantum systems.

In quantum information theory, it has been shown that one can perform
gate decompositions with the help of neural. In lattice quantum chromodynamics,
DNNs have been used to learn action parameters in regions of parameter
space where PCA fails.

The applications are not limited to the natural sciences. There is a plethora of applications in essentially all disciplines, from the humanities to life science and medicine.

## Neural network types

An artificial neural network (NN), is a computational model that
consists of layers of connected neurons, or *nodes*.  It is supposed
to mimic a biological nervous system by letting each neuron interact
with other neurons by sending signals in the form of mathematical
functions between layers.  A wide variety of different NNs have been
developed, but most of them consist of an input layer, an output layer
and eventual layers in-between, called *hidden layers*. All layers can
contain an arbitrary number of nodes, and each connection between two
nodes is associated with a weight variable.

Neural networks (also called neural nets) are neural-inspired
nonlinear models for supervised learning.  As we will see, neural nets
can be viewed as natural, more powerful extensions of supervised
learning methods such as linear and logistic regression and soft-max
methods.

## Feed-forward neural networks

The feed-forward neural network (FFNN) was the first and simplest type of NN devised. In this network,
the information moves in only one direction: forward through the layers.

Nodes are represented by circles, while the arrows display the connections between the nodes, including the
direction of information flow. Additionally, each arrow corresponds to a weight variable, not displayed here.
We observe that each node in a layer is connected to *all* nodes in the subsequent layer,
making this a so-called *fully-connected* FFNN.

A different variant of FFNNs are *convolutional neural networks* (CNNs), which have a connectivity pattern
inspired by the animal visual cortex. Individual neurons in the visual cortex only respond to stimuli from
small sub-regions of the visual field, called a receptive field. This makes the neurons well-suited to exploit the strong
spatially local correlation present in natural images. The response of each neuron can be approximated mathematically
as a convolution operation.

CNNs emulate the behaviour of neurons in the visual cortex by enforcing a *local* connectivity pattern
between nodes of adjacent layers: Each node
in a convolutional layer is connected only to a subset of the nodes in the previous layer,
in contrast to the fully-connected FFNN.
Often, CNNs
consist of several convolutional layers that learn local features of the input, with a fully-connected layer at the end,
which gathers all the local data and produces the outputs. They have wide applications in image and video recognition

## Recurrent neural networks

So far we have only mentioned NNs where information flows in one direction: forward. *Recurrent neural networks* on
the other hand, have connections between nodes that form directed *cycles*. This creates a form of
internal memory which are able to capture information on what has been calculated before; the output is dependent
on the previous computations. Recurrent NNs make use of sequential information by performing the same task for
every element in a sequence, where each element depends on previous elements. An example of such information is
sentences, making recurrent NNs especially well-suited for handwriting and speech recognition.

## Other types of networks

There are many other kinds of NNs that have been developed. One type that is specifically designed for interpolation
in multidimensional space is the radial basis function (RBF) network. RBFs are typically made up of three layers:
an input layer, a hidden layer with non-linear radial symmetric activation functions and a linear output layer (''linear'' here
means that each node in the output layer has a linear activation function). The layers are normally fully-connected and
there are no cycles, thus RBFs can be viewed as a type of fully-connected FFNN. They are however usually treated as
a separate type of NN due the unusual activation functions.

## Multilayer perceptrons

One uses often so-called fully-connected feed-forward neural networks
with three or more layers (an input layer, one or more hidden layers
and an output layer) consisting of neurons that have non-linear
activation functions.

Such networks are often called *multilayer perceptrons* (MLPs)

## Why multilayer perceptrons?

According to the *Universal approximation theorem*, a feed-forward neural network with just a single hidden layer containing
a finite number of neurons can approximate a continuous multidimensional function to arbitrary accuracy,
assuming the activation function for the hidden layer is a **non-constant, bounded and monotonically-increasing continuous function**.

Note that the requirements on the activation function only applies to
the hidden layer, the output nodes are always assumed to be linear, so
as to not restrict the range of output values.

We note that this theorem is only applicable to an NN with *one* hidden
layer.  Therefore, we can easily construct an NN that employs
activation functions which do not satisfy the above requirements, as
long as we have at least one layer with activation functions that
*do*. Furthermore, although the universal approximation theorem lays
the theoretical foundation for regression with neural networks, it
does not say anything about how things work in practice: A neural
network can still be able to approximate a given function reasonably
well without having the flexibility to fit *all other* functions.

## Mathematical model

$$
\begin{equation}
 y = f\left(\sum_{i=1}^n w_ix_i + b_i\right) = f(u)
\label{artificialNeuron2} \tag{2}
\end{equation}
$$

In an FFNN of such neurons, the *inputs* $x_i$ are the *outputs* of
the neurons in the preceding layer. Furthermore, an MLP is
fully-connected, which means that each neuron receives a weighted sum
of the outputs of *all* neurons in the previous layer.

## Mathematical model

First, for each node $i$ in the first hidden layer, we calculate a weighted sum $u_i^1$ of the input coordinates $x_j$,

$$
\begin{equation}
 u_i^1 = \sum_{j=1}^2 w_{ij}^1 x_j  + b_i^1
\label{_auto1} \tag{3}
\end{equation}
$$

This value is the argument to the activation function $f_1$ of each neuron $i$,
producing the output $y_i^1$ of all neurons in layer 1,

$$
\begin{equation}
 y_i^1 = f_1(u_i^1) = f_1\left(\sum_{j=1}^2 w_{ij}^1 x_j  + b_i^1\right)
\label{outputLayer1} \tag{4}
\end{equation}
$$

where we assume that all nodes in the same layer have identical
activation functions, hence the notation $f_l$

$$
\begin{equation}
 y_i^l = f_l(u_i^l) = f_l\left(\sum_{j=1}^{N_{l-1}} w_{ij}^l y_j^{l-1} + b_i^l\right)
\label{generalLayer} \tag{5}
\end{equation}
$$

where $N_l$ is the number of nodes in layer $l$. When the output of
all the nodes in the first hidden layer are computed, the values of
the subsequent layer can be calculated and so forth until the output
is obtained.

## Mathematical model

The output of neuron $i$ in layer 2 is thus,

$$
\begin{equation}
 y_i^2 = f_2\left(\sum_{j=1}^3 w_{ij}^2 y_j^1 + b_i^2\right)
\label{_auto2} \tag{6}
\end{equation}
$$

$$
\begin{equation}
 = f_2\left[\sum_{j=1}^3 w_{ij}^2f_1\left(\sum_{k=1}^2 w_{jk}^1 x_k + b_j^1\right) + b_i^2\right]
\label{outputLayer2} \tag{7}
\end{equation}
$$

where we have substituted $y_m^1$ with. Finally, the NN output yields,

$$
\begin{equation}
 y_1^3 = f_3\left(\sum_{j=1}^3 w_{1m}^3 y_j^2 + b_1^3\right)
\label{_auto3} \tag{8}
\end{equation}
$$

$$
\begin{equation}
 = f_3\left[\sum_{j=1}^3 w_{1j}^3 f_2\left(\sum_{k=1}^3 w_{jk}^2 f_1\left(\sum_{m=1}^2 w_{km}^1 x_m + b_k^1\right) + b_j^2\right)
  + b_1^3\right]
\label{_auto4} \tag{9}
\end{equation}
$$

## Mathematical model

We can generalize this expression to an MLP with $l$ hidden
layers. The complete functional form is,

$$
\begin{equation}
y^{l+1}_1\! = \!f_{l+1}\!\left[\!\sum_{j=1}^{N_l}\! w_{1j}^3 f_l\!\left(\!\sum_{k=1}^{N_{l-1}}\! w_{jk}^2 f_{l-1}\!\left(\!
 \dots \!f_1\!\left(\!\sum_{n=1}^{N_0} \!w_{mn}^1 x_n\! + \!b_m^1\!\right)
 \!\dots \!\right) \!+ \!b_k^2\!\right)
 \!+ \!b_1^3\!\right]
\label{completeNN} \tag{10}
\end{equation}
$$

which illustrates a basic property of MLPs: The only independent
variables are the input values $x_n$.

## Mathematical model

This confirms that an MLP, despite its quite convoluted mathematical
form, is nothing more than an analytic function, specifically a
mapping of real-valued vectors $\vec{x} \in \mathbb{R}^n \rightarrow
\vec{y} \in \mathbb{R}^m$.  In our example, $n=2$ and
$m=1$. Consequentially, the number of input and output values of the
function we want to fit must be equal to the number of inputs and
outputs of our MLP.

Furthermore, the flexibility and universality of a MLP can be
illustrated by realizing that the expression is essentially a nested
sum of scaled activation functions of the form

$$
\begin{equation}
 h(x) = c_1 f(c_2 x + c_3) + c_4
\label{_auto5} \tag{11}
\end{equation}
$$

where the parameters $c_i$ are weights and biases. By adjusting these
parameters, the activation functions can be shifted up and down or
left and right, change slope or be rescaled which is the key to the
flexibility of a neural network.

### Matrix-vector notation

We can introduce a more convenient notation for the activations in a NN.

Additionally, we can represent the biases and activations
as layer-wise column vectors $\vec{b}_l$ and $\vec{y}_l$, so that the $i$-th element of each vector
is the bias $b_i^l$ and activation $y_i^l$ of node $i$ in layer $l$ respectively.

We have that $\mathrm{W}_l$ is a $N_{l-1} \times N_l$ matrix, while $\vec{b}_l$ and $\vec{y}_l$ are $N_l \times 1$ column vectors.
With this notation, the sum in  becomes a matrix-vector multiplication, and we can write
the equation for the activations of hidden layer 2 in

$$
\begin{equation}
 \vec{y}_2 = f_2(\mathrm{W}_2 \vec{y}_{1} + \vec{b}_{2}) =
 f_2\left(\left[\begin{array}{ccc}
    w^2_{11} &w^2_{12} &w^2_{13} \\
    w^2_{21} &w^2_{22} &w^2_{23} \\
    w^2_{31} &w^2_{32} &w^2_{33} \\
    \end{array} \right] \cdot
    \left[\begin{array}{c}
           y^1_1 \\
           y^1_2 \\
           y^1_3 \\
          \end{array}\right] +
    \left[\begin{array}{c}
           b^2_1 \\
           b^2_2 \\
           b^2_3 \\
          \end{array}\right]\right).
\label{_auto6} \tag{12}
\end{equation}
$$

### Matrix-vector notation  and activation

The activation of node $i$ in layer 2 is

$$
\begin{equation}
 y^2_i = f_2\Bigr(w^2_{i1}y^1_1 + w^2_{i2}y^1_2 + w^2_{i3}y^1_3 + b^2_i\Bigr) =
 f_2\left(\sum_{j=1}^3 w^2_{ij} y_j^1 + b^2_i\right).
\label{_auto7} \tag{13}
\end{equation}
$$

This is not just a convenient and compact notation, but also a useful
and intuitive way to think about MLPs: The output is calculated by a
series of matrix-vector multiplications and vector additions that are
used as input to the activation functions. For each operation
$\mathrm{W}_l \vec{y}_{l-1}$ we move forward one layer.

### Activation functions

A property that characterizes a neural network, other than its
connectivity, is the choice of activation function(s).  As described
in, the following restrictions are imposed on an activation function
for a FFNN to fulfill the universal approximation theorem

  * Non-constant

  * Bounded

  * Monotonically-increasing

  * Continuous

### Activation functions, Logistic and Hyperbolic ones

The second requirement excludes all linear functions. Furthermore, in
a MLP with only linear activation functions, each layer simply
performs a linear transformation of its inputs.

Regardless of the number of layers, the output of the NN will be
nothing but a linear function of the inputs. Thus we need to introduce
some kind of non-linearity to the NN to be able to fit non-linear
functions Typical examples are the logistic *Sigmoid*

$$
\begin{equation}
 f(x) = \frac{1}{1 + e^{-x}},
\label{sigmoidActivationFunction} \tag{14}
\end{equation}
$$

and the *hyperbolic tangent* function

$$
\begin{equation}
 f(x) = \tanh(x)
\label{tanhActivationFunction} \tag{15}
\end{equation}
$$

### Relevance

The *sigmoid* function are more biologically plausible because the
output of inactive neurons are zero. Such activation function are
called *one-sided*. However, it has been shown that the hyperbolic
tangent performs better than the sigmoid for training MLPs.  has
become the most popular for *deep neural networks*

```python
%matplotlib inline

"""The sigmoid function (or the logistic curve) is a
function that takes any real number, z, and outputs a number (0,1).
It is useful in neural networks for assigning weights on a relative scale.
The value z is the weighted sum of parameters involved in the learning algorithm."""

import numpy
import matplotlib.pyplot as plt
import math as mt

z = numpy.arange(-5, 5, .1)
sigma_fn = numpy.vectorize(lambda z: 1/(1+numpy.exp(-z)))
sigma = sigma_fn(z)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(z, sigma)
ax.set_ylim([-0.1, 1.1])
ax.set_xlim([-5,5])
ax.grid(True)
ax.set_xlabel('z')
ax.set_title('sigmoid function')

plt.show()

"""Step Function"""
z = numpy.arange(-5, 5, .02)
step_fn = numpy.vectorize(lambda z: 1.0 if z >= 0.0 else 0.0)
step = step_fn(z)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(z, step)
ax.set_ylim([-0.5, 1.5])
ax.set_xlim([-5,5])
ax.grid(True)
ax.set_xlabel('z')
ax.set_title('step function')

plt.show()

"""Sine Function"""
z = numpy.arange(-2*mt.pi, 2*mt.pi, 0.1)
t = numpy.sin(z)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(z, t)
ax.set_ylim([-1.0, 1.0])
ax.set_xlim([-2*mt.pi,2*mt.pi])
ax.grid(True)
ax.set_xlabel('z')
ax.set_title('sine function')

plt.show()

"""Plots a graph of the squashing function used by a rectified linear
unit"""
z = numpy.arange(-2, 2, .1)
zero = numpy.zeros(len(z))
y = numpy.max([zero, z], axis=0)

fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(z, y)
ax.set_ylim([-2.0, 2.0])
ax.set_xlim([-2.0, 2.0])
ax.grid(True)
ax.set_xlabel('z')
ax.set_title('Rectified linear unit')

plt.show()
```

## Setting up a Multi-layer perceptron model

```python
from scipy import optimize

class Neural_Network(object):
    def __init__(self, Lambda=0):
        #Define Hyperparameters
        self.inputLayerSize = 2
        self.outputLayerSize = 1
        self.hiddenLayerSize = 3

        #Weights (parameters)
        self.W1 = np.random.randn(self.inputLayerSize,self.hiddenLayerSize)
        self.W2 = np.random.randn(self.hiddenLayerSize,self.outputLayerSize)

        #Regularization Parameter:
        self.Lambda = Lambda

    def forward(self, X):
        #Propogate inputs though network
        self.z2 = np.dot(X, self.W1)
        self.a2 = self.sigmoid(self.z2)
        self.z3 = np.dot(self.a2, self.W2)
        yHat = self.sigmoid(self.z3)
        return yHat

    def sigmoid(self, z):
        #Apply sigmoid activation function to scalar, vector, or matrix
        return 1/(1+np.exp(-z))

    def sigmoidPrime(self,z):
        #Gradient of sigmoid
        return np.exp(-z)/((1+np.exp(-z))**2)

    def costFunction(self, X, y):
        #Compute cost for given X,y, use weights already stored in class.
        self.yHat = self.forward(X)
        J = 0.5*sum((y-self.yHat)**2)/X.shape[0] + (self.Lambda/2)*(np.sum(self.W1**2)+np.sum(self.W2**2))
        return J

    def costFunctionPrime(self, X, y):
        #Compute derivative with respect to W and W2 for a given X and y:
        self.yHat = self.forward(X)

        delta3 = np.multiply(-(y-self.yHat), self.sigmoidPrime(self.z3))
        #Add gradient of regularization term:
        dJdW2 = np.dot(self.a2.T, delta3)/X.shape[0] + self.Lambda*self.W2

        delta2 = np.dot(delta3, self.W2.T)*self.sigmoidPrime(self.z2)
        #Add gradient of regularization term:
        dJdW1 = np.dot(X.T, delta2)/X.shape[0] + self.Lambda*self.W1

        return dJdW1, dJdW2

    #Helper functions for interacting with other methods/classes
    def getParams(self):
        #Get W1 and W2 Rolled into vector:
        params = np.concatenate((self.W1.ravel(), self.W2.ravel()))
        return params

    def setParams(self, params):
        #Set W1 and W2 using single parameter vector:
        W1_start = 0
        W1_end = self.hiddenLayerSize*self.inputLayerSize
        self.W1 = np.reshape(params[W1_start:W1_end], \
                             (self.inputLayerSize, self.hiddenLayerSize))
        W2_end = W1_end + self.hiddenLayerSize*self.outputLayerSize
        self.W2 = np.reshape(params[W1_end:W2_end], \
                             (self.hiddenLayerSize, self.outputLayerSize))

    def computeGradients(self, X, y):
        dJdW1, dJdW2 = self.costFunctionPrime(X, y)
        return np.concatenate((dJdW1.ravel(), dJdW2.ravel()))

class trainer(object):
    def __init__(self, N):
        #Make Local reference to network:
        self.N = N

    def callbackF(self, params):
        self.N.setParams(params)
        self.J.append(self.N.costFunction(self.X, self.y))
        self.testJ.append(self.N.costFunction(self.testX, self.testY))

    def costFunctionWrapper(self, params, X, y):
        self.N.setParams(params)
        cost = self.N.costFunction(X, y)
        grad = self.N.computeGradients(X,y)
        return cost, grad

    def train(self, trainX, trainY, testX, testY):
        #Make an internal variable for the callback function:
        self.X = trainX
        self.y = trainY

        self.testX = testX
        self.testY = testY

        #Make empty list to store training costs:
        self.J = []
        self.testJ = []

        params0 = self.N.getParams()

        options = {'maxiter': 200, 'disp' : True}
        _res = optimize.minimize(self.costFunctionWrapper, params0, jac=True, method='BFGS', \
                                 args=(trainX, trainY), options=options, callback=self.callbackF)

        self.N.setParams(_res.x)
        self.optimizationResults = _res
```

## Two-layer Neural Network

```python
import numpy as np

#sigmoid
def nonlin(x, deriv=False):
    if (deriv==True):
        return x*(1-x)
    return 1/(1+np.exp(-x))

#input data
x=np.array([[0,0,1],[0,1,1],[1,0,1],[1,1,1]])

#output data
y=np.array([0,1,1,0]).T

#seed random numbers to make calculation
np.random.seed(1)

#initialize weights with mean=0
syn0=2*np.random.random((3,4))-1

for iter in range(10000):
    #forward propogation
    l0=x
    l1=nonlin(np.dot(l0,syn0))
    l1_error=y-l1
    #multiply error by slope of sigmoid at values of l1
    l1_delta=l1_error*nonlin(l1,True)
    #update weights
    syn0+=np.dot(l0.T, l1_delta)

print("Output after training: ",l1 )
```

    Output after training:  [[6.55109972e-03 9.93684857e-01 9.93925710e-01 6.62304973e-03]
     [1.71082162e-03 9.97516440e-01 9.97766376e-01 1.82685927e-03]
     [2.05800960e-03 9.98268211e-01 9.97548919e-01 1.77362990e-03]
     [5.35659849e-04 9.99320839e-01 9.99100767e-01 4.87503198e-04]]

```python
import numpy as np
import random
class Network(object):

    def _init_(self, sizes):
        self.num_layers=len(sizes)
        self.sizes=sizes
        self.biases=[np.random.randn(y,1) for y in sizes[1:]]
        self.weights=[np.random.randn(y,x) for x,y in zip(sizes[:-1], sizes[1:])]

#sizes is the number of neurons in each layer
#for example, say n_1st_layer=3, n_2nd_layer=3, n_3rd_layer=1, then net=Network([3,3,1])

#The biases and weights are initialized randomly, using Gaussian distributions of mean=0, stdev=1
#z is a vector (or a np.array)

    def feedforward(self,a):
        #returns output w/ 'a' as an input
        for b, w in zip(self.biases, self.weights):
            a=sigmoid(np.dot(w,b)+b)
        return a

#Apply a Stochastic Gradient Descent (SGD) method:
    def SGD(self, training_data, epochs, mini_batch_size, eta, test_data=None):
        """Trains network using batches incorporating SGD. The network will be evaluated against the
        test data after each epoch, with partial progress being printed out (this is useful for tracking,
        but slows the process.)"""
        if test_data: n_test=len(test_data)
        n=len(training_data)
        for j in xrange(epochs):
            random.shuffle(training_data)
            mini_batches=[training_data[k:k+mini_batch_size] for k in xrange(o,n,mini_batch_size)]
            for mini_batch in mini_batches:
                self.update_mini_batch(mini_batch, eta)
            if test_data:
                print ("Epoch {0}: {1}/{2}".format(j, self.evaluate(test_data), n_test))
            else:
                print ("Epoch {0} complete".format(j))

    def update_mini_batch(self, mini_batch, eta):
        #updates w and b using backpropagation to a single mini batch. eta is the learning rate."
        nabla_b=[np.zeros(b.shape) for b in self.biases]
        nabla_w=[np.zeros(w.shape) for w in self.weights]
        for x,y in mini_batch:
            delta_nabla_b, delta_nabla_w=self.backprop(x,y)
            nabla_b=[nb+dnb for nb, dnb in zip(nabla_b, delta_nabla_b)]
            nabla_w=[nw+dnw for nw, dnw in zip(nabla_w, delta_nabla_w)]
        self.weights=[w-(eta/len(mini_batch))*nw for w, nw in zip(self.weights, nabla_w)]
        self.biases=[b-(eta/len(mini_batch))*nb for b, nb in zip(self.biases, nabla_b)]

    def backprop(self, x, y):
        """Return a tuple ``(nabla_b, nabla_w)`` representing the
        gradient for the cost function C_x.  ``nabla_b`` and
        ``nabla_w`` are layer-by-layer lists of numpy arrays, similar
        to ``self.biases`` and ``self.weights``."""
        nabla_b = [np.zeros(b.shape) for b in self.biases]
        nabla_w = [np.zeros(w.shape) for w in self.weights]
        # feedforward
        activation = x
        activations = [x] # list to store all the activations, layer by layer
        zs = [] # list to store all the z vectors, layer by layer
        for b, w in zip(self.biases, self.weights):
            z = np.dot(w, activation)+b
            zs.append(z)
            activation = sigmoid(z)
            activations.append(activation)
        # backward pass
        delta = self.cost_derivative(activations[-1], y) * \
            sigmoid_prime(zs[-1])
        nabla_b[-1] = delta
        nabla_w[-1] = np.dot(delta, activations[-2].transpose())
        # Note that the variable l in the loop below is used a little
        # differently to the notation in Chapter 2 of the book.  Here,
        # l = 1 means the last layer of neurons, l = 2 is the
        # second-last layer, and so on.  It's a renumbering of the
        # scheme in the book, used here to take advantage of the fact
        # that Python can use negative indices in lists.
        for l in xrange(2, self.num_layers):
            z = zs[-l]
            sp = sigmoid_prime(z)
            delta = np.dot(self.weights[-l+1].transpose(), delta) * sp
            nabla_b[-l] = delta
            nabla_w[-l] = np.dot(delta, activations[-l-1].transpose())
        return (nabla_b, nabla_w)

    def evaluate(self, test_data):
        """Return the number of test inputs for which the neural
        network outputs the correct result. Note that the neural
        network's output is assumed to be the index of whichever
        neuron in the final layer has the highest activation."""
        test_results = [(np.argmax(self.feedforward(x)), y)
                        for (x, y) in test_data]
        return sum(int(x == y) for (x, y) in test_results)

    def cost_derivative(self, output_activations, y):
        """Return the vector of partial derivatives \partial C_x /
        \partial a for the output activations."""
        return (output_activations-y)

#Functions
def sigmoid(z):
        return 1.0/(1.0+np.exp(-z))

def sigmoid_prime(z):
    return sigmoid(z)*(1-sigmoid(z))

network=Network()
```

```python
# %load neural-networks-and-deep-learning/src/mnist_loader.py
"""
mnist_loader
~~~~~~~~~~~~

A library to load the MNIST image data.  For details of the data
structures that are returned, see the doc strings for ``load_data``
and ``load_data_wrapper``.  In practice, ``load_data_wrapper`` is the
function usually called by our neural network code.
"""

#### Libraries
# Standard library
import pickle
import gzip

# Third-party libraries
import numpy as np

def load_data():
    """Return the MNIST data as a tuple containing the training data,
    the validation data, and the test data.

    The ``training_data`` is returned as a tuple with two entries.
    The first entry contains the actual training images.  This is a
    numpy ndarray with 50,000 entries.  Each entry is, in turn, a
    numpy ndarray with 784 values, representing the 28 * 28 = 784
    pixels in a single MNIST image.

    The second entry in the ``training_data`` tuple is a numpy ndarray
    containing 50,000 entries.  Those entries are just the digit
    values (0...9) for the corresponding images contained in the first
    entry of the tuple.

    The ``validation_data`` and ``test_data`` are similar, except
    each contains only 10,000 images.

    This is a nice data format, but for use in neural networks it's
    helpful to modify the format of the ``training_data`` a little.
    That's done in the wrapper function ``load_data_wrapper()``, see
    below.
    """
    f = gzip.open('../data/mnist.pkl.gz', 'rb')
    training_data, validation_data, test_data = cPickle.load(f)
    f.close()
    return (training_data, validation_data, test_data)

def load_data_wrapper():
    """Return a tuple containing ``(training_data, validation_data,
    test_data)``. Based on ``load_data``, but the format is more
    convenient for use in our implementation of neural networks.

    In particular, ``training_data`` is a list containing 50,000
    2-tuples ``(x, y)``.  ``x`` is a 784-dimensional numpy.ndarray
    containing the input image.  ``y`` is a 10-dimensional
    numpy.ndarray representing the unit vector corresponding to the
    correct digit for ``x``.

    ``validation_data`` and ``test_data`` are lists containing 10,000
    2-tuples ``(x, y)``.  In each case, ``x`` is a 784-dimensional
    numpy.ndarry containing the input image, and ``y`` is the
    corresponding classification, i.e., the digit values (integers)
    corresponding to ``x``.

    Obviously, this means we're using slightly different formats for
    the training data and the validation / test data.  These formats
    turn out to be the most convenient for use in our neural network
    code."""
    tr_d, va_d, te_d = load_data()
    training_inputs = [np.reshape(x, (784, 1)) for x in tr_d[0]]
    training_results = [vectorized_result(y) for y in tr_d[1]]
    training_data = zip(training_inputs, training_results)
    validation_inputs = [np.reshape(x, (784, 1)) for x in va_d[0]]
    validation_data = zip(validation_inputs, va_d[1])
    test_inputs = [np.reshape(x, (784, 1)) for x in te_d[0]]
    test_data = zip(test_inputs, te_d[1])
    return (training_data, validation_data, test_data)

def vectorized_result(j):
    """Return a 10-dimensional unit vector with a 1.0 in the jth
    position and zeroes elsewhere.  This is used to convert a digit
    (0...9) into a corresponding desired output from the neural
    network."""
    e = np.zeros((10, 1))
    e[j] = 1.0
    return e

net=network.Network([784,30,30])
net.SGD(training_data,30,10,3,test_data=test_data)
```

```python

```

max_stars_repo_name: anacost/MachineLearning

max_issues_repo_name: bernharl/MachineLearning

max_forks_repo_name: bernharl/MachineLearning