text:
[Sebastian Raschka](http://sebastianraschka.com)

- [Open in IPython nbviewer](http://nbviewer.ipython.org/github/rasbt/python_reference/blob/master/tutorials/python_data_entry_point.ipynb?create=1)

- [Link to this IPython notebook on Github](https://github.com/rasbt/python_reference/blob/master/tutorials/python_data_entry_point.ipynb)

- [Link to the GitHub Repository pattern_classification](http://nbviewer.ipython.org/github/rasbt/pattern_classification/blob/master/python_howtos/)

```python
%load_ext watermark
```

```python
%watermark -a 'Gopala KR' -v -d -p numpy,scipy,matplotlib,scikit-learn
```

    Gopala KR 2018-01-28

    CPython 3.6.3
    IPython 6.2.1

    numpy 1.13.1
    scipy 0.19.1
    matplotlib 2.0.2
    scikit-learn 0.19.0

    /srv/venv/lib/python3.6/site-packages/watermark/watermark.py:155: DeprecationWarning: Importing scikit-learn as `scikit-learn` has been depracated and will not be supported anymore in v1.7.0. Please use the package name `sklearn` instead.
      DeprecationWarning)

[More information](http://nbviewer.ipython.org/github/rasbt/python_reference/blob/master/ipython_magic/watermark.ipynb) about the `watermark` magic command extension.

I would be happy to hear your comments and suggestions.
Please feel free to drop me a note via
[twitter](https://twitter.com/rasbt), [email](mailto:bluewoodtree@gmail.com), or [google+](https://plus.google.com/+SebastianRaschka).

# Entry point: Data

### - Using Python's sci-packages to prepare data for Machine Learning tasks and other data analyses

In this short tutorial I want to provide a short overview of some of my favorite Python tools for common procedures as entry points for general pattern classification and machine learning tasks, and various other data analyses.

# Sections

- [Installing Python packages](#Installing-Python-packages)

- [About the dataset](#About-the-dataset)

- [Downloading and saving CSV data files from the web](#Downloading-and-savin-CSV-data-files-from-the-web)

- [Reading in a dataset from a CSV file](#Reading-in-a-dataset-from-a-CSV-file)

- [Visualizating of a dataset](#Visualizating-of-a-data)

    - [Histograms](#Histograms)

    - [Scatterplots](#Scatterplots)

- [Splitting into training and test dataset](#Splitting-into-training-and-test-dataset)

- [Feature Scaling](#Feature-Scaling)

    - [Standardization](#Standardization)

    - [Min-Max scaling (Normalization)](#Min-Max-scaling-Normalization)

- [Linear Transformation: Principal Component Analysis (PCA)](#PCA)

- [Linear Transformation: Linear Discrciminant Analysis (LDA)](#MDA)

- [Simple Supervised Classification](#Simple-Supervised-Classification)

    - [Linear Discriminant Analysis as simple linear classifier](#Linear-Discriminant-Analysis-as-simple-linear-classifier)

    - [Classification Stochastic Gradient Descent (SGD)](#SGD)

- [Saving the processed datasets](#Saving-the-processed-datasets)

    - [Pickle](#Pickle)

    - [Comma Separated Values (CSV)](#Comma-Separated-Values)

## Installing Python packages

[[back to top]](#Sections)

**In this section want to recommend a way for installing the required Python-packages packages if you have not done so, yet. Otherwise you can skip this part.**

The packages we will be using in this tutorial are:

- [NumPy](http://www.numpy.org)
- [SciPy](http://www.scipy.org)
- [matplotlib](http://matplotlib.org)
- [scikit-learn](http://scikit-learn.org/stable/)

Although they can be installed step-by-step "manually", but I highly recommend you to take a look at the [Anaconda](https://store.continuum.io/cshop/anaconda/) Python distribution for scientific computing.

Anaconda is distributed by Continuum Analytics, but it is completely free and includes more than 195+ packages for science and data analysis as of today.
The installation procedure is nicely summarized here: http://docs.continuum.io/anaconda/install.html

If this is too much, the [Miniconda](http://conda.pydata.org/miniconda.html) might be right for you. Miniconda is basically just a Python distribution with the Conda package manager, which let's us install a list of Python packages into a specified `conda` environment from the Shell terminal, e.g.,

$[bash]> conda create -n myenv python=3
$[bash]> source activate myenv
$[bash]> conda install -n myenv numpy scipy matplotlib scikit-learn

When we start "python" in your current shell session now, it will use the Python distribution in the virtual environment "myenv" that we have just created. To un-attach the virtual environment, you can just use

$[bash]> source deactivate myenv

**Note:** environments will be created in ROOT_DIR/envs by default, you can use the `-p` instead of the `-n` flag in the conda commands above in order to specify a custom path.

**I find this procedure very convenient, especially if you are working with different distributions and versions of Python with different modules and packages installed and it is extremely useful for testing your own modules.**

## About the dataset

[[back to top]](#Sections)

For the following tutorial, we will be working with the free "Wine" Dataset that is deposited on the UCI machine learning repository
(http://archive.ics.uci.edu/ml/datasets/Wine).

**Reference:**
Forina, M. et al, PARVUS - An Extendible Package for Data
Exploration, Classification and Correlation. Institute of Pharmaceutical
and Food Analysis and Technologies, Via Brigata Salerno,
16147 Genoa, Italy.

Bache, K. & Lichman, M. (2013). UCI Machine Learning Repository [http://archive.ics.uci.edu/ml]. Irvine, CA: University of California, School of Information and Computer Science.

The Wine dataset consists of 3 different classes where each row correspond to a particular wine sample.

The class labels (1, 2, 3) are listed in the first column, and the columns 2-14 correspond to the following 13 attributes (features):

1) Alcohol
2) Malic acid
3) Ash
4) Alcalinity of ash
5) Magnesium
6) Total phenols
7) Flavanoids
8) Nonflavanoid phenols
9) Proanthocyanins
10) Color intensity
11) Hue
12) OD280/OD315 of diluted wines
13) Proline

An excerpt from the wine_data.csv dataset:

1,14.23,1.71,2.43,15.6,127,2.8,3.06,.28,2.29,5.64,1.04,3.92,1065
1,13.2,1.78,2.14,11.2,100,2.65,2.76,.26,1.28,4.38,1.05,3.4,1050
[...]
2,12.37,.94,1.36,10.6,88,1.98,.57,.28,.42,1.95,1.05,1.82,520
2,12.33,1.1,2.28,16,101,2.05,1.09,.63,.41,3.27,1.25,1.67,680
[...]
3,12.86,1.35,2.32,18,122,1.51,1.25,.21,.94,4.1,.76,1.29,630
3,12.88,2.99,2.4,20,104,1.3,1.22,.24,.83,5.4,.74,1.42,530

## Downloading and saving CSV data files from the web

[[back to top]](#Sections)

Usually, we have our data stored locally on our disk in as a common text (or CSV) file with comma-, tab-, or whitespace-separated rows. Below is just an example for how you can CSV datafile from a HTML website directly into Python and optionally save it locally.

```python
import csv
import urllib

url = 'https://raw.githubusercontent.com/rasbt/pattern_classification/master/data/wine_data.csv'
csv_cont = urllib.request.urlopen(url)
csv_cont = csv_cont.read() #.decode('utf-8')

# Optional: saving the data to your local drive
with open('./wine_data.csv', 'wb') as out:
    out.write(csv_cont)
```

**Note:** If you'd rather like to work with the data directly in `str`ing format, you could just apply the `.decode('utf-8')` method to the data that was read in byte-format by default.

## Reading in a dataset from a CSV file

[[back to top]](#Sections)

Since it is quite typical to have the input data stored locally, as mentioned above, we will use the [`numpy.loadtxt`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html) function  now to read in the data from the CSV file.
(alternatively [`np.genfromtxt()`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.genfromtxt.html) could be used in similar way, it provides some additional options)

```python
import numpy as np

# reading in all data into a NumPy array
all_data = np.loadtxt(open("./wine_data.csv","r"),
        delimiter=",",
        skiprows=0,
        dtype=np.float64
        )

# load class labels from column 1
y_wine = all_data[:,0]

# conversion of the class labels to integer-type array
y_wine = y_wine.astype(np.int64, copy=False)

# load the 14 features
X_wine = all_data[:,1:]

# printing some general information about the data
print('\ntotal number of samples (rows):', X_wine.shape[0])
print('total number of features (columns):', X_wine.shape[1])

# printing the 1st wine sample
float_formatter = lambda x: '{:.2f}'.format(x)
np.set_printoptions(formatter={'float_kind':float_formatter})
print('\n1st sample (i.e., 1st row):\nClass label: {:d}\n{:}\n'
          .format(int(y_wine[0]), X_wine[0]))

# printing the rel.frequency of the class labels
print('Class label frequencies')
print('Class 1 samples: {:.2%}'.format(list(y_wine).count(1)/y_wine.shape[0]))
print('Class 2 samples: {:.2%}'.format(list(y_wine).count(2)/y_wine.shape[0]))
print('Class 3 samples: {:.2%}'.format(list(y_wine).count(3)/y_wine.shape[0]))
```

    total number of samples (rows): 178
    total number of features (columns): 13

    1st sample (i.e., 1st row):
    Class label: 1
    [14.23 1.71 2.43 15.60 127.00 2.80 3.06 0.28 2.29 5.64 1.04 3.92 1065.00]

    Class label frequencies
    Class 1 samples: 33.15%
    Class 2 samples: 39.89%
    Class 3 samples: 26.97%

## Visualizating of a dataset

[[back to top]](#Sections)

There are endless way to visualize datasets for get an initial idea of how the data looks like. The most common ones are probably histograms and scatter plots.

### Histograms

[[back to top]](#Sections)

Histograms are a useful data to explore the distribution of each feature across the different classes. This could provide us with intuitive insights which features have a good and not-so-good inter-class separation. Below, we will plot a sample histogram for the "Alcohol content" feature for the three wine classes.

```python
%matplotlib inline
```

```python
from matplotlib import pyplot as plt
from math import floor, ceil # for rounding up and down

plt.figure(figsize=(10,8))

# bin width of the histogram in steps of 0.15
bins = np.arange(floor(min(X_wine[:,0])), ceil(max(X_wine[:,0])), 0.15)

# get the max count for a particular bin for all classes combined
max_bin = max(np.histogram(X_wine[:,0], bins=bins)[0])

# the order of the colors for each histogram
colors = ('blue', 'red', 'green')

for label,color in zip(
        range(1,4), colors):

    mean = np.mean(X_wine[:,0][y_wine == label]) # class sample mean
    stdev = np.std(X_wine[:,0][y_wine == label]) # class standard deviation
    plt.hist(X_wine[:,0][y_wine == label],
             bins=bins,
             alpha=0.3, # opacity level
             label='class {} ($\mu={:.2f}$, $\sigma={:.2f}$)'.format(label, mean, stdev),
             color=color)

plt.ylim([0, max_bin*1.3])
plt.title('Wine data set - Distribution of alocohol contents')
plt.xlabel('alcohol by volume', fontsize=14)
plt.ylabel('count', fontsize=14)
plt.legend(loc='upper right')

plt.show()
```

### Scatterplots

[[back to top]](#Sections)

Scatter plots are useful for visualizing features in more than just one dimension, for example to get a feeling for the correlation between particular features.
Unfortunately, we can't plot all 13 features here at once, since the visual cortex of us humans is limited to a maximum of three dimensions.

Below, we will create an example 2D-Scatter plot from the features "Alcohol content" and "Malic acid content".
Additionally, we will use the [`scipy.stats.pearsonr`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.pearsonr.html) function to calculate a Pearson correlation coefficient between these two features.

```python
from scipy.stats import pearsonr

plt.figure(figsize=(10,8))

for label,marker,color in zip(
        range(1,4),('x', 'o', '^'),('blue', 'red', 'green')):

    # Calculate Pearson correlation coefficient
    R = pearsonr(X_wine[:,0][y_wine == label], X_wine[:,1][y_wine == label])
    plt.scatter(x=X_wine[:,0][y_wine == label], # x-axis: feat. from col. 1
                y=X_wine[:,1][y_wine == label], # y-axis: feat. from col. 2
                marker=marker, # data point symbol for the scatter plot
                color=color,
                alpha=0.7,
                label='class {:}, R={:.2f}'.format(label, R[0]) # label for the legend
                )

plt.title('Wine Dataset')
plt.xlabel('alcohol by volume in percent')
plt.ylabel('malic acid in g/l')
plt.legend(loc='upper right')

plt.show()
```

If we want to pack 3 different features into one scatter plot at once, we can also do the same thing in 3D:

```python
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')

for label,marker,color in zip(
        range(1,4),('x', 'o', '^'),('blue','red','green')):

    ax.scatter(X_wine[:,0][y_wine == label],
               X_wine[:,1][y_wine == label],
               X_wine[:,2][y_wine == label],
               marker=marker,
               color=color,
               s=40,
               alpha=0.7,
               label='class {}'.format(label))

ax.set_xlabel('alcohol by volume in percent')
ax.set_ylabel('malic acid in g/l')
ax.set_zlabel('ash content in g/l')

plt.title('Wine dataset')

plt.show()
```

## Splitting into training and test dataset

[[back to top]](#Sections)

It is a typical procedure for machine learning and pattern classification tasks to split one dataset into two: a training dataset and a test dataset.
The training dataset is henceforth used to train our algorithms or classifier, and the test dataset is a way to validate the outcome quite objectively before we apply it to "new, real world data".

Here, we will split the dataset randomly so that 70% of the total dataset will become our training dataset, and 30% will become our test dataset, respectively.

```python
from sklearn.cross_validation import train_test_split
from sklearn import preprocessing

X_train, X_test, y_train, y_test = train_test_split(X_wine, y_wine,
     test_size=0.30, random_state=123)
```

    /srv/venv/lib/python3.6/site-packages/sklearn/cross_validation.py:41: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.
      "This module will be removed in 0.20.", DeprecationWarning)

Note that since this a random assignment, the original relative frequencies for each class label are not maintained.

```python
print('Class label frequencies')

print('\nTraining Dataset:')
for l in range(1,4):
    print('Class {:} samples: {:.2%}'.format(l, list(y_train).count(l)/y_train.shape[0]))

print('\nTest Dataset:')
for l in range(1,4):
    print('Class {:} samples: {:.2%}'.format(l, list(y_test).count(l)/y_test.shape[0]))
```

    Class label frequencies

    Training Dataset:
    Class 1 samples: 36.29%
    Class 2 samples: 42.74%
    Class 3 samples: 20.97%

    Test Dataset:
    Class 1 samples: 25.93%
    Class 2 samples: 33.33%
    Class 3 samples: 40.74%

## Feature Scaling

[[back to top]](#Sections)

### Standardization

Another important procedure is to standardize the data prior to fitting the model and other analyses so that the features will have the properties of a standard normal distribution with

$\mu = 0$ and $\sigma = 1$

where $\mu$ is the mean (average) and $\sigma$ is the standard deviation from the mean; standard scores (also called ***z*** scores) of the samples are calculated as follows:

\begin{equation} z = \frac{x - \mu}{\sigma}\end{equation}

Standardizing the features so that they are centered around 0 with a standard deviation of 1 is especially important if we are comparing measurements that have different units, e.g., in our "wine data" example, where the alcohol content is measured in volume percent, and the malic acid content in g/l.

```python
std_scale = preprocessing.StandardScaler().fit(X_train)
X_train = std_scale.transform(X_train)
X_test = std_scale.transform(X_test)
```

```python
f, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(10,5))

for a,x_dat, y_lab in zip(ax, (X_train, X_test), (y_train, y_test)):

    for label,marker,color in zip(
        range(1,4),('x', 'o', '^'),('blue','red','green')):

        a.scatter(x=x_dat[:,0][y_lab == label],
                y=x_dat[:,1][y_lab == label],
                marker=marker,
                color=color,
                alpha=0.7,
                label='class {}'.format(label)
                )

    a.legend(loc='upper left')

ax[0].set_title('Training Dataset')
ax[1].set_title('Test Dataset')
f.text(0.5, 0.04, 'malic acid (standardized)', ha='center', va='center')
f.text(0.08, 0.5, 'alcohol (standardized)', ha='center', va='center', rotation='vertical')

plt.show()
```

### Min-Max scaling (Normalization)

[[back to top]](#Sections)

An alternative approach to standardization is the so-called Min-Max scaling (sometimes also referred to as "normalization").
In this approach, the data is scaled to a fixed range - usually 0 to 1.
The cost of having this bounded range - in contrast to standardization - is that we will end up with small standard deviations, for example in the case where outliers are present.

The equation to calculate the "normalized" score is:

\begin{equation} X' = \frac{X - X_{min}}{X_{max}-X_{min}} \end{equation}

```python
minmax_scale = preprocessing.MinMaxScaler(feature_range=(0, 1)).fit(X_train)
X_train_minmax = minmax_scale.transform(X_train)
X_test_minmax = minmax_scale.transform(X_test)
```

```python
f, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(10,5))

for a,x_dat, y_lab in zip(ax, (X_train_minmax, X_test_minmax), (y_train, y_test)):

    for label,marker,color in zip(
        range(1,4),('x', 'o', '^'),('blue','red','green')):

        a.scatter(x=x_dat[:,0][y_lab == label],
                y=x_dat[:,1][y_lab == label],
                marker=marker,
                color=color,
                alpha=0.7,
                label='class {}'.format(label)
                )

    a.legend(loc='upper left')

ax[0].set_title('Training Dataset')
ax[1].set_title('Test Dataset')
f.text(0.5, 0.04, 'malic acid (normalized)', ha='center', va='center')
f.text(0.08, 0.5, 'alcohol (normalized)', ha='center', va='center', rotation='vertical')

plt.show()
```

## Linear Transformation: Principal Component Analysis (PCA)

[[back to top]](#Sections)

The main purposes of a principal component analysis are the analysis of data to identify patterns and finding patterns to reduce the dimensions of the dataset with minimal loss of information.

Here, our desired outcome of the principal component analysis is to project a feature space (our dataset consisting of n x d-dimensional samples) onto a smaller subspace that represents our data "well". A possible application would be a pattern classification task, where we want to reduce the computational costs and the error of parameter estimation by reducing the number of dimensions of our feature space by extracting a subspace that describes our data "best".

If you are interested in the Principal Component Analysis in more detail, I have outlined the procedure in a separate article
["Implementing a Principal Component Analysis (PCA) in Python step by step](http://sebastianraschka.com/Articles/2014_pca_step_by_step.html).

Here, we will use the [`sklearn.decomposition.PCA`](http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html) to transform our training data onto 2 dimensional subspace:

```python
from sklearn.decomposition import PCA
sklearn_pca = PCA(n_components=2) # number of components to keep
sklearn_transf = sklearn_pca.fit_transform(X_train)

plt.figure(figsize=(10,8))

for label,marker,color in zip(
        range(1,4),('x', 'o', '^'),('blue', 'red', 'green')):

    plt.scatter(x=sklearn_transf[:,0][y_train == label],
                y=sklearn_transf[:,1][y_train == label],
                marker=marker,
                color=color,
                alpha=0.7,
                label='class {}'.format(label)
                )

plt.xlabel('vector 1')
plt.ylabel('vector 2')

plt.legend()
plt.title('Most significant singular vectors after linear transformation via PCA')

plt.show()
```

#### PCA for feature extraction

As mentioned in the short introduction above (and in more detail in my separate [PCA article](http://sebastianraschka.com/Articles/2014_pca_step_by_step.html)), PCA is commonly used in the field of pattern classification for feature selection (or dimensionality reduction).
By default, the transformed data will be ordered by the components with the maximum variance (in descending order).

In the example above, I only kept the top 2 components (the 2 components with the maximum variance along the axes): The sample space of projected onto a 2-dimensional subspace, which was basically sufficient for plotting the data onto a 2D scatter plot.

However, if we want to use PCA for feature selection, we probably don't want to reduce the dimensionality that drastically. By default, the `PCA` function (`PCA(n_components=None)`) keeps all the components in ranked order. So we could basically either set the number `n_components` to a smaller size then the input dataset, or we could extract the top **n** components later from the returned NumPy array.

To get an idea about how well each components (relatively) "explains" the variance, we can use `explained_variance_ratio_` instant method, which also confirms that the components are ordered from most explanatory to least explanatory (the ratios sum up to 1.0).

```python
sklearn_pca = PCA(n_components=None)
sklearn_transf = sklearn_pca.fit_transform(X_train)
sklearn_pca.explained_variance_ratio_
```

    array([0.36, 0.21, 0.10, 0.08, 0.06, 0.05, 0.04, 0.03, 0.02, 0.02, 0.01,
           0.01, 0.01])

## Linear Transformation: Linear Discriminant Analysis (MDA)

[[back to top]](#Sections)

The main purposes of a Linear Discriminant Analysis (LDA) is to analyze the data to identify patterns to project it onto a subspace that yields a better separation of the classes. Also, the dimensionality of the dataset shall be reduced with minimal loss of information.

**The approach is very similar to a Principal Component Analysis (PCA), but in addition to finding the component axes that maximize the variance of our data, we are additionally interested in the axes that maximize the separation of our classes (e.g., in a supervised pattern classification problem)**

Here, our desired outcome of the Linear discriminant analysis is to project a feature space (our dataset consisting of n d-dimensional samples) onto a smaller subspace that represents our data "well" and has a good class separation. A possible application would be a pattern classification task, where we want to reduce the computational costs and the error of parameter estimation by reducing the number of dimensions of our feature space by extracting a subspace that describes our data "best".

#### Principal Component Analysis (PCA) Vs. Linear Discriminant Analysis (LDA)

Both Linear Discriminant Analysis and Principal Component Analysis are linear transformation methods and closely related to each other. In PCA, we are interested to find the directions (components) that maximize the variance in our dataset, where in LDA, we are additionally interested to find the directions that maximize the separation (or discrimination) between different classes (for example, in pattern classification problems where our dataset consists of multiple classes. In contrast two PCA, which ignores the class labels).

**In other words, via PCA, we are projecting the entire set of data (without class labels) onto a different subspace, and in LDA, we are trying to determine a suitable subspace to distinguish between patterns that belong to different classes. Or, roughly speaking in PCA we are trying to find the axes with maximum variances where the data is most spread (within a class, since PCA treats the whole data set as one class), and in LDA we are additionally maximizing the spread between classes.**

In typical pattern recognition problems, a PCA is often followed by an LDA.

If you are interested, you can find more information about the LDA in my IPython notebook
[Stepping through a Linear Discriminant Analysis - using Python's NumPy and matplotlib](http://nbviewer.ipython.org/github/rasbt/pattern_classification/blob/master/dimensionality_reduction/projection/minear_discriminant_analysis.ipynb?create=1).

Like we did in the PCA section above, we will use a `scikit-learn` funcion, [`sklearn.lda.LDA`](http://scikit-learn.org/stable/modules/generated/sklearn.lda.LDA.html) in order to transform our training data onto 2 dimensional subspace:

```python
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
sklearn_lda = LDA(n_components=2)
transf_lda = sklearn_lda.fit_transform(X_train, y_train)

plt.figure(figsize=(10,8))

for label,marker,color in zip(
        range(1,4),('x', 'o', '^'),('blue', 'red', 'green')):

    plt.scatter(x=transf_lda[:,0][y_train == label],
                y=transf_lda[:,1][y_train == label],
                marker=marker,
                color=color,
                alpha=0.7,
                label='class {}'.format(label)
                )

plt.xlabel('vector 1')
plt.ylabel('vector 2')

plt.legend()
plt.title('Most significant singular vectors after linear transformation via LDA')

plt.show()
```

#### LDA for feature extraction

If we want to use LDA for projecting our data onto a smaller subspace (i.e., for dimensionality reduction), we can directly set the number of components to keep via `LDA(n_components=...)`; this is analogous to the [PCA function](#PCA-for-feature-extraction), which we have seen above.

## Simple Supervised Classification

[[back to top]](#Sections)

### Linear Discriminant Analysis as simple linear classifier

[[back to top]](#Sections)

The LDA that we've just used in the section above can also be used as a simple linear classifier.

```python
# fit model
lda_clf = LDA()
lda_clf.fit(X_train, y_train)
LDA(n_components=None, priors=None)

# prediction
print('1st sample from test dataset classified as:', lda_clf.predict(X_test[0:]))
print('actual class label:', y_test[0])
```

    1st sample from test dataset classified as: [3 2 3 2 2 3 1 3 3 2 3 3 3 1 1 3 2 2 1 2 3 3 3 3 2 3 3 2 1 1 1 1 3 2 3 2 3
     1 2 2 3 3 1 1 1 1 1 2 1 2 2 3 3 2]
    actual class label: 3

Another handy subpackage of sklearn is `metrics`. The [`metrics.accuracy_score`](http://scikit-learn.org/stable/modules/generated/sklearn.metrics.accuracy_score.html), for example, is quite useful to evaluate how many samples can be classified correctly:

```python
from sklearn import metrics
pred_train_lda = lda_clf.predict(X_train)

print('Prediction accuracy for the training dataset')
print('{:.2%}'.format(metrics.accuracy_score(y_train, pred_train_lda)))
```

    Prediction accuracy for the training dataset
    100.00%

To verify that over model was not overfitted to the training dataset, let us evaluate the classifier's accuracy on the test dataset:

```python
pred_test_lda = lda_clf.predict(X_test)

print('Prediction accuracy for the test dataset')
print('{:.2%}'.format(metrics.accuracy_score(y_test, pred_test_lda)))
```

    Prediction accuracy for the test dataset
    98.15%

**Confusion Matrix**
As we can see above, there was a very low misclassification rate when we'd apply the classifier on the test data set. A confusion matrix can tell us in more detail which particular classes could not classified correctly.

	<colgroup width="60"></colgroup>
	<colgroup span="4" width="82"></colgroup>

		predicted class

		class 1
		class 2
		class 3

		actual class
		class 1
		True positives

		class 2

		True positives

		class 3

		True positives

```python
print('Confusion Matrix of the LDA-classifier')
print(metrics.confusion_matrix(y_test, lda_clf.predict(X_test)))
```

    Confusion Matrix of the LDA-classifier
    [[14  0  0]
     [ 1 17  0]
     [ 0  0 22]]

As we can see, one sample from class 2 was incorrectly labeled as class 1, from the perspective of class 1, this would be 1 "False Negative" or a "False Postive" from the perspective of class 2, respectively

### Classification Stochastic Gradient Descent (SGD)

[[back to top]](#Sections)

Let us now compare the classification accuracy of the LDA classifier with a simple classification (we also use the probably not ideal default settings here) via stochastic gradient descent, an algorithm that minimizes a linear objective function.
More information about the `sklearn.linear_model.SGDClassifier` can be found [here](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.SGDClassifier.html).

```python
from sklearn.linear_model import SGDClassifier

sgd_clf = SGDClassifier()
sgd_clf.fit(X_train, y_train)

pred_train_sgd = sgd_clf.predict(X_train)
pred_test_sgd = sgd_clf.predict(X_test)

print('\nPrediction accuracy for the training dataset')
print('{:.2%}\n'.format(metrics.accuracy_score(y_train, pred_train_sgd)))

print('Prediction accuracy for the test dataset')
print('{:.2%}\n'.format(metrics.accuracy_score(y_test, pred_test_sgd)))

print('Confusion Matrix of the SGD-classifier')
print(metrics.confusion_matrix(y_test, sgd_clf.predict(X_test)))
```

    Prediction accuracy for the training dataset
    100.00%

    Prediction accuracy for the test dataset
    94.44%

    Confusion Matrix of the SGD-classifier
    [[14  0  0]
     [ 1 17  0]
     [ 0  2 20]]

    /srv/venv/lib/python3.6/site-packages/sklearn/linear_model/stochastic_gradient.py:84: FutureWarning: max_iter and tol parameters have been added in <class 'sklearn.linear_model.stochastic_gradient.SGDClassifier'> in 0.19. If both are left unset, they default to max_iter=5 and tol=None. If tol is not None, max_iter defaults to max_iter=1000. From 0.21, default max_iter will be 1000, and default tol will be 1e-3.
      "and default tol will be 1e-3." % type(self), FutureWarning)

Quite impressively, we achieved a 100% prediction accuracy on the test dataset without any additional efforts of tweaking any parameters and settings.

### Decision Regions

```python
sgd_clf2 = SGDClassifier()
sgd_clf2.fit(X_train[:, :2], y_train)

x_min = X_test[:, 0].min()
x_max = X_test[:, 0].max()
y_min = X_test[:, 1].min()
y_max = X_test[:, 1].max()

step = 0.01
X, Y = np.meshgrid(np.arange(x_min, x_max, step), np.arange(y_min, y_max, step))

Z = sgd_clf2.predict(np.c_[X.ravel(), Y.ravel()])
Z = Z.reshape(X.shape)

# Plots decision regions
plt.contourf(X, Y, Z)

# Plots samples from training data set
plt.scatter(X_train[:, 0], X_train[:, 1], c=y_train)
plt.show()
```

## Saving the processed datasets

[[back to top]](#Sections)

### Pickle

[[back to top]](#Sections)

The in-built [`pickle`](https://docs.python.org/3.4/library/pickle.html) module is a convenient tool in Python's standard library to save Python objects in byte format. This allows us, for example, to save our NumPy arrays and classifiers so that we can load them in a later or different Python session to continue working with our data, e.g., to train a classifier.

```python
# export objects via pickle

import pickle

pickle_out = open('standardized_data.pkl', 'wb')
pickle.dump([X_train, X_test, y_train, y_test], pickle_out)
pickle_out.close()

pickle_out = open('classifiers.pkl', 'wb')
pickle.dump([lda_clf, sgd_clf], pickle_out)
pickle_out.close()
```

```python
# import objects via pickle

my_object_file = open('standardized_data.pkl', 'rb')
X_train, X_test, y_train, y_test = pickle.load(my_object_file)
my_object_file.close()

my_object_file = open('classifiers.pkl', 'rb')
lda_clf, sgd_clf = pickle.load(my_object_file)
my_object_file.close()

print('Confusion Matrix of the SGD-classifier')
print(metrics.confusion_matrix(y_test, sgd_clf.predict(X_test)))
```

    Confusion Matrix of the SGD-classifier
    [[14  0  0]
     [ 1 17  0]
     [ 0  2 20]]

### Comma-Separated-Values

[[back to top]](#Sections)

And it is also always a good idea to save your data in common text formats, such as the CSV format that we started with. But first, let us add back the class labels to the front column of the test and training data sets.

```python
training_data = np.hstack((y_train.reshape(y_train.shape[0], 1), X_train))
test_data = np.hstack((y_test.reshape(y_test.shape[0], 1), X_test))
```

Now, we can save our test and training datasets as 2 separate CSV files using the [`numpy.savetxt`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.savetxt.html) function.

```python
np.savetxt('./training_set.csv', training_data, delimiter=',')
np.savetxt('./test_set.csv', test_data, delimiter=',')
```

```python

```

```python

```

max_stars_repo_name: gopala-kr/ds-notebooks

max_issues_repo_name: gopala-kr/ds-notebooks

max_forks_repo_name: gopala-kr/ds-notebooks


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

Autograd: 자동 미분
===================================

PyTorch의 모든 신경망의 중심에는 ``autograd`` 패키지가 있습니다.
먼저 이것을 가볍게 살펴본 뒤, 첫번째 신경망을 학습시켜보겠습니다.

``autograd`` 패키지는 Tensor의 모든 연산에 대해 자동 미분을 제공합니다.
이는 실행-기반-정의(define-by-run) 프레임워크로, 이는 코드를 어떻게 작성하여
실행하느냐에 따라 역전파가 정의된다는 뜻이며, 역전파는 학습 과정의 매 단계마다
달라집니다.

더 간단한 용어로 몇 가지 예를 살펴보겠습니다.

Tensor
--------

패키지의 중심에는 ``torch.Tensor`` 클래스가 있습니다. 만약 ``.requires_grad``
속성을 ``True`` 로 설정하면, 그 tensor에서 이뤄진 모든 연산들을 추적(track)하기
시작합니다. 계산이 완료된 후 ``.backward()`` 를 호출하여 모든 변화도(gradient)를
자동으로 계산할 수 있습니다. 이 Tensor의 변화도는 ``.grad`` 속성에 누적됩니다.

Tensor가 기록을 추적하는 것을 중단하게 하려면, ``.detach()`` 를 호출하여 연산
기록으로부터 분리(detach)하여 이후 연산들이 추적되는 것을 방지할 수 있습니다.

기록을 추적하는 것(과 메모리를 사용하는 것)을 방지하기 위해, 코드 블럭을
``with torch.no_grad():`` 로 감쌀 수 있습니다. 이는 특히 변화도(gradient)는
필요없지만, `requires_grad=True` 가 설정되어 학습 가능한 매개변수를 갖는 모델을
평가(evaluate)할 때 유용합니다.

Autograd 구현에서 매우 중요한 클래스가 하나 더 있는데, 이것은 바로 ``Function``
클래스입니다.

``Tensor`` 와 ``Function`` 은 서로 연결되어 있으며, 모든 연산 과정을
부호화(encode)하여 순환하지 않는 그래프(acyclic graph)를 생성합니다. 각 tensor는
``.grad_fn`` 속성을 갖고 있는데, 이는 ``Tensor`` 를 생성한 ``Function`` 을
참조하고 있습니다. (단, 사용자가 만든 Tensor는 예외로, 이 때 ``grad_fn`` 은
``None`` 입니다.)

도함수를 계산하기 위해서는 ``Tensor`` 의 ``.backward()`` 를 호출하면
됩니다. 만약 ``Tensor`` 가 스칼라(scalar)인 경우(예. 하나의 요소 값만 갖는 등)에는
``backward`` 에 인자를 정해줄 필요가 없습니다. 하지만 여러 개의 요소를 갖고 있을
때는 tensor의 모양을 ``gradient`` 의 인자로 지정할 필요가 있습니다.

```python
import torch
```

tensor를 생성하고 ``requires_grad=True`` 를 설정하여 연산을 기록합니다.

```python
x = torch.ones(2, 2, requires_grad=True)
print(x)
```

tensor에 연산을 수행합니다:

```python
y = x + 2
print(y)
```

``y`` 는 연산의 결과로 생성된 것이므로 ``grad_fn`` 을 갖습니다.

```python
print(y.grad_fn)
```

``y`` 에 다른 연산을 수행합니다.

```python
z = y * y * 3
out = z.mean()

print(z, out)
```

``.requires_grad_( ... )`` 는 기존 Tensor의 ``requires_grad`` 값을 바꿔치기
(in-place)하여 변경합니다. 입력값이 지정되지 않으면 기본값은 ``False`` 입니다.

```python
a = torch.randn(2, 2)
a = ((a * 3) / (a - 1))
print(a.requires_grad)
a.requires_grad_(True)
print(a.requires_grad)
b = (a * a).sum()
print(b.grad_fn)
```

변화도(Gradient)
-----------------
이제 역전파(backprop)를 해보겠습니다.
``out`` 은 하나의 스칼라 값만 갖고 있기 때문에, ``out.backward()`` 는
``out.backward(torch.tensor(1.))`` 과 동일합니다.

```python
out.backward()
```

변화도 d(out)/dx를 출력합니다.

```python
print(x.grad)
```

``4.5`` 로 이루어진 행렬을 확인할 수 있습니다. ``out`` 을 *Tensor* “$o$”
라고 하면, 다음과 같이 구할 수 있습니다.
$o = \frac{1}{4}\sum_i z_i$ 이고,
$z_i = 3(x_i+2)^2$ 이므로 $z_i\bigr\rvert_{x_i=1} = 27$ 입니다.
따라서,
$\frac{\partial o}{\partial x_i} = \frac{3}{2}(x_i+2)$ 이므로,
$\frac{\partial o}{\partial x_i}\bigr\rvert_{x_i=1} = \frac{9}{2} = 4.5$ 입니다.

수학적으로 벡터 함수 $\vec{y}=f(\vec{x})$ 에서 $\vec{x}$ 에
대한 $\vec{y}$ 의 변화도는 야코비안 행렬(Jacobian Matrix)입니다:

\begin{align}J=\left(\begin{array}{ccc}
   \frac{\partial y_{1}}{\partial x_{1}} & \cdots & \frac{\partial y_{1}}{\partial x_{n}}\\
   \vdots & \ddots & \vdots\\
   \frac{\partial y_{m}}{\partial x_{1}} & \cdots & \frac{\partial y_{m}}{\partial x_{n}}
   \end{array}\right)\end{align}

일반적으로, ``torch.autograd`` 는 벡터-야코비안 곱을 계산하는 엔진입니다. 즉,
어떤 벡터 $v=\left(\begin{array}{cccc} v_{1} & v_{2} & \cdots & v_{m}\end{array}\right)^{T}$
에 대해 $v^{T}\cdot J$ 을 연산합니다. 만약 $v$ 가 스칼라 함수
$l=g\left(\vec{y}\right)$ 의 기울기인 경우,
$v=\left(\begin{array}{ccc}\frac{\partial l}{\partial y_{1}} & \cdots & \frac{\partial l}{\partial y_{m}}\end{array}\right)^{T}$
이며, 연쇄법칙(chain rule)에 따라 벡터-야코비안 곱은 $\vec{x}$ 에 대한
$l$ 의 기울기가 됩니다:

\begin{align}J^{T}\cdot v=\left(\begin{array}{ccc}
   \frac{\partial y_{1}}{\partial x_{1}} & \cdots & \frac{\partial y_{m}}{\partial x_{1}}\\
   \vdots & \ddots & \vdots\\
   \frac{\partial y_{1}}{\partial x_{n}} & \cdots & \frac{\partial y_{m}}{\partial x_{n}}
   \end{array}\right)\left(\begin{array}{c}
   \frac{\partial l}{\partial y_{1}}\\
   \vdots\\
   \frac{\partial l}{\partial y_{m}}
   \end{array}\right)=\left(\begin{array}{c}
   \frac{\partial l}{\partial x_{1}}\\
   \vdots\\
   \frac{\partial l}{\partial x_{n}}
   \end{array}\right)\end{align}

(여기서 $v^{T}\cdot J$ 은 $J^{T}\cdot v$ 를 취했을 때의 열 벡터로
취급할 수 있는 행 벡터를 갖습니다.)

벡터-야코비안 곱의 이러한 특성은 스칼라가 아닌 출력을 갖는 모델에 외부 변화도를
제공(feed)하는 것을 매우 편리하게 해줍니다.

이제 벡터-야코비안 곱의 예제를 살펴보도록 하겠습니다:

```python
x = torch.randn(3, requires_grad=True)

y = x * 2
while y.data.norm() < 1000:
    y = y * 2

print(y)
```

이 경우 ``y`` 는 더 이상 스칼라 값이 아닙니다. ``torch.autograd`` 는
전체 야코비안을 직접 계산할수는 없지만, 벡터-야코비안 곱은 간단히
``backward`` 에 해당 벡터를 인자로 제공하여 얻을 수 있습니다:

```python
v = torch.tensor([0.1, 1.0, 0.0001], dtype=torch.float)
y.backward(v)

print(x.grad)
```

또한 ``with torch.no_grad():`` 로 코드 블럭을 감싸서 autograd가
``.requires_grad=True`` 인 Tensor들의 연산 기록을 추적하는 것을 멈출 수 있습니다.

```python
print(x.requires_grad)
print((x ** 2).requires_grad)

with torch.no_grad():
	print((x ** 2).requires_grad)
```

또는 ``.detach()`` 를 호출하여 내용물(content)은 같지만 require_grad가 다른
새로운 Tensor를 가져옵니다:

```python
print(x.requires_grad)
y = x.detach()
print(y.requires_grad)
print(x.eq(y).all())
```

**더 읽을거리:**

``autograd.Function`` 관련 문서는 https://pytorch.org/docs/stable/autograd.html#function
에서 찾아볼 수 있습니다.

ext: ipynb

max_stars_repo_name: pleiades-s/PyTorch-tutorials-kr

max_issues_repo_name: pleiades-s/PyTorch-tutorials-kr

max_forks_repo_name: pleiades-s/PyTorch-tutorials-kr


text:
```python
# Installation of PM4Py Python package
#Source: https://pm4py.fit.fraunhofer.de/install-page
!pip install pm4py
```

    Collecting pm4py
      Downloading https://files.pythonhosted.org/packages/d2/83/dcfa7bfa85d7fda35417f4bf76c8e5980919de92aa8d458c28d1518cc162/pm4py-2.2.10.1-py3-none-any.whl (1.4MB)
         |████████████████████████████████| 1.5MB 5.1MB/s
    Requirement already satisfied: pytz in /usr/local/lib/python3.7/dist-packages (from pm4py) (2018.9)
    Collecting pulp<=2.1
      Downloading https://files.pythonhosted.org/packages/41/34/757c88c320f80ce602199603afe63aed1e0bc11180b9a9fb6018fb2ce7ef/PuLP-2.1-py3-none-any.whl (40.6MB)
         |████████████████████████████████| 40.6MB 73kB/s
    Requirement already satisfied: pydotplus in /usr/local/lib/python3.7/dist-packages (from pm4py) (2.0.2)
    Requirement already satisfied: matplotlib in /usr/local/lib/python3.7/dist-packages (from pm4py) (3.2.2)
    Requirement already satisfied: intervaltree in /usr/local/lib/python3.7/dist-packages (from pm4py) (2.1.0)
    Requirement already satisfied: cvxopt in /usr/local/lib/python3.7/dist-packages (from pm4py) (1.2.6)
    Collecting pyvis
      Downloading https://files.pythonhosted.org/packages/07/d1/e87844ec86e96df7364f21af2263ad6030c0d727660ae89935c7af56a540/pyvis-0.1.9-py3-none-any.whl
    Requirement already satisfied: graphviz in /usr/local/lib/python3.7/dist-packages (from pm4py) (0.10.1)
    Requirement already satisfied: pandas in /usr/local/lib/python3.7/dist-packages (from pm4py) (1.1.5)
    Collecting deprecation
      Downloading https://files.pythonhosted.org/packages/02/c3/253a89ee03fc9b9682f1541728eb66db7db22148cd94f89ab22528cd1e1b/deprecation-2.1.0-py2.py3-none-any.whl
    Requirement already satisfied: scikit-learn in /usr/local/lib/python3.7/dist-packages (from pm4py) (0.22.2.post1)
    Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (from pm4py) (1.4.1)
    Requirement already satisfied: lxml in /usr/local/lib/python3.7/dist-packages (from pm4py) (4.2.6)
    Requirement already satisfied: sympy in /usr/local/lib/python3.7/dist-packages (from pm4py) (1.7.1)
    Collecting jsonpickle
      Downloading https://files.pythonhosted.org/packages/bb/1a/f2db026d4d682303793559f1c2bb425ba3ec0d6fd7ac63397790443f2461/jsonpickle-2.0.0-py2.py3-none-any.whl
    Requirement already satisfied: networkx in /usr/local/lib/python3.7/dist-packages (from pm4py) (2.5.1)
    Requirement already satisfied: numpy!=1.19.4 in /usr/local/lib/python3.7/dist-packages (from pm4py) (1.19.5)
    Requirement already satisfied: tqdm in /usr/local/lib/python3.7/dist-packages (from pm4py) (4.41.1)
    Collecting stringdist
      Downloading https://files.pythonhosted.org/packages/85/f0/c56cbe92b4b06fbc7adaa81917ad34d7027834e166fff2d2db73961c67fa/StringDist-1.0.9.tar.gz
    Requirement already satisfied: pyparsing>=2.0.1 in /usr/local/lib/python3.7/dist-packages (from pulp<=2.1->pm4py) (2.4.7)
    Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.7/dist-packages (from matplotlib->pm4py) (0.10.0)
    Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->pm4py) (1.3.1)
    Requirement already satisfied: python-dateutil>=2.1 in /usr/local/lib/python3.7/dist-packages (from matplotlib->pm4py) (2.8.1)
    Requirement already satisfied: sortedcontainers in /usr/local/lib/python3.7/dist-packages (from intervaltree->pm4py) (2.4.0)
    Requirement already satisfied: jinja2>=2.9.6 in /usr/local/lib/python3.7/dist-packages (from pyvis->pm4py) (2.11.3)
    Requirement already satisfied: ipython>=5.3.0 in /usr/local/lib/python3.7/dist-packages (from pyvis->pm4py) (5.5.0)
    Requirement already satisfied: packaging in /usr/local/lib/python3.7/dist-packages (from deprecation->pm4py) (21.0)
    Requirement already satisfied: joblib>=0.11 in /usr/local/lib/python3.7/dist-packages (from scikit-learn->pm4py) (1.0.1)
    Requirement already satisfied: mpmath>=0.19 in /usr/local/lib/python3.7/dist-packages (from sympy->pm4py) (1.2.1)
    Requirement already satisfied: importlib-metadata; python_version < "3.8" in /usr/local/lib/python3.7/dist-packages (from jsonpickle->pm4py) (4.6.1)
    Requirement already satisfied: decorator<5,>=4.3 in /usr/local/lib/python3.7/dist-packages (from networkx->pm4py) (4.4.2)
    Requirement already satisfied: six in /usr/local/lib/python3.7/dist-packages (from cycler>=0.10->matplotlib->pm4py) (1.15.0)
    Requirement already satisfied: MarkupSafe>=0.23 in /usr/local/lib/python3.7/dist-packages (from jinja2>=2.9.6->pyvis->pm4py) (2.0.1)
    Requirement already satisfied: pexpect; sys_platform != "win32" in /usr/local/lib/python3.7/dist-packages (from ipython>=5.3.0->pyvis->pm4py) (4.8.0)
    Requirement already satisfied: simplegeneric>0.8 in /usr/local/lib/python3.7/dist-packages (from ipython>=5.3.0->pyvis->pm4py) (0.8.1)
    Requirement already satisfied: pygments in /usr/local/lib/python3.7/dist-packages (from ipython>=5.3.0->pyvis->pm4py) (2.6.1)
    Requirement already satisfied: setuptools>=18.5 in /usr/local/lib/python3.7/dist-packages (from ipython>=5.3.0->pyvis->pm4py) (57.2.0)
    Requirement already satisfied: traitlets>=4.2 in /usr/local/lib/python3.7/dist-packages (from ipython>=5.3.0->pyvis->pm4py) (5.0.5)
    Requirement already satisfied: pickleshare in /usr/local/lib/python3.7/dist-packages (from ipython>=5.3.0->pyvis->pm4py) (0.7.5)
    Requirement already satisfied: prompt-toolkit<2.0.0,>=1.0.4 in /usr/local/lib/python3.7/dist-packages (from ipython>=5.3.0->pyvis->pm4py) (1.0.18)
    Requirement already satisfied: typing-extensions>=3.6.4; python_version < "3.8" in /usr/local/lib/python3.7/dist-packages (from importlib-metadata; python_version < "3.8"->jsonpickle->pm4py) (3.7.4.3)
    Requirement already satisfied: zipp>=0.5 in /usr/local/lib/python3.7/dist-packages (from importlib-metadata; python_version < "3.8"->jsonpickle->pm4py) (3.5.0)
    Requirement already satisfied: ptyprocess>=0.5 in /usr/local/lib/python3.7/dist-packages (from pexpect; sys_platform != "win32"->ipython>=5.3.0->pyvis->pm4py) (0.7.0)
    Requirement already satisfied: ipython-genutils in /usr/local/lib/python3.7/dist-packages (from traitlets>=4.2->ipython>=5.3.0->pyvis->pm4py) (0.2.0)
    Requirement already satisfied: wcwidth in /usr/local/lib/python3.7/dist-packages (from prompt-toolkit<2.0.0,>=1.0.4->ipython>=5.3.0->pyvis->pm4py) (0.2.5)
    Building wheels for collected packages: stringdist
      Building wheel for stringdist (setup.py) ... done
      Created wheel for stringdist: filename=StringDist-1.0.9-cp37-cp37m-linux_x86_64.whl size=23562 sha256=1fe20d631016b567e747cf0f7a152dc67bfcfdff3279bc22ac1c98ba0d2af1a4
      Stored in directory: /root/.cache/pip/wheels/c3/1b/bb/bf0de4d64d8ca38759811fe3353e441f12feb606fb3d1d2e11
    Successfully built stringdist
    Installing collected packages: pulp, jsonpickle, pyvis, deprecation, stringdist, pm4py
    Successfully installed deprecation-2.1.0 jsonpickle-2.0.0 pm4py-2.2.10.1 pulp-2.1 pyvis-0.1.9 stringdist-1.0.9

```python
from google.colab import drive
drive.mount('/content/drive')
```

```python
# Import dependencies and set the filepath
import pm4py
import os
import pandas as pd
import numpy as np
import datetime as dt
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px         # needed for plotly express in general
import plotly.graph_objects as go   # needed for plotly graph objects e.g. pie
```

```python
# Set the file path
# Please upload the local file sent via email to colab before running this
file_path = '/content/Hospital Billing - Event Log.xes' # if using the original file
#file_path = '/content/log.xes' # if renamed
```

```python
# Read the XES file and create the log
event_log = pm4py.read_xes(file_path)
```

    HBox(children=(FloatProgress(value=0.0, description='parsing log, completed traces :: ', max=85269.0, style=Pr…

```python
def import_xes(file_path):
  event_log = pm4py.read_xes(file_path)
  start_activities = pm4py.get_start_activities(event_log)
  end_activities = pm4py.get_end_activities(event_log)
print("Start activities: {}\nEnd activities: {}".format(start_activities, end_activities))
```

```python
# Obtaining a Process Model
# event_log = pm4py.read_xes('C:/Users/demo/Downloads/running-example.xes')

process_tree = pm4py.discover_process_tree_inductive(event_log) # discover_tree_inductive is deprecated
bpmn_model = pm4py.convert_to_bpmn(process_tree)
pm4py.view_bpmn(bpmn_model)
```

```python
#obtain process map
dfg, start_activities, end_activities = pm4py.discover_dfg(event_log)
pm4py.view_dfg(dfg, start_activities, end_activities)
```

```python
#heuristics map
map = pm4py.discover_heuristics_net(event_log)
pm4py.view_heuristics_net(map)
```

```python
# View Process Tree
process_tree = pm4py.discover_process_tree_inductive(event_log) # discover_tree_inductive is deprecated
pm4py.view_process_tree(process_tree)
```

```python
print("Start activities: {}\nEnd activities: {}".format(start_activities, end_activities))
```

## Miner Algorithms
Please refer to the documentation here:
https://pm4py.fit.fraunhofer.de/documentation#item-3-1

```python
# Alpha Miner implementation

from pm4py.algo.discovery.alpha import algorithm as alpha_miner
net, initial_marking, final_marking = alpha_miner.apply(event_log)

# Inductive Miner implementation
from pm4py.algo.discovery.inductive import algorithm as inductive_miner
net, initial_marking, final_marking = inductive_miner.apply(event_log)
```

## Visualisations

```python
# Visualisation of IM

from pm4py.visualization.process_tree import visualizer as pt_visualizer
tree = inductive_miner.apply_tree(event_log)

gviz = pt_visualizer.apply(tree)
pt_visualizer.view(gviz)
```

    '/tmp/tmpn43cpox3.gv.png'

```python

```

```python
# Converting process tree in petri net

from pm4py.objects.conversion.process_tree import converter as pt_converter
net, initial_marking, final_marking = pt_converter.apply(tree, variant=pt_converter.Variants.TO_PETRI_NET)
gviz = pt_visualizer.apply(tree)
pt_visualizer.view(gviz)
```

    '/tmp/tmpgzmi0vfv.gv.png'

```python
# Dotted chart
pm4py.view_dotted_chart(event_log, format="svg")
```

Decision tree & Root Cause Analysis

```python
from pm4py.objects.log.util import get_class_representation , get_log_representation
from sklearn import tree
data, feature_names = get_log_representation.get_default_representation(event_log)

#changing into pandas data frame
dataframe = pd.DataFrame(data, columns=feature_names)

#creating decision tree
target, classes = get_class_representation.get_class_representation_by_trace_duration(event_log, 11007360)
clf = tree.DecisionTreeClassifier()
clf.fit(data, target)
```

    /usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:3: DeprecatedWarning:

    get_default_representation is deprecated as of 2.2.8 and will be removed in 3.0.0. please use pm4py.algo.transformation.log_to_features instead

    DecisionTreeClassifier(ccp_alpha=0.0, class_weight=None, criterion='gini',
                           max_depth=None, max_features=None, max_leaf_nodes=None,
                           min_impurity_decrease=0.0, min_impurity_split=None,
                           min_samples_leaf=1, min_samples_split=2,
                           min_weight_fraction_leaf=0.0, presort='deprecated',
                           random_state=None, splitter='best')

```python
#visualizing decision tree

from pm4py.visualization.decisiontree import visualizer as dectree_visualizer
from pm4py.objects.log.util import get_class_representation
gviz = dectree_visualizer.apply(clf, feature_names, classes)
target, classes = get_class_representation.get_class_representation_by_trace_duration(event_log, 11007360)
dectree_visualizer.view(gviz)
```

    '/tmp/tmpdu8qp6ui.gv.png'


text:
# Exercise 13.1 JT: Improved governance

The setup in this exercise:

* Two period model, $t=0,1$.

* Continuum of entrepreneurs with initial wealth $(A)$ distributed with the CDF $G(A)$ with support on $[0,I]$.

* At time $t=0$: Entrepreneurs can invest $I$ and borrow $(I-A)$ or invest in other firms.

* Entrepreneurs are risk neutral and thus have utility $u_e=c_0+c_1$.

* As before: Receive $R$ is investment is successful, with probability $p_H$ if entrepreneur behaves or $p_L$ if he shirks (misbehaves). Shirking yields private benefits of $B$.

* Now add consumers: Behavior results in an increasing savings function $S(r)$, with $S(r)=0$ for $r<0$.

## Q1: When are projects financed?

Consider the entrepreneur again (similar to the exercises above):

* He behaves if (IC) holds:

$$\begin{align}
    (\Delta p)R_b\geq B. \tag{IC}
\end{align}$$

* If (IC) holds profits are given by:

$$\begin{align}
    p_H\left(R-\dfrac{B}{\Delta p}\right). \tag{IC, $\pi$}
\end{align}$$

* Lenders will finance project as long as profits are larger than alternative cost:

$$\begin{align}
    (1+r)(I-A)\leq p_H\left(R-\dfrac{B}{\Delta p}\right), \tag{IR}
\end{align}$$

where the LHS is the profit from saving (alternative cost w. rent r) and the RHS is the profit from investing with a loan contract such that the (IC) constraint holds with equality.

* As we saw earlier the solution depends on the level $A$. All $A>A^*$ such that (IR) holds are financed.

## Q2: Give the expression determining the interest rate r. How does it change if investor protection improves (interpreted as fall in B)?

The **supply** side of the financial market consists of consumers with savings function $S(r)$ and entrepreneurs who do not get extended credit to carry out the project, i.e. all $A<A^*$. Thus:

$$\begin{align}
    \tilde{S}(r) = S(r)+\int_0^{A^*(r)}A\mbox{ }dG(A)
\end{align}$$

The **demand** side of the financial market consists of entrepreneurs with $A\geq A^*$:

$$\begin{align}
    D(r) = \int_{A^*(r)}^I (I-A)\mbox{ }dG(A)
\end{align}$$

Equilibrium balances demand and supply, i.e. $D(r)=\tilde{S}(r)$.

**Investor protection improves:** ($B$ falls)

* Motive for not behaving drops for all entrepreneurs. Thus entrepreneurs that were before marginally unprofitable to invest in, now increase the *pledgeable income* they can promise investors. In other words; the IC constraint is easier to satisfy, $R_b$ is lower.

* This lowers the threshold $A^*$, i.e. more entrepreneurs are extended credit. Demand increases.

* For supply to increase and adjust to the higher demand, the interest rate increaes.

**An extra comment on the model, The balance-sheet channel:**
Note that in this model, a marginal change in r has the effect on supply:

$$\begin{align}
    S'(r)+A^*(r)g(A^*(r))\dfrac{\partial A^*(r)}{\partial r}>0.
\end{align}$$
If $r$ increaess, less investment projects are profitable $\Rightarrow$ supply increases. The effet on demand is given by

$$\begin{align}
    -(I-A^*(r))dg(A^*(r))\leq 0
\end{align}$$

The implication is that *shocks that shift demand upwards (downwards) and supply downwards (upwards) can lead to large fluctuations in investment, without significant changes in the interest rate (r)*. Large fluctuations in investment, but small changes in r, is a stylized fact; something that we actually observe.

max_stars_repo_name: ChampionApe/FinancialFriction

max_issues_repo_name: ChampionApe/FinancialFriction

max_forks_repo_name: ChampionApe/FinancialFriction


text:
# Regularisation

In this lab we will explore different methods for regularising networks to reduce overfitting and improve generalisation. This uses the material covered in the [fifth lecture slides](http://www.inf.ed.ac.uk/teaching/courses/mlp/2016/mlp05-hid.pdf).

# Exercise 1: L1 and L2 penalties

In the previous lab notebook we explored the issue of overfitting. There we saw that this arises when the model is 'too complex' ($\sim$ has too many degrees of freedom / parameters) for the amount of data we have available.

One method for trying to reduce overfitting is therefore to try to decrease the flexibility of the model. We can do this by simply reducing the number of free parameters in the model (e.g. by using a shallower model with fewer layers or layers with smaller dimensionality). More generally however we might want some way of more continuously varying the effective flexibility of a model with a fixed architecture.

A common method for doing this is to add an additional term to the objective function being minimised during training which penalises some measure of the complexity of a model as a function of the model parameters. The aim of training is then to minimise with respect to the model parameters the sum $E^\star$ of the data-driven error function term $\bar{E}$ and a model complexity term $C$.

\begin{equation}
   E^\star =
   \underbrace{\bar{E}}_{\textrm{data term}} + \underbrace{C}_{\textrm{complexity term}}
\end{equation}

We need the complexity term $C$ to be easy to compute and differentiable with respect to the model parameters. A common choice is to use terms involving the *norms* ($\sim$ a measure of size) of the parameters. This penalises models with large parameter values. Two commonly used norms are the L1 and L2 norms.

For a $D$ dimensional vector $\mathbf{v}$ the L1 norm is defined as

\begin{equation}
\| \boldsymbol{v} \|_1 = \sum_{d=1}^D \left| v_d \right|,
\end{equation}

and the L2 norm is defined as

\begin{equation}
\| \boldsymbol{v} \|_2 = \left[ \sum_{d=1}^D \left( v_d^2 \right) \right]^{\frac{1}{2}}.
\end{equation}

For a $K \times D$ matrix $\mathbf{M}$, we will define norms by collapsing the matrix to a vector $\boldsymbol{m} = \mathrm{vec}\left[\mathbf{M}\right] = \left[ M_{1,1} \dots M_{1,D} ~ M_{2,1} \dots M_{K,D} \right]^{\rm T}$ and then taking the norm as defined above of this resulting vector (practically this just results in summing over two sets of indices rather than one).

The overall complexity penalty term $C$ is defined as a sum over individual complexity terms for each of the $P$ parameters of the model

\begin{equation}
    C = \sum_{i=1}^P \left[ C^{(i)} \right]
\end{equation}

Some of these per-parameter penalty terms $C^{(i)}$ may be set to zero if we do not wish to penalise the size of the corresponding parameter.

To enable us to tradeoff between the model complexity and data error terms, it is typical to introduce positive scalar coefficients $\beta_i$ to scale the penalty term on the $i$th parameter. A *L1 penalty* on the $i$th vector parameter $\boldsymbol{p}^{(i)}$ (or matrix parameter collapsed to a vector) is then commonly defined as

\begin{equation}
  C^{(i)}_{\textrm{L1}} =
  \beta_i \left\| \boldsymbol{p}^{(i)} \right\|_1 =
  \beta_i  \sum_{d=1}^D \left| p^{(i)}_d \right|.
\end{equation}

This has a gradient with respect to the parameter vector

\begin{equation}
  \frac{\partial C^{(i)}_{\textrm{L1}}}{\partial p^{(i)}_d} = \beta_i \, \textrm{sgn}\left( p^{(i)}_d \right)
\end{equation}

where $\textrm{sgn}(u) = +1$ if $u > 0$, $\textrm{sgn}(u) = -1$ if $u < 0$ (and is not well defined for $u=0$ though a common convention is to have $\textrm{sgn}(0) = 0$).

Similarly a *L2 penalty* on the $i$th vector parameter $\boldsymbol{p}^{(i)}$ (or matrix parameter collapsed to a vector) is commonly defined as

\begin{equation}
  C^{(i)}_{\textrm{L2}} =
  \frac{1}{2} \beta_i \left\| \boldsymbol{p}^{(i)} \right\|_2^2 =
  \frac{1}{2} \beta_i \sum_{d=1}^D \left[ \left( p^{(i)}_d \right)^2 \right].
\end{equation}

Somewhat confusingly this is proportional to the square of the L2 norm rather than the L2 norm itself, however it is an almost universal convention to call this an L2 penalty so we will stick with this nomenclature here. The $\frac{1}{2}$ term is less universal and is sometimes not included; we include it here for consistency with how we defined the sum of squared errors cost. Similarly to that case, the $\frac{1}{2}$ cancels when calculating the gradient with respect to the parameter

\begin{equation}
  \frac{\partial C^{(i)}_{\textrm{L2}}}{\partial p^{(i)}_d} = \beta_i p^{(i)}_d
\end{equation}

Use the above definitions for the L1 and L2 penalties for a parameter and corresponding gradients to implement the `__call__` and `grad` methods respectively for the skeleton `L1Penalty` and `L2Penalty` class definitions below. The `coefficient` propert of these classes should be used as the $\beta_i$ value in the equations above. The parameter the penalty term (or gradient) is being evaluated for will be either a one or two-dimensional NumPy array (corresponding to a vector or matrix parameter respectively) and your implementations should be able to cope with both cases.

```python
import numpy as np
seed = 22102017
rng = np.random.RandomState(seed)
class L1Penalty(object):
    """L1 parameter penalty.

    Term to add to the objective function penalising parameters
    based on their L1 norm.
    """

    def __init__(self, coefficient):
        """Create a new L1 penalty object.

        Args:
            coefficient: Positive constant to scale penalty term by.
        """
        assert coefficient > 0., 'Penalty coefficient must be positive.'
        self.coefficient = coefficient

    def __call__(self, parameter):
        """Calculate L1 penalty value for a parameter.

        Args:
            parameter: Array corresponding to a model parameter.

        Returns:
            Value of penalty term.
        """
        return self.coefficient * abs(parameter).sum()

    def grad(self, parameter):
        """Calculate the penalty gradient with respect to the parameter.

        Args:
            parameter: Array corresponding to a model parameter.

        Returns:
            Value of penalty gradient with respect to parameter. This
            should be an array of the same shape as the parameter.
        """
        return self.coefficient * np.sign(parameter)

    def __repr__(self):
        return 'L1Penalty({0})'.format(self.coefficient)

class L2Penalty(object):
    """L1 parameter penalty.

    Term to add to the objective function penalising parameters
    based on their L2 norm.
    """

    def __init__(self, coefficient):
        """Create a new L2 penalty object.

        Args:
            coefficient: Positive constant to scale penalty term by.
        """
        assert coefficient > 0., 'Penalty coefficient must be positive.'
        self.coefficient = coefficient

    def __call__(self, parameter):
        """Calculate L2 penalty value for a parameter.

        Args:
            parameter: Array corresponding to a model parameter.

        Returns:
            Value of penalty term.
        """
        return 0.5 * self.coefficient * (parameter**2).sum()

    def grad(self, parameter):
        """Calculate the penalty gradient with respect to the parameter.

        Args:
            parameter: Array corresponding to a model parameter.

        Returns:
            Value of penalty gradient with respect to parameter. This
            should be an array of the same shape as the parameter.
        """
        return self.coefficient * parameter

    def __repr__(self):
        return 'L2Penalty({0})'.format(self.coefficient)
```

Test your implementations by running the cells below.

```python
test_params_1 = np.array([[0.5, 0.3, -1.2, 5.8], [0.2, -3.1, 4.9, -5.0]])
test_params_2 = np.array([0.8, -0.6, -0.3, 1.5, 2.8])
true_l1_cost_1 = 10.5
true_l1_grad_1 = np.array([[0.5, 0.5, -0.5, 0.5], [0.5, -0.5, 0.5, -0.5]])
true_l1_cost_2 = 3.
true_l1_grad_2 = np.array([0.5, -0.5, -0.5, 0.5, 0.5])
l1 = L1Penalty(0.5)
if (not np.allclose(l1(test_params_1), true_l1_cost_1) or
        not np.allclose(l1(test_params_2), true_l1_cost_2)):
    print('L1Penalty.__call__ giving incorrect value(s).')
elif (not np.allclose(l1.grad(test_params_1), true_l1_grad_1) or
          not np.allclose(l1.grad(test_params_2), true_l1_grad_2)):
    print('L1Penalty.grad giving incorrect value(s).')
else:
    print('All test values calculated correctly for L1Penalty')
```

    All test values calculated correctly for L1Penalty

```python
test_params_1 = np.array([[0.5, 0.3, -1.2, 5.8], [0.2, -3.1, 4.9, -5.0]])
test_params_2 = np.array([0.8, -0.6, -0.3, 1.5, 2.8])
true_l2_cost_1 = 23.52
true_l2_grad_1 = np.array([[0.25, 0.15, -0.6, 2.9], [0.1, -1.55, 2.45, -2.5]])
true_l2_cost_2 = 2.795
true_l2_grad_2 = np.array([0.4, -0.3, -0.15, 0.75, 1.4])
l2 = L2Penalty(0.5)
if (not np.allclose(l2(test_params_1), true_l2_cost_1) or
        not np.allclose(l2(test_params_2), true_l2_cost_2)):
    print('L2Penalty.__call__ giving incorrect value(s).')
elif (not np.allclose(l2.grad(test_params_1), true_l2_grad_1) or
          not np.allclose(l2.grad(test_params_2), true_l2_grad_2)):
    print('L2Penalty.grad giving incorrect value(s).')
else:
    print('All test values calculated correctly for L2Penalty')
```

    All test values calculated correctly for L2Penalty

## Exercise 2: Training with regularisation

Previously in the second laboratory you implemented a function `grads_wrt_params`  to calculate the gradient of an error function with respect to the parameters of an affine model (layer), given gradients of the error function with respect to the model (layer) outputs.

If we are training a model using a regularised objective function, we need to additionally calculate the gradients of the regularisation penalty terms with respect to the parameters and add these to the error function gradient terms. Following from the definition of the regularised objective $E^\star$ above we have that the gradient of the overall objective with respect to the $d$th element of the $i$th model parameter is

\begin{equation}
    \frac{\partial E^\star}{\partial p^{(i)}_d} =
    \frac{\partial \bar{E}}{\partial p^{(i)}_d} +
    \frac{\partial C}{\partial p^{(i)}_d}
\end{equation}

We have already discussed in the second lab notebook how to calculate the error function gradient term $\frac{\partial \bar{E}}{\partial p^{(i)}_d}$. As the model complexity term is composed of a sum of per parameter terms and only one of these will depend on the $i$th parameter we can write

\begin{equation}
    \frac{\partial C}{\partial p^{(i)}_d} = \frac{\partial C^{(i)}}{\partial p^{(i)}_d}
\end{equation}

which corresponds to the penalty term gradients you implemented above. To enable us to use the same `Optimiser` implementation that we have previously used to train models without regularisation, we have altered the implementation of the `AffineLayer` class (this being the only layer we currently have defined with parameters) to allow us to specify penalty terms on the weight matrix and bias vector when creating an instance of the class and to add the corresponding penalty gradients to the returned value from the `grads_wrt_params` method.

The penalty terms need to be specified as a class matching the interface of the `L1Penalty` and `L2Penalty` classes you implemented above, defining both a `__call__` method to calculate the penalty value for a parameter and a `grad` method to calculate the gradient of the penalty with respect to the parameter. Separate penalties can be specified for the weight and bias parameters, with it common to only regularise the weight parameters.

The penalty terms for a layer are specifed using the `weights_penalty` and `biases_penalty` arguments to the `__init__` method of the `AffineLayer` class. If either (or both) ofthese are set to `None` (the default) no regularisation is applied to the corresponding parameter.

Using the `L1Penalty` and `L2Penalty` classes you implemented in the previous exercise, train models to classify MNIST digit images with

  * no regularisation
  * an L1 penalty with coefficient $10^{-5}$ on the all of the weight matrix parameters
  * an L1 penalty with coefficient $10^{-3}$ on the all of the weight matrix parameters
  * an L2 penalty with coefficient $10^{-4}$ on the all of the weight matrix parameters
  * an L2 penalty with coefficient $10^{-2}$ on the all of the weight matrix parameters

The models should all have three affine layers interspersed with rectified linear layers (as implemented in the first exercise) and intermediate layers between the input and output should have dimensionalities of 100. The final output layer should be an `AffineLayer` (the model outputting the logarithms of the non-normalised class probabilties) and you should use the `CrossEntropySoftmaxError` as the error function (which calculates the softmax of the model outputs to convert to normalised class probabilities before calculating the corresponding multi-class cross entropy error).

Use the `GlorotInit` class introduced in the first coursework to initialise the weights in all layers, using a gain of 0.5 (this adjusts for the fact that the rectified linear sets zeros all negative inputs), and initialises the biases to zero with a `ConstantInit` object.

As an example the necessary parameter initialisers, model and error can be defined using

```python
weights_init = GlorotUniformInit(0.5, rng)
biases_init = ConstantInit(0.)
input_dim, output_dim, hidden_dim = 784, 10, 100
model = MultipleLayerModel([
    AffineLayer(input_dim, hidden_dim, weights_init,
                biases_init, weights_penalty=weights_penalty),
    ReluLayer(),
    AffineLayer(hidden_dim, hidden_dim, weights_init,
                biases_init, weights_penalty=weights_penalty),
    ReluLayer(),
    AffineLayer(hidden_dim, output_dim, weights_init,
                biases_init, weights_penalty=weights_penalty)
])
error = CrossEntropySoftmaxError()
```

This assumes all the relevant classes have been imported from their modules, a penalty object has been assigned to `weights_penalty` and a seeded random number generator assigned to `rng`.

For each regularisation scheme, train the model for 100 epochs with a batch size of 50 and using a gradient descent with momentum learning rule with learning rate 0.01 and momentum coefficient 0.9. For each regularisation scheme you should store the run statistics (output of `Optimiser.train`) and the final values of the first layer weights for each of the trained models.

```python
import numpy as np
from collections import OrderedDict
import logging
from mlp.layers import AffineLayer, ReluLayer
from mlp.errors import CrossEntropySoftmaxError
from mlp.models import MultipleLayerModel
from mlp.initialisers import GlorotUniformInit, ConstantInit
from mlp.learning_rules import MomentumLearningRule
from mlp.data_providers import MNISTDataProvider
from mlp.optimisers import Optimiser

# Seed a random number generator
seed = 24102016
rng = np.random.RandomState(seed)

num_epochs = 100
stats_interval = 5
batch_size = 50
learning_rate = 0.01
mom_coeff = 0.9
weights_init_gain = 0.5
biases_init = 0.

# Set up a logger object to print info about the training run to stdout
logger = logging.getLogger()
logger.setLevel(logging.INFO)
logger.handlers = [logging.StreamHandler()]

# Create data provider objects for the MNIST data set
train_data = MNISTDataProvider('train', batch_size, rng=rng)
valid_data = MNISTDataProvider('valid', batch_size, rng=rng)
input_dim, output_dim, hidden_dim = 784, 10, 100

weights_init = GlorotUniformInit(weights_init_gain, rng)
biases_init = ConstantInit(biases_init)
error = CrossEntropySoftmaxError()
learning_rule = MomentumLearningRule(learning_rate, mom_coeff)
data_monitors = {'acc': lambda y, t: (y.argmax(-1) == t.argmax(-1)).mean()}

```

```python
weights_penalties = [
    None,
    L1Penalty(1e-5),
    L1Penalty(1e-3),
    L2Penalty(1e-4),
    L2Penalty(1e-2)
]
run_info = OrderedDict()
models = OrderedDict()
for weights_penalty in weights_penalties:
    # Reset random number generator and data provider states on each run
    # to ensure reproducibility of results
    rng.seed(seed)
    train_data.reset()
    valid_data.reset()
    print('Regularisation: {0}'.format(weights_penalty))
    model = MultipleLayerModel([
        AffineLayer(input_dim, hidden_dim, weights_init,
                    biases_init, weights_penalty),
        ReluLayer(),
        AffineLayer(hidden_dim, hidden_dim, weights_init,
                    biases_init, weights_penalty),
        ReluLayer(),
        AffineLayer(hidden_dim, output_dim, weights_init,
                    biases_init, weights_penalty)
    ])
    optimiser = Optimiser(model, error, learning_rule, train_data, valid_data, data_monitors)
    run_info[weights_penalty] = optimiser.train(num_epochs, stats_interval)
    models[weights_penalty] = model
```

    Regularisation: None

    Epoch 5: 10.9s to complete
        error(train)=6.34e-02, acc(train)=9.82e-01, error(valid)=9.67e-02, acc(valid)=9.72e-01
    Epoch 10: 10.1s to complete
        error(train)=2.36e-02, acc(train)=9.94e-01, error(valid)=8.34e-02, acc(valid)=9.77e-01
    Epoch 15: 11.7s to complete
        error(train)=8.49e-03, acc(train)=9.98e-01, error(valid)=8.34e-02, acc(valid)=9.78e-01
    Epoch 20: 10.5s to complete
        error(train)=3.42e-03, acc(train)=1.00e+00, error(valid)=9.58e-02, acc(valid)=9.78e-01
    Epoch 25: 10.7s to complete
        error(train)=1.42e-03, acc(train)=1.00e+00, error(valid)=9.48e-02, acc(valid)=9.79e-01
    Epoch 30: 11.4s to complete
        error(train)=9.06e-04, acc(train)=1.00e+00, error(valid)=1.01e-01, acc(valid)=9.79e-01
    Epoch 35: 12.0s to complete
        error(train)=6.40e-04, acc(train)=1.00e+00, error(valid)=1.04e-01, acc(valid)=9.79e-01
    Epoch 40: 10.9s to complete
        error(train)=5.06e-04, acc(train)=1.00e+00, error(valid)=1.06e-01, acc(valid)=9.79e-01
    Epoch 45: 15.1s to complete
        error(train)=4.23e-04, acc(train)=1.00e+00, error(valid)=1.08e-01, acc(valid)=9.79e-01
    Epoch 50: 13.2s to complete
        error(train)=3.59e-04, acc(train)=1.00e+00, error(valid)=1.10e-01, acc(valid)=9.79e-01
    Epoch 55: 10.1s to complete
        error(train)=3.07e-04, acc(train)=1.00e+00, error(valid)=1.11e-01, acc(valid)=9.79e-01
    Epoch 60: 10.3s to complete
        error(train)=2.71e-04, acc(train)=1.00e+00, error(valid)=1.13e-01, acc(valid)=9.79e-01
    Epoch 65: 11.4s to complete
        error(train)=2.42e-04, acc(train)=1.00e+00, error(valid)=1.15e-01, acc(valid)=9.79e-01
    Epoch 70: 10.8s to complete
        error(train)=2.19e-04, acc(train)=1.00e+00, error(valid)=1.16e-01, acc(valid)=9.79e-01
    Epoch 75: 13.8s to complete
        error(train)=1.97e-04, acc(train)=1.00e+00, error(valid)=1.17e-01, acc(valid)=9.79e-01
    Epoch 80: 11.4s to complete
        error(train)=1.82e-04, acc(train)=1.00e+00, error(valid)=1.18e-01, acc(valid)=9.79e-01
    Epoch 85: 15.5s to complete
        error(train)=1.67e-04, acc(train)=1.00e+00, error(valid)=1.19e-01, acc(valid)=9.79e-01
    Epoch 90: 15.0s to complete
        error(train)=1.54e-04, acc(train)=1.00e+00, error(valid)=1.20e-01, acc(valid)=9.79e-01
    Epoch 95: 11.7s to complete
        error(train)=1.44e-04, acc(train)=1.00e+00, error(valid)=1.22e-01, acc(valid)=9.79e-01
    Epoch 100: 12.6s to complete
        error(train)=1.34e-04, acc(train)=1.00e+00, error(valid)=1.22e-01, acc(valid)=9.79e-01

    Regularisation: L1Penalty(1e-05)

    Epoch 5: 25.2s to complete
        error(train)=6.46e-02, acc(train)=9.82e-01, error(valid)=9.39e-02, acc(valid)=9.73e-01
    Epoch 10: 15.4s to complete
        error(train)=2.68e-02, acc(train)=9.93e-01, error(valid)=8.04e-02, acc(valid)=9.78e-01
    Epoch 15: 17.1s to complete
        error(train)=1.42e-02, acc(train)=9.97e-01, error(valid)=8.12e-02, acc(valid)=9.78e-01
    Epoch 20: 15.4s to complete
        error(train)=6.73e-03, acc(train)=9.99e-01, error(valid)=8.81e-02, acc(valid)=9.77e-01
    Epoch 25: 43.1s to complete
        error(train)=3.55e-03, acc(train)=1.00e+00, error(valid)=8.79e-02, acc(valid)=9.78e-01
    Epoch 30: 18.9s to complete
        error(train)=2.77e-03, acc(train)=1.00e+00, error(valid)=9.02e-02, acc(valid)=9.78e-01
    Epoch 35: 23.9s to complete
        error(train)=2.19e-03, acc(train)=1.00e+00, error(valid)=8.69e-02, acc(valid)=9.79e-01
    Epoch 40: 22.3s to complete
        error(train)=2.16e-03, acc(train)=1.00e+00, error(valid)=8.94e-02, acc(valid)=9.79e-01
    Epoch 45: 17.4s to complete
        error(train)=1.91e-03, acc(train)=1.00e+00, error(valid)=8.66e-02, acc(valid)=9.78e-01
    Epoch 50: 22.0s to complete
        error(train)=1.76e-03, acc(train)=1.00e+00, error(valid)=8.73e-02, acc(valid)=9.79e-01
    Epoch 55: 16.9s to complete
        error(train)=2.14e-03, acc(train)=1.00e+00, error(valid)=9.00e-02, acc(valid)=9.79e-01
    Epoch 60: 21.0s to complete
        error(train)=1.98e-03, acc(train)=1.00e+00, error(valid)=8.77e-02, acc(valid)=9.80e-01
    Epoch 65: 15.8s to complete
        error(train)=1.09e-02, acc(train)=9.96e-01, error(valid)=1.07e-01, acc(valid)=9.75e-01
    Epoch 70: 14.3s to complete
        error(train)=9.57e-04, acc(train)=1.00e+00, error(valid)=8.86e-02, acc(valid)=9.80e-01
    Epoch 75: 18.5s to complete
        error(train)=9.22e-04, acc(train)=1.00e+00, error(valid)=8.57e-02, acc(valid)=9.80e-01
    Epoch 80: 20.6s to complete
        error(train)=1.03e-03, acc(train)=1.00e+00, error(valid)=8.58e-02, acc(valid)=9.81e-01
    Epoch 85: 15.3s to complete
        error(train)=1.20e-03, acc(train)=1.00e+00, error(valid)=8.28e-02, acc(valid)=9.80e-01
    Epoch 90: 14.0s to complete
        error(train)=1.24e-03, acc(train)=1.00e+00, error(valid)=8.32e-02, acc(valid)=9.81e-01
    Epoch 95: 18.4s to complete
        error(train)=1.33e-03, acc(train)=1.00e+00, error(valid)=8.41e-02, acc(valid)=9.80e-01
    Epoch 100: 14.4s to complete
        error(train)=1.70e-03, acc(train)=1.00e+00, error(valid)=8.21e-02, acc(valid)=9.80e-01

    Regularisation: L1Penalty(0.001)

    Epoch 5: 18.9s to complete
        error(train)=2.53e-01, acc(train)=9.31e-01, error(valid)=2.40e-01, acc(valid)=9.35e-01
    Epoch 10: 17.4s to complete
        error(train)=2.20e-01, acc(train)=9.36e-01, error(valid)=2.13e-01, acc(valid)=9.38e-01
    Epoch 15: 22.7s to complete
        error(train)=2.26e-01, acc(train)=9.34e-01, error(valid)=2.14e-01, acc(valid)=9.37e-01
    Epoch 20: 14.2s to complete
        error(train)=2.75e-01, acc(train)=9.15e-01, error(valid)=2.58e-01, acc(valid)=9.22e-01
    Epoch 25: 17.6s to complete
        error(train)=2.31e-01, acc(train)=9.30e-01, error(valid)=2.29e-01, acc(valid)=9.30e-01
    Epoch 30: 21.7s to complete
        error(train)=2.20e-01, acc(train)=9.36e-01, error(valid)=2.12e-01, acc(valid)=9.40e-01
    Epoch 35: 14.0s to complete
        error(train)=2.11e-01, acc(train)=9.37e-01, error(valid)=2.03e-01, acc(valid)=9.40e-01
    Epoch 40: 15.9s to complete
        error(train)=2.06e-01, acc(train)=9.37e-01, error(valid)=2.03e-01, acc(valid)=9.39e-01
    Epoch 45: 17.8s to complete
        error(train)=2.02e-01, acc(train)=9.40e-01, error(valid)=1.99e-01, acc(valid)=9.42e-01
    Epoch 50: 15.7s to complete
        error(train)=2.01e-01, acc(train)=9.40e-01, error(valid)=1.95e-01, acc(valid)=9.45e-01
    Epoch 55: 14.7s to complete
        error(train)=2.07e-01, acc(train)=9.36e-01, error(valid)=2.00e-01, acc(valid)=9.39e-01
    Epoch 60: 21.3s to complete
        error(train)=2.38e-01, acc(train)=9.28e-01, error(valid)=2.34e-01, acc(valid)=9.31e-01
    Epoch 65: 17.6s to complete
        error(train)=2.07e-01, acc(train)=9.38e-01, error(valid)=2.02e-01, acc(valid)=9.41e-01
    Epoch 70: 21.5s to complete
        error(train)=1.98e-01, acc(train)=9.40e-01, error(valid)=1.92e-01, acc(valid)=9.44e-01
    Epoch 75: 14.2s to complete
        error(train)=2.09e-01, acc(train)=9.37e-01, error(valid)=2.08e-01, acc(valid)=9.39e-01
    Epoch 80: 16.7s to complete
        error(train)=2.03e-01, acc(train)=9.40e-01, error(valid)=1.98e-01, acc(valid)=9.42e-01
    Epoch 85: 16.8s to complete
        error(train)=1.98e-01, acc(train)=9.42e-01, error(valid)=1.93e-01, acc(valid)=9.43e-01
    Epoch 90: 23.6s to complete
        error(train)=2.00e-01, acc(train)=9.38e-01, error(valid)=1.97e-01, acc(valid)=9.40e-01
    Epoch 95: 15.5s to complete
        error(train)=1.88e-01, acc(train)=9.44e-01, error(valid)=1.86e-01, acc(valid)=9.46e-01
    Epoch 100: 17.7s to complete
        error(train)=1.97e-01, acc(train)=9.42e-01, error(valid)=1.93e-01, acc(valid)=9.46e-01

    Regularisation: L2Penalty(0.0001)

    Epoch 5: 10.8s to complete
        error(train)=6.49e-02, acc(train)=9.81e-01, error(valid)=9.48e-02, acc(valid)=9.72e-01
    Epoch 10: 15.5s to complete
        error(train)=2.65e-02, acc(train)=9.93e-01, error(valid)=8.14e-02, acc(valid)=9.78e-01
    Epoch 15: 13.5s to complete
        error(train)=1.41e-02, acc(train)=9.97e-01, error(valid)=8.20e-02, acc(valid)=9.78e-01
    Epoch 20: 14.2s to complete
        error(train)=6.66e-03, acc(train)=9.99e-01, error(valid)=8.53e-02, acc(valid)=9.78e-01
    Epoch 25: 11.3s to complete
        error(train)=3.56e-03, acc(train)=1.00e+00, error(valid)=8.34e-02, acc(valid)=9.79e-01
    Epoch 30: 12.6s to complete
        error(train)=2.80e-03, acc(train)=1.00e+00, error(valid)=8.27e-02, acc(valid)=9.79e-01
    Epoch 35: 10.2s to complete
        error(train)=2.53e-03, acc(train)=1.00e+00, error(valid)=8.15e-02, acc(valid)=9.80e-01
    Epoch 40: 13.2s to complete
        error(train)=2.44e-03, acc(train)=1.00e+00, error(valid)=8.21e-02, acc(valid)=9.79e-01
    Epoch 45: 17.8s to complete
        error(train)=2.29e-03, acc(train)=1.00e+00, error(valid)=7.91e-02, acc(valid)=9.80e-01
    Epoch 50: 13.7s to complete
        error(train)=2.24e-03, acc(train)=1.00e+00, error(valid)=8.01e-02, acc(valid)=9.79e-01
    Epoch 55: 12.4s to complete
        error(train)=2.20e-03, acc(train)=1.00e+00, error(valid)=7.90e-02, acc(valid)=9.80e-01
    Epoch 60: 13.7s to complete
        error(train)=2.39e-03, acc(train)=1.00e+00, error(valid)=7.88e-02, acc(valid)=9.80e-01
    Epoch 65: 11.1s to complete
        error(train)=2.12e-03, acc(train)=1.00e+00, error(valid)=7.79e-02, acc(valid)=9.80e-01
    Epoch 70: 11.4s to complete
        error(train)=2.22e-03, acc(train)=1.00e+00, error(valid)=7.74e-02, acc(valid)=9.80e-01
    Epoch 75: 14.5s to complete
        error(train)=2.16e-03, acc(train)=1.00e+00, error(valid)=7.53e-02, acc(valid)=9.80e-01
    Epoch 80: 12.6s to complete
        error(train)=2.19e-03, acc(train)=1.00e+00, error(valid)=7.74e-02, acc(valid)=9.81e-01
    Epoch 85: 11.5s to complete
        error(train)=2.21e-03, acc(train)=1.00e+00, error(valid)=7.49e-02, acc(valid)=9.81e-01
    Epoch 90: 11.4s to complete
        error(train)=2.11e-03, acc(train)=1.00e+00, error(valid)=7.49e-02, acc(valid)=9.81e-01
    Epoch 95: 11.2s to complete
        error(train)=2.12e-03, acc(train)=1.00e+00, error(valid)=7.54e-02, acc(valid)=9.80e-01
    Epoch 100: 15.5s to complete
        error(train)=2.32e-03, acc(train)=1.00e+00, error(valid)=7.57e-02, acc(valid)=9.81e-01

    Regularisation: L2Penalty(0.01)

    Epoch 5: 10.6s to complete
        error(train)=2.36e-01, acc(train)=9.41e-01, error(valid)=2.21e-01, acc(valid)=9.45e-01
    Epoch 10: 10.9s to complete
        error(train)=2.13e-01, acc(train)=9.45e-01, error(valid)=1.99e-01, acc(valid)=9.52e-01
    Epoch 15: 11.8s to complete
        error(train)=2.22e-01, acc(train)=9.44e-01, error(valid)=2.08e-01, acc(valid)=9.48e-01
    Epoch 20: 11.0s to complete
        error(train)=2.15e-01, acc(train)=9.46e-01, error(valid)=2.03e-01, acc(valid)=9.49e-01
    Epoch 25: 15.5s to complete
        error(train)=2.14e-01, acc(train)=9.47e-01, error(valid)=2.03e-01, acc(valid)=9.50e-01
    Epoch 30: 11.8s to complete
        error(train)=2.02e-01, acc(train)=9.51e-01, error(valid)=1.91e-01, acc(valid)=9.55e-01
    Epoch 35: 14.7s to complete
        error(train)=2.11e-01, acc(train)=9.47e-01, error(valid)=2.00e-01, acc(valid)=9.50e-01
    Epoch 40: 12.3s to complete
        error(train)=2.01e-01, acc(train)=9.50e-01, error(valid)=1.92e-01, acc(valid)=9.53e-01
    Epoch 45: 11.9s to complete
        error(train)=2.01e-01, acc(train)=9.52e-01, error(valid)=1.91e-01, acc(valid)=9.55e-01
    Epoch 50: 12.1s to complete
        error(train)=2.05e-01, acc(train)=9.50e-01, error(valid)=1.93e-01, acc(valid)=9.55e-01
    Epoch 55: 11.4s to complete
        error(train)=2.03e-01, acc(train)=9.50e-01, error(valid)=1.92e-01, acc(valid)=9.54e-01
    Epoch 60: 12.8s to complete
        error(train)=2.07e-01, acc(train)=9.49e-01, error(valid)=1.95e-01, acc(valid)=9.54e-01
    Epoch 65: 11.7s to complete
        error(train)=2.10e-01, acc(train)=9.46e-01, error(valid)=1.98e-01, acc(valid)=9.49e-01
    Epoch 70: 12.8s to complete
        error(train)=2.03e-01, acc(train)=9.50e-01, error(valid)=1.91e-01, acc(valid)=9.55e-01
    Epoch 75: 15.1s to complete
        error(train)=2.06e-01, acc(train)=9.50e-01, error(valid)=1.95e-01, acc(valid)=9.54e-01
    Epoch 80: 11.1s to complete
        error(train)=1.98e-01, acc(train)=9.52e-01, error(valid)=1.86e-01, acc(valid)=9.56e-01
    Epoch 85: 12.4s to complete
        error(train)=2.02e-01, acc(train)=9.51e-01, error(valid)=1.91e-01, acc(valid)=9.55e-01
    Epoch 90: 13.2s to complete
        error(train)=2.02e-01, acc(train)=9.51e-01, error(valid)=1.91e-01, acc(valid)=9.55e-01
    Epoch 95: 10.6s to complete
        error(train)=2.00e-01, acc(train)=9.52e-01, error(valid)=1.87e-01, acc(valid)=9.57e-01
    Epoch 100: 11.3s to complete
        error(train)=1.98e-01, acc(train)=9.53e-01, error(valid)=1.86e-01, acc(valid)=9.58e-01

Plot the training set error against epoch number for all different regularisation schemes on the same axis. On a second axis plot the validation set error against epoch number for all the different regularisation schemes. Interpret and comment on what you see.

```python
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use('ggplot')

fig = plt.figure(figsize=(12, 6))
ax1 = fig.add_subplot(1, 2, 1)
ax2 = fig.add_subplot(1, 2, 2)
for weight_penalty, run in run_info.items():
    stats, keys, run_time = run
    ax1.plot(np.arange(1, stats.shape[0]) * stats_interval,
             stats[1:, keys['error(train)']], label=str(weight_penalty))
    ax2.plot(np.arange(1, stats.shape[0]) * stats_interval,
             stats[1:, keys['error(valid)']], label=str(weight_penalty))
ax1.legend(loc=0)
ax1.set_xlabel('Epoch number')
ax1.set_ylabel('Training set error')
ax1.set_yscale('log')
ax2.legend(loc=0)
ax2.set_xlabel('Epoch number')
ax2.set_ylabel('Validation set error')
ax2.set_yscale('log')
```

The cell below defines two functions for visualising the first layer weights of the models trained above. The first plots a histogram of the weight values and the second plots the first layer weights as feature maps, i.e. each row of the first layer weight matrix (corresponding to the weights going from the input MNIST image to a particular first layer output) is visualised as a $28\times 28$ image. In these feature maps white corresponds to negative weights, black to positive weights and grey to weights close to zero.

Use these functions to plot a histogram and feature map visualisation for the first layer weights of each model trained above. You should try to interpret the plots in the context of what you were told in the lecture about the behaviour of L1 versus L2 regularisation.

```python
def plot_param_histogram(param, fig_size=(6, 3), interval=[-1.5, 1.5]):
    """Plots a normalised histogram of an array of parameter values."""
    fig = plt.figure(figsize=fig_size)
    ax = fig.add_subplot(111)
    ax.hist(param.flatten(), 50, interval, normed=True)
    ax.set_xlabel('Parameter value')
    ax.set_ylabel('Normalised frequency density')
    return fig, ax

def visualise_first_layer_weights(weights, fig_size=(5, 5)):
    """Plots a grid of first layer weights as feature maps."""
    fig = plt.figure(figsize=fig_size)
    num_feature_maps = weights.shape[0]
    grid_size = int(num_feature_maps**0.5)
    max_abs = np.abs(model.params[0]).max()
    tiled = -np.ones((30 * grid_size,
                      30 * num_feature_maps // grid_size)) * max_abs
    for i, fm in enumerate(model.params[0]):
        r, c = i % grid_size, i // grid_size
        tiled[1 + r * 30:(r + 1) * 30 - 1,
              1 + c * 30:(c + 1) * 30 - 1] = fm.reshape((28, 28))
    ax = fig.add_subplot(111)
    max_abs = np.abs(tiled).max()
    ax.imshow(tiled, cmap='Greys', vmin=-max_abs, vmax=max_abs)
    ax.axis('off')
    fig.tight_layout()
    plt.show()
    return fig, ax
```

```python
for weight_penalty, model in models.items():
    print('Regularisation: {0}'.format(weight_penalty))
    _ = plot_param_histogram(model.params[0], interval=[-0.5, 0.5])
    _ = visualise_first_layer_weights(model.params[0])

```

## Exercise 4: Random data augmentation

Another technique mentioned in the lectures for trying to reduce overfitting is to artificially augment the training data set by performing random transformations to the original training data inputs. The idea is to produce further artificial inputs corresponding to the same target class as the original input. The closer the artificially generated inputs are to the appearing like the true inputs the better as they provide more realistic additional examples for the model to learn from.

For the handwritten image inputs in the MNIST dataset, an obvious way to considering augmenting the dataset is to apply small rotations to the original images. Providing the rotations are small we would generally expect that what we would identify as the class of a digit image will remain the same.

Implement a function which given a batch of MNIST images as 784-dimensional vectors, i.e. an array of shape `(batch_size, 784)`

  * chooses 25% of the images in the batch at random
  * for each image in the 25% chosen, rotates the image by a random angle in $\left[-30^\circ,\,30^\circ\right]$
  * returns a new array of size `(batch_size, 784)` in which the rows corresponding to the 25% chosen images are the vectors corresponding to the new randomly rotated images, while the remaining rows correspond to the original images.

You will need to make use of the [`scipy.ndimage.interpolation.rotate`](https://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.ndimage.interpolation.rotate.html#scipy.ndimage.interpolation.rotate) function which is imported below for you. For computational efficiency you should use bilinear interpolation by setting `order=1` as a keyword argument to this function rather than using the default of bicubic interpolation. Additionally you should make sure the original shape of the images is maintained by passing a `reshape=False` keyword argument.

```python
from scipy.ndimage.interpolation import rotate

def random_rotation(inputs, rng):
    """Randomly rotates a subset of images in a batch.

    Args:
        inputs: Input image batch, an array of shape (batch_size, 784).
        rng: A seeded random number generator.

    Returns:
        An array of shape (batch_size, 784) corresponding to a copy
        of the original `inputs` array with the randomly selected
        images rotated by a random angle. The original `inputs`
        array should not be modified.
    """
    orig_ims = inputs.reshape((-1, 28, 28))
    new_ims = orig_ims.copy()
    indices = rng.choice(orig_ims.shape[0], orig_ims.shape[0] // 4, False)
    angles = rng.uniform(-1., 1., size=indices.shape[0]) * 30.
    for i, j in enumerate(indices):
        new_ims[j] = rotate(orig_ims[j], angles[i], order=1, reshape=False)
    return new_ims.reshape((-1, 784))
```

Use the cell below to test your implementation. This uses the `show_batch_of_images` function we implemented in the first lab notebook to visualise the images in a batch before and after application of the random rotation transformation.

```python
from mlp.data_providers import MNISTDataProvider
import matplotlib.pyplot as plt
def show_batch_of_images(img_batch, fig_size=(3, 3)):
    fig = plt.figure(figsize=fig_size)
    batch_size, im_height, im_width = img_batch.shape
    # calculate no. columns per grid row to give square grid
    grid_size = int(batch_size**0.5)
    # intialise empty array to tile image grid into
    tiled = np.empty((im_height * grid_size,
                      im_width * batch_size // grid_size))
    # iterate over images in batch + indexes within batch
    for i, img in enumerate(img_batch):
        # calculate grid row and column indices
        r, c = i % grid_size, i // grid_size
        tiled[r * im_height:(r + 1) * im_height,
              c * im_height:(c + 1) * im_height] = img
    ax = fig.add_subplot(111)
    ax.imshow(tiled, cmap='Greys') #, vmin=0., vmax=1.)
    ax.axis('off')
    fig.tight_layout()
    plt.show()
    return fig, ax

test_data = MNISTDataProvider('test', 100, rng=rng)
inputs, targets = test_data.next()
_ = show_batch_of_images(inputs.reshape((-1, 28, 28)))
transformed_inputs = random_rotation(inputs, rng)
_ = show_batch_of_images(transformed_inputs.reshape((-1, 28, 28)))
```

## Exercise 5: Training with data augmentation

One simple way to use data augmentation is to just statically augment the training data set - for example we could iterate through the training dataset applying a transformation function like that implemented above to generate new artificial inputs, and use both the original and newly generated data in a new data provider object. We are quite limited however in how far we can augment the dataset by with a static method like this however - if we wanted to apply 9 random rotations to each image in the original datase, we would end up with a dataset with 10 times the memory requirements and that would take 10 times as long to run through each epoch.

An alternative is to randomly augment the data on the fly as we iterate through the data provider in each epoch. In this method a new data provider class can be defined that inherits from the original data provider to be augmented, and provides a new `next` method which applies a random transformation function like that implemented in the previous exercise to each input batch before returning it. This method means that on every epoch a different set of training examples are provided to the model and so in some ways corresponds to an 'infinite' data set (although the amount of variability in the dataset will still be significantly less than the variability in all possible digit images). Compared to static augmentation, this dynamic augmentation scheme comes at the computational cost of having to apply the random transformation each time a new batch is provided. We can vary this overhead by changing the proportion of images in a batch randomly transformed.

An implementation of this scheme has been provided for the MNIST data set in the `AugmentedMNISTDataProvider` object in the `mlp.data_providers` module. In addition to the arguments of the original `MNISTDataProvider.__init__` method, this additional takes a `transformer` argument, which should be a function which takes as arguments an inputs batch array and a random number generator object, and returns an array corresponding to a random transformation of the inputs.

Train a model with the same architecture as in exercise 3 and with no L1 / L2 regularisation using a training data provider which randomly augments the training images using your `random_rotation` transformer function. Plot the training and validation set errors over the training epochs and compare this plot to your previous results from exercise 3.

```python
from mlp.data_providers import AugmentedMNISTDataProvider

aug_train_data = AugmentedMNISTDataProvider('train', rng=rng, transformer=random_rotation)
```

```python
batch_size = 100
num_epochs = 100
learning_rate = 0.01
mom_coeff = 0.9
stats_interval = 5

rng.seed(seed)
aug_train_data.reset()
valid_data.reset()
aug_train_data.batch_size = batch_size
valid_data.batch_size = batch_size

weights_init = GlorotUniformInit(0.5, rng=rng)
biases_init = ConstantInit(0.)

model = MultipleLayerModel([
    AffineLayer(input_dim, hidden_dim, weights_init, biases_init),
    ReluLayer(),
    AffineLayer(hidden_dim, hidden_dim, weights_init, biases_init),
    ReluLayer(),
    AffineLayer(hidden_dim, output_dim, weights_init, biases_init)
])

error = CrossEntropySoftmaxError()
learning_rule = MomentumLearningRule(learning_rate=learning_rate, mom_coeff=mom_coeff)
data_monitors={'acc': lambda y, t: (y.argmax(-1) == t.argmax(-1)).mean()}
optimiser = Optimiser(
    model, error, learning_rule, aug_train_data, valid_data, data_monitors)

aug_stats, aug_keys, aug_run_time = optimiser.train(
    num_epochs=num_epochs, stats_interval=stats_interval)
```

    Epoch 5: 18.6s to complete
        error(train)=1.45e-01, acc(train)=9.57e-01, error(valid)=1.23e-01, acc(valid)=9.65e-01
    Epoch 10: 19.2s to complete
        error(train)=8.53e-02, acc(train)=9.74e-01, error(valid)=8.48e-02, acc(valid)=9.74e-01
    Epoch 15: 19.3s to complete
        error(train)=6.26e-02, acc(train)=9.81e-01, error(valid)=7.69e-02, acc(valid)=9.77e-01
    Epoch 20: 22.9s to complete
        error(train)=5.26e-02, acc(train)=9.85e-01, error(valid)=7.79e-02, acc(valid)=9.76e-01
    Epoch 25: 16.4s to complete
        error(train)=4.33e-02, acc(train)=9.88e-01, error(valid)=7.29e-02, acc(valid)=9.78e-01
    Epoch 30: 19.1s to complete
        error(train)=3.72e-02, acc(train)=9.89e-01, error(valid)=7.23e-02, acc(valid)=9.78e-01
    Epoch 35: 19.4s to complete
        error(train)=3.56e-02, acc(train)=9.90e-01, error(valid)=7.50e-02, acc(valid)=9.79e-01
    Epoch 40: 18.6s to complete
        error(train)=2.90e-02, acc(train)=9.93e-01, error(valid)=7.00e-02, acc(valid)=9.79e-01
    Epoch 45: 18.8s to complete
        error(train)=2.77e-02, acc(train)=9.93e-01, error(valid)=6.72e-02, acc(valid)=9.80e-01
    Epoch 50: 19.0s to complete
        error(train)=3.12e-02, acc(train)=9.91e-01, error(valid)=7.63e-02, acc(valid)=9.80e-01
    Epoch 55: 20.1s to complete
        error(train)=2.23e-02, acc(train)=9.94e-01, error(valid)=6.63e-02, acc(valid)=9.83e-01
    Epoch 60: 25.8s to complete
        error(train)=2.10e-02, acc(train)=9.94e-01, error(valid)=6.99e-02, acc(valid)=9.81e-01
    Epoch 65: 20.6s to complete
        error(train)=1.95e-02, acc(train)=9.95e-01, error(valid)=6.53e-02, acc(valid)=9.83e-01
    Epoch 70: 17.1s to complete
        error(train)=2.11e-02, acc(train)=9.94e-01, error(valid)=6.88e-02, acc(valid)=9.81e-01
    Epoch 75: 24.2s to complete
        error(train)=1.93e-02, acc(train)=9.95e-01, error(valid)=6.44e-02, acc(valid)=9.83e-01
    Epoch 80: 20.2s to complete
        error(train)=1.78e-02, acc(train)=9.95e-01, error(valid)=6.52e-02, acc(valid)=9.83e-01
    Epoch 85: 20.1s to complete
        error(train)=1.68e-02, acc(train)=9.96e-01, error(valid)=6.30e-02, acc(valid)=9.84e-01
    Epoch 90: 19.9s to complete
        error(train)=2.17e-02, acc(train)=9.94e-01, error(valid)=7.03e-02, acc(valid)=9.82e-01
    Epoch 95: 18.5s to complete
        error(train)=1.72e-02, acc(train)=9.96e-01, error(valid)=6.59e-02, acc(valid)=9.82e-01
    Epoch 100: 20.6s to complete
        error(train)=2.21e-02, acc(train)=9.93e-01, error(valid)=8.22e-02, acc(valid)=9.80e-01

```python
fig = plt.figure(figsize=(12, 12))
ax1 = fig.add_subplot(2, 1, 1)
ax2 = fig.add_subplot(2, 1, 2)
for weight_penalty, run in run_info.items():
    stats, keys, run_time = run
    ax1.plot(np.arange(1, stats.shape[0]) * stats_interval,
             stats[1:, keys['error(train)']], label=str(weight_penalty))
    ax2.plot(np.arange(1, stats.shape[0]) * stats_interval,
             stats[1:, keys['error(valid)']], label=str(weight_penalty))
ax1.plot(np.arange(1, aug_stats.shape[0]) * stats_interval,
         aug_stats[1:, aug_keys['error(train)']], label='Data augmentation')
ax2.plot(np.arange(1, aug_stats.shape[0]) * stats_interval,
         aug_stats[1:, aug_keys['error(valid)']], label='Data augmentation')
ax1.legend(loc=0)
ax1.set_xlabel('Epoch number')
ax1.set_ylabel('Training set error')
ax1.set_yscale('log')
ax2.legend(loc=0)
ax2.set_xlabel('Epoch number')
ax2.set_ylabel('Validation set error')
ax2.set_yscale('log')
```

```python

```

```python

```

```python

```

max_stars_repo_name: Amritds/mlpractical

max_issues_repo_name: Amritds/mlpractical

max_forks_repo_name: Amritds/mlpractical


text:
```python
import sympy as sym
from sympy.plotting import plot
from sympy import simplify

x = sym.Symbol('x')
ex = (x**2-2*x-6) / (x-1)
a = sym.diff(ex,x)
# print(a)
# plot(a)
result = a.subs(x,0)
print("답", result)
simplify(a)

```

    답 8

$\displaystyle \frac{x^{2} - 2 x + 8}{x^{2} - 2 x + 1}$

답 8

```python
n = 0.0001
r = ((n**2-2*n-6) / (n - 1)  - 6) / n
print(r)
```

    8.000700070001088

```python

```

max_stars_repo_name: AramMoon/Math

max_issues_repo_name: AramMoon/Math

max_forks_repo_name: AramMoon/Math