<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Sebastian Raschka
- Link to the containing GitHub Repository
Step1:  <font size="1.5em">More information about the watermark magic command extension.</font>
<hr>
I would be happy to hear your comments and suggestions. 
Please feel free to drop me a note via
twitter, email, or google+.
<hr>
Linear Discriminant Analysis bit by bit
<br>
<br>
Sections
Introduction
Principal Component Analysis vs. Linear Discriminant Analysis
What is a "good" feature subspace?
Summarizing the LDA approach in 5 steps
Preparing the sample data set
About the Iris dataset
Reading in the dataset
Histograms and feature selection
Standardization
Normality assumptions
LDA in 5 steps
Step 1
Step2:  <a name="sample_data"></a>
<br>
<br>
Reading in the dataset
[back to top]
Step3:  $\pmb X = \begin{bmatrix} x_{1_{\text{sepal length}}} & x_{1_{\text{sepal width}}} & x_{1_{\text{petal length}}} & x_{1_{\text{petal width}}}\ 
x_{2_{\text{sepal length}}} & x_{2_{\text{sepal width}}} & x_{2_{\text{petal length}}} & x_{2_{\text{petal width}}}\ 
...  \
x_{150_{\text{sepal length}}} & x_{150_{\text{sepal width}}} & x_{150_{\text{petal length}}} & x_{150_{\text{petal width}}}\ 
\end{bmatrix}, \;\;
\pmb y = \begin{bmatrix} \omega_{\text{setosa}}\ 
\omega_{\text{setosa}}\ 
...  \
\omega_{\text{virginica}}\end{bmatrix}$
<a name="sample_data"></a>
<br>
<br>
Since it is more convenient to work with numerical values, we will use the LabelEncode from the scikit-learn library to convert the class labels into numbers
Step4:  $\pmb y = \begin{bmatrix}{\text{setosa}}\ 
{\text{setosa}}\ 
...  \
{\text{virginica}}\end{bmatrix} \quad \Rightarrow
\begin{bmatrix} {\text{1}}\ 
{\text{1}}\ 
...  \
{\text{3}}\end{bmatrix}$
<a name="sample_data"></a>
<br>
<br>
Histograms and feature selection
[back to top]
Just to get a rough idea how the samples of our three classes $\omega_1$, $\omega_2$ and $\omega_3$ are distributed, let us visualize the distributions of the four different features in 1-dimensional histograms.
Step5:  From just looking at these simple graphical representations of the features, we can already tell that the petal lengths and widths are likely better suited as potential features two separate between the three flower classes. In practice, instead of reducing the dimensionality via a projection (here
Step6:  <a name="sample_data"></a>
<br>
<br>
Normality assumptions
[back to top]
It should be mentioned that LDA assumes normal distributed data, features that are statistically independent, and identical covariance matrices for every class. However, this only applies for LDA as classifier and LDA for dimensionality reduction can also work reasonably well if those assumptions are violated. And even for classification tasks LDA seems can be quite robust to the distribution of the data
Step7:  <a name="sample_data"></a>
<br>
<br>
<a name="sc_matrix"></a>
Step 2
Step8:  <br>
2.1 b
Alternatively, we could also compute the class-covariance matrices by adding the scaling factor $\frac{1}{N-1}$ to the within-class scatter matrix, so that our equation becomes
$\Sigma_i = \frac{1}{N_{i}-1} \sum\limits_{\pmb x \in D_i}^n (\pmb x - \pmb m_i)\;(\pmb x - \pmb m_i)^T$.
and $S_W = \sum\limits_{i=1}^{c} (N_{i}-1) \Sigma_i$
where $N_{i}$ is the sample size of the respective class (here
Step9:  <br>
<br>
Step 3
Step10:  <br>
<br>
After this decomposition of our square matrix into eigenvectors and eigenvalues, let us briefly recapitulate how we can interpret those results. As we remember from our first linear algebra class in high school or college, both eigenvectors and eigenvalues are providing us with information about the distortion of a linear transformation
Step11:  <br>
<br>
Step 4
Step12:  <br>
<br>
If we take a look at the eigenvalues, we can already see that 2 eigenvalues are close to 0 and conclude that the eigenpairs are less informative than the other two. Let's express the "explained variance" as percentage
Step13:  <br>
<br>
The first eigenpair is by far the most informative one, and we won't loose much information if we would form a 1D-feature spaced based on this eigenpair.
<br>
<br>
4.2. Choosing k eigenvectors with the largest eigenvalues
[back to top]
After sorting the eigenpairs by decreasing eigenvalues, it is now time to construct our $k \times d$-dimensional eigenvector matrix $\pmb W$ (here $4 \times 2$
Step14:  <br>
<br>
Step 5
Step15:  The scatter plot above represents our new feature subspace that we constructed via LDA. We can see that the first linear discriminant "LD1" separates the classes quite nicely. However, the second discriminant, "LD2", does not add much valuable information, which we've already concluded when we looked at the ranked eigenvalues is step 4.
<br>
<br>
A comparison of PCA and LDA
[back to top]
In order to compare the feature subspace that we obtained via the Linear Discriminant Analysis, we will use the PCA class from the scikit-learn machine-learning library. The documentation can be found here
Step16:  <br>
<br>
The two plots above nicely confirm what we have discussed before
<ASSISTANT_TASK:>
Python Code:
%load_ext watermark
%watermark -v -d -u -p pandas,scikit-learn,numpy,matplotlib 
Explanation: Sebastian Raschka
- Link to the containing GitHub Repository: https://github.com/rasbt/pattern_classification
- Link to this IPython Notebook on GitHub: linear_discriminant_analysis.ipynb
End of explanation
feature_dict = {i:label for i,label in zip(
                range(4),
                  ('sepal length in cm', 
                  'sepal width in cm', 
                  'petal length in cm', 
                  'petal width in cm', ))}
Explanation: <font size="1.5em">More information about the watermark magic command extension.</font>
<hr>
I would be happy to hear your comments and suggestions. 
Please feel free to drop me a note via
twitter, email, or google+.
<hr>
Linear Discriminant Analysis bit by bit
<br>
<br>
Sections
Introduction
Principal Component Analysis vs. Linear Discriminant Analysis
What is a "good" feature subspace?
Summarizing the LDA approach in 5 steps
Preparing the sample data set
About the Iris dataset
Reading in the dataset
Histograms and feature selection
Standardization
Normality assumptions
LDA in 5 steps
Step 1: Computing the d-dimensional mean vectors
Step 2: Computing the Scatter Matrices
Step 3: Solving the generalized eigenvalue problem for the matrix $S_{W}^{-1}S_B$
Step 4: Selecting linear discriminants for the new feature subspace
Step 5: Transforming the samples onto the new subspace
A comparison of PCA and LDA
LDA via scikit-learn
Update-scikit-learn-0.15.2
<br>
<br>
Introduction
[back to top]
Linear Discriminant Analysis (LDA) is most commonly used as dimensionality reduction technique in the pre-processing step for pattern-classification and machine learning applications. 
The goal is to project a dataset onto a lower-dimensional space with good class-separability in order avoid overfitting ("curse of dimensionality") and also reduce computational costs.
Ronald A. Fisher formulated the Linear Discriminant in 1936 (The Use of Multiple  Measurements in Taxonomic Problems), and it also has some practical uses as classifier. The original Linear discriminant was described for a 2-class problem, and it was then later generalized as "multi-class Linear Discriminant Analysis" or "Multiple Discriminant Analysis" by C. R. Rao in 1948 (The utilization of multiple measurements in problems of biological classification)
The general LDA approach is very similar to a Principal Component Analysis (for more information about the PCA, see the previous article Implementing a Principal Component Analysis (PCA) in Python step by step), but in addition to finding the component axes that maximize the variance of our data (PCA), we are additionally interested in the axes that maximize the separation between multiple classes (LDA).
So, in a nutshell, often the goal of an LDA is to project a feature space (a dataset n-dimensional samples) onto a smaller subspace $k$ (where $k \leq n-1$) while maintaining the class-discriminatory information. 
In general, dimensionality reduction does not only help reducing computational costs for a given classification task, but it can also be helpful to avoid overfitting by minimizing the error in parameter estimation ("curse of dimensionality").
<br>
<br>
Principal Component Analysis vs. Linear Discriminant Analysis
[back to top]
Both Linear Discriminant Analysis (LDA) and Principal Component Analysis (PCA) are linear transformation techniques that are commonly used for dimensionality reduction. PCA can be described as an "unsupervised" algorithm, since it "ignores" class labels and its goal is to find the directions (the so-called principal components) that maximize the variance in a dataset.
In contrast to PCA, LDA is "supervised" and computes the directions ("linear discriminants") that  will represent the axes that that maximize the separation between multiple classes.
Although it might sound intuitive that LDA is superior to PCA for a multi-class classification task where the class labels are known, this might not always the case.
For example, comparisons between classification accuracies for image recognition after using PCA or LDA show that PCA tends to outperform LDA if the number of samples per class is relatively small (PCA vs. LDA, A.M. Martinez et al., 2001).
In practice, it is also not uncommon to use both LDA and PCA in combination: E.g., PCA for dimensionality reduction followed by an LDA.
<br>
<br>
<br>
<br>
What is a "good" feature subspace?
[back to top]
Let's assume that our goal is to reduce the dimensions of a $d$-dimensional dataset by projecting it onto a $(k)$-dimensional subspace (where $k\;<\;d$). 
So, how do we know what size we should choose for $k$ ($k$ = the number of dimensions of the new feature subspace), and how do we know if we have a feature space that represents our data "well"?  
Later, we will compute eigenvectors (the components) from our data set and collect them in a so-called scatter-matrices (i.e., the in-between-class scatter matrix and within-class scatter matrix).
Each of these eigenvectors is associated with an eigenvalue, which tells us about the "length" or "magnitude" of the eigenvectors. 
If we would observe that all eigenvalues have a similar magnitude, then this may be a good indicator that our data is already projected on a "good" feature space.  
And in the other scenario, if some of the eigenvalues are much much larger than others, we might be interested in keeping only those eigenvectors with the highest eigenvalues, since they contain more information about our data distribution. Vice versa, eigenvalues that are close to 0 are less informative and we might consider dropping those for constructing the new feature subspace.
<br>
<br>
Summarizing the LDA approach in 5 steps
[back to top]
Listed below are the 5 general steps for performing a linear discriminant analysis; we will explore them in more detail in the following sections.
Compute the $d$-dimensional mean vectors for the different classes from the dataset.
Compute the scatter matrices (in-between-class and within-class scatter matrix).
Compute the eigenvectors ($\pmb e_1, \; \pmb e_2, \; ..., \; \pmb e_d$) and corresponding eigenvalues ($\pmb \lambda_1, \; \pmb \lambda_2, \; ..., \; \pmb \lambda_d$) for the scatter matrices.
Sort the eigenvectors by decreasing eigenvalues and choose $k$ eigenvectors with the largest eigenvalues to form a $k \times d$ dimensional matrix $\pmb W\;$ (where every column represents an eigenvector).
Use this $k \times d$ eigenvector matrix to transform the samples onto the new subspace. This can be summarized by the mathematical equation: $\pmb Y = \pmb X \times \pmb W$ (where $\pmb X$ is a $n \times d$-dimensional matrix representing the $n$ samples, and $\pmb y$ are the transformed $n \times k$-dimensional samples in the new subspace).
<a name="sample_data"></a>
<br>
<br>
Preparing the sample data set
[back to top]
<a name="sample_data"></a>
<br>
<br>
About the Iris dataset
[back to top]
For the following tutorial, we will be working with the famous "Iris" dataset that has been deposited on the UCI machine learning repository
(https://archive.ics.uci.edu/ml/datasets/Iris).
<font size="1">
Reference:
Bache, K. & Lichman, M. (2013). UCI Machine Learning Repository. Irvine, CA: University of California, School of Information and Computer Science.</font>
The iris dataset contains measurements for 150 iris flowers from  three different species.
The three classes in the Iris dataset:
Iris-setosa (n=50)
Iris-versicolor (n=50)
Iris-virginica (n=50)
The four features of the Iris dataset:
sepal length in cm
sepal width in cm
petal length in cm
petal width in cm
End of explanation
import pandas as pd
df = pd.io.parsers.read_csv(
    filepath_or_buffer='https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', 
    header=None, 
    sep=',', 
    )
df.columns = [l for i,l in sorted(feature_dict.items())] + ['class label']
df.dropna(how="all", inplace=True) # to drop the empty line at file-end
df.tail()
Explanation: <a name="sample_data"></a>
<br>
<br>
Reading in the dataset
[back to top]
End of explanation
from sklearn.preprocessing import LabelEncoder
X = df[[0,1,2,3]].values 
y = df['class label'].values 
enc = LabelEncoder()
label_encoder = enc.fit(y)
y = label_encoder.transform(y) + 1
label_dict = {1: 'Setosa', 2: 'Versicolor', 3:'Virginica'}
Explanation: $\pmb X = \begin{bmatrix} x_{1_{\text{sepal length}}} & x_{1_{\text{sepal width}}} & x_{1_{\text{petal length}}} & x_{1_{\text{petal width}}}\ 
x_{2_{\text{sepal length}}} & x_{2_{\text{sepal width}}} & x_{2_{\text{petal length}}} & x_{2_{\text{petal width}}}\ 
...  \
x_{150_{\text{sepal length}}} & x_{150_{\text{sepal width}}} & x_{150_{\text{petal length}}} & x_{150_{\text{petal width}}}\ 
\end{bmatrix}, \;\;
\pmb y = \begin{bmatrix} \omega_{\text{setosa}}\ 
\omega_{\text{setosa}}\ 
...  \
\omega_{\text{virginica}}\end{bmatrix}$
<a name="sample_data"></a>
<br>
<br>
Since it is more convenient to work with numerical values, we will use the LabelEncode from the scikit-learn library to convert the class labels into numbers: 1, 2, and 3.
End of explanation
%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
import math
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12,6))
for ax,cnt in zip(axes.ravel(), range(4)):  
    
    # set bin sizes
    min_b = math.floor(np.min(X[:,cnt]))
    max_b = math.ceil(np.max(X[:,cnt]))
    bins = np.linspace(min_b, max_b, 25)
    
    # plottling the histograms
    for lab,col in zip(range(1,4), ('blue', 'red', 'green')):
        ax.hist(X[y==lab, cnt],
                   color=col, 
                   label='class %s' %label_dict[lab], 
                   bins=bins,
                   alpha=0.5,)
    ylims = ax.get_ylim()
    
    # plot annotation
    leg = ax.legend(loc='upper right', fancybox=True, fontsize=8)
    leg.get_frame().set_alpha(0.5)
    ax.set_ylim([0, max(ylims)+2])
    ax.set_xlabel(feature_dict[cnt])
    ax.set_title('Iris histogram #%s' %str(cnt+1))
    
    # hide axis ticks
    ax.tick_params(axis="both", which="both", bottom="off", top="off",  
            labelbottom="on", left="off", right="off", labelleft="on")
    # remove axis spines
    ax.spines["top"].set_visible(False)  
    ax.spines["right"].set_visible(False) 
    ax.spines["bottom"].set_visible(False) 
    ax.spines["left"].set_visible(False)    
 
axes[0][0].set_ylabel('count')
axes[1][0].set_ylabel('count')
    
fig.tight_layout()       
        
plt.show()
Explanation: $\pmb y = \begin{bmatrix}{\text{setosa}}\ 
{\text{setosa}}\ 
...  \
{\text{virginica}}\end{bmatrix} \quad \Rightarrow
\begin{bmatrix} {\text{1}}\ 
{\text{1}}\ 
...  \
{\text{3}}\end{bmatrix}$
<a name="sample_data"></a>
<br>
<br>
Histograms and feature selection
[back to top]
Just to get a rough idea how the samples of our three classes $\omega_1$, $\omega_2$ and $\omega_3$ are distributed, let us visualize the distributions of the four different features in 1-dimensional histograms.
End of explanation
from sklearn import preprocessing
preprocessing.scale(X, axis=0, with_mean=True, with_std=True, copy=False)
print()
Explanation: From just looking at these simple graphical representations of the features, we can already tell that the petal lengths and widths are likely better suited as potential features two separate between the three flower classes. In practice, instead of reducing the dimensionality via a projection (here: LDA), a good alternative would be a feature selection technique. For low-dimensional datasets like Iris, a glance at those histograms would already be very informative. Another simple, but very useful technique would be to use feature selection algorithms, which I have described in more detail in another article: Feature Selection Algorithms in Python
<a name="sample_data"></a>
<br>
<br>
Standardization
[back to top]
Normalization is one important part of every data pre-processing step and typically a requirement for best performances of many machine learning algorithms.
The two most popular approaches for data normalization are the so-called "standardization" and "min-max scaling".
Standardization (or Z-score normalization): Rescaling of the features so that they'll have the properties of a standard normal distribution with &mu;=0 and &sigma;=1 (i.e., unit variance centered around the mean).  
Min-max scaling: Rescaling of the features to unit range, typically a range between 0 and 1. Quite often, min-max scaling is also just called "normalization", which can be quite confusing depending on the context where the term is being used. Via Min-max scaling, 
Both are very important procedures, so that I have also a separate article about it with more details: About Feature Scaling and Normalization.
In our case, although the features are already on the same scale (measured in centimeters), we still want to scale the features to unit variance (&sigma;=1, &mu;=0).
End of explanation
np.set_printoptions(precision=4)
mean_vectors = []
for cl in range(1,4):
    mean_vectors.append(np.mean(X[y==cl], axis=0))
    print('Mean Vector class %s: %s\n' %(cl, mean_vectors[cl-1]))
Explanation: <a name="sample_data"></a>
<br>
<br>
Normality assumptions
[back to top]
It should be mentioned that LDA assumes normal distributed data, features that are statistically independent, and identical covariance matrices for every class. However, this only applies for LDA as classifier and LDA for dimensionality reduction can also work reasonably well if those assumptions are violated. And even for classification tasks LDA seems can be quite robust to the distribution of the data: 
"linear discriminant analysis frequently achieves good performances in
the tasks of face and object recognition, even though the assumptions
of common covariance matrix among groups and normality are often
violated (Duda, et al., 2001)" (Tao Li, et al., 2006).
<br>
<font size="1">References:  
Tao Li, Shenghuo Zhu, and Mitsunori Ogihara. “Using Discriminant Analysis for Multi-Class Classification: An Experimental Investigation.” Knowledge and Information Systems 10, no. 4 (2006): 453–72.)   
Duda, Richard O, Peter E Hart, and David G Stork. 2001. Pattern Classification. New York: Wiley.</font>
<a name="sample_data"></a>
<br>
<br>
LDA in 5 steps
[back to top]
After we went through several preparation steps, our data is finally ready for the actual LDA. In practice, LDA for dimensionality reduction would be just another preprocessing step for a typical machine learning or pattern classification task.
<a name="sample_data"></a>
<br>
<br>
Step 1: Computing the d-dimensional mean vectors
[back to top]
In this first step, we will start off with a simple computation of the mean vectors $\pmb m_i$, $(i = 1,2,3)$ of the 3 different flower classes:
$\pmb m_i = \begin{bmatrix} 
\mu_{\omega_i (\text{sepal length)}}\ 
\mu_{\omega_i (\text{sepal width})}\ 
\mu_{\omega_i (\text{petal length)}}\
\mu_{\omega_i (\text{petal width})}\
\end{bmatrix} \; , \quad \text{with} \quad i = 1,2,3$
End of explanation
S_W = np.zeros((4,4))
for cl,mv in zip(range(1,4), mean_vectors):
    class_sc_mat = np.zeros((4,4))                  # scatter matrix for every class
    for row in X[y == cl]:
        row, mv = row.reshape(4,1), mv.reshape(4,1) # make column vectors
        class_sc_mat += (row-mv).dot((row-mv).T)
    S_W += class_sc_mat                             # sum class scatter matrices
print('within-class Scatter Matrix:\n', S_W)
Explanation: <a name="sample_data"></a>
<br>
<br>
<a name="sc_matrix"></a>
Step 2: Computing the Scatter Matrices
[back to top]
Now, we will compute the two 4x4-dimensional matrices: The within-class and the between-class scatter matrix.
<br>
<br>
2.1 Within-class scatter matrix $S_W$
[back to top]
The within-class scatter matrix $S_W$ is computed by the following equation:  
$S_W = \sum\limits_{i=1}^{c}  S_i$
where
$S_i = \sum\limits_{\pmb x \in D_i}^n (\pmb x - \pmb m_i)\;(\pmb x - \pmb m_i)^T$
(scatter matrix for every class) 
and $\pmb m_i$ is the mean vector  
$\pmb m_i = \frac{1}{n_i} \sum\limits_{\pmb x \in D_i}^n \; \pmb x_k$
End of explanation
overall_mean = np.mean(X, axis=0)
S_B = np.zeros((4,4))
for i,mean_vec in enumerate(mean_vectors):  
    n = X[y==i+1,:].shape[0]
    mean_vec = mean_vec.reshape(4,1) # make column vector
    overall_mean = overall_mean.reshape(4,1) # make column vector
    S_B += n * (mean_vec - overall_mean).dot((mean_vec - overall_mean).T)
    
print('between-class Scatter Matrix:\n', S_B)
Explanation: <br>
2.1 b
Alternatively, we could also compute the class-covariance matrices by adding the scaling factor $\frac{1}{N-1}$ to the within-class scatter matrix, so that our equation becomes
$\Sigma_i = \frac{1}{N_{i}-1} \sum\limits_{\pmb x \in D_i}^n (\pmb x - \pmb m_i)\;(\pmb x - \pmb m_i)^T$.
and $S_W = \sum\limits_{i=1}^{c} (N_{i}-1) \Sigma_i$
where $N_{i}$ is the sample size of the respective class (here: 50), and in this particular case, we can drop the term ($N_{i}-1)$ 
since all classes have the same sample size.
However, the resulting eigenspaces will be identical (identical eigenvectors, only the eigenvalues are scaled differently by a constant factor).
<br>
<br>
2.2 Between-class scatter matrix $S_B$
[back to top]
The between-class scatter matrix $S_B$ is computed by the following equation:  
$S_B =  \sum\limits_{i=1}^{c} N_{i} (\pmb m_i - \pmb m) (\pmb m_i - \pmb m)^T$
where
 $\pmb m$ is the overall mean, and $\pmb m_{i}$ and $N_{i}$ are the sample mean and sizes of the respective classes.
End of explanation
eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_B))
for i in range(len(eig_vals)):
    eigvec_sc = eig_vecs[:,i].reshape(4,1)   
    print('\nEigenvector {}: \n{}'.format(i+1, eigvec_sc.real))
    print('Eigenvalue {:}: {:.2e}'.format(i+1, eig_vals[i].real))
Explanation: <br>
<br>
Step 3: Solving the generalized eigenvalue problem for the matrix $S_{W}^{-1}S_B$
[back to top]
Next, we will solve the generalized eigenvalue problem for the matrix $S_{W}^{-1}S_B$ to obtain the linear discriminants.
End of explanation
for i in range(len(eig_vals)):
    eigv = eig_vecs[:,i].reshape(4,1) 
    np.testing.assert_array_almost_equal(np.linalg.inv(S_W).dot(S_B).dot(eigv), 
                                         eig_vals[i] * eigv, 
                                         decimal=6, err_msg='', verbose=True)
print('ok')
Explanation: <br>
<br>
After this decomposition of our square matrix into eigenvectors and eigenvalues, let us briefly recapitulate how we can interpret those results. As we remember from our first linear algebra class in high school or college, both eigenvectors and eigenvalues are providing us with information about the distortion of a linear transformation: The eigenvectors are basically the direction of this distortion, and the eigenvalues are the scaling factor for the eigenvectors that describing the magnitude of the distortion.  
If we are performing the LDA for dimensionality reduction, the eigenvectors are important since they will form the new axes of our new feature subspace; the associated eigenvalues are of particular interest since they will tell us how "informative" the new "axes" are.  
Let us briefly double-check our calculation and talk more about the eigenvalues in the next section.
<br>
<br>
Checking the eigenvector-eigenvalue calculation
[back to top]
A quick check that the eigenvector-eigenvalue calculation is correct and satisfy the equation:
$\pmb A\pmb{v} =  \lambda\pmb{v}$  
<br>
where
$\pmb A = S_{W}^{-1}S_B\
\pmb{v} = \; \text{Eigenvector}\
\lambda = \; \text{Eigenvalue}$
End of explanation
# Make a list of (eigenvalue, eigenvector) tuples
eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]
# Sort the (eigenvalue, eigenvector) tuples from high to low
eig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True)
# Visually confirm that the list is correctly sorted by decreasing eigenvalues
print('Eigenvalues in decreasing order:\n')
for i in eig_pairs:
    print(i[0])
Explanation: <br>
<br>
Step 4: Selecting linear discriminants for the new feature subspace
[back to top]
<br>
<br>
4.1. Sorting the eigenvectors by decreasing eigenvalues
[back to top]
Remember from the introduction that we are not only interested in merely projecting the data into a subspace that improves the class separability, but also reduces the dimensionality of our feature space, (where the eigenvectors will form the axes of this new feature subspace). 
However, the eigenvectors only define the directions of the new axis, since they have all the same unit length 1.
So, in order to decide which eigenvector(s) we want to drop for our lower-dimensional subspace, we have to take a look at the corresponding eigenvalues of the eigenvectors. Roughly speaking, the eigenvectors with the lowest eigenvalues bear the least information about the distribution of the data, and those are the ones we want to drop.
The common approach is to rank the eigenvectors from highest to lowest corresponding eigenvalue and choose the top $k$ eigenvectors.
End of explanation
print('Variance explained:\n')
eigv_sum = sum(eig_vals)
for i,j in enumerate(eig_pairs):
    print('eigenvalue {0:}: {1:.2%}'.format(i+1, (j[0]/eigv_sum).real))
Explanation: <br>
<br>
If we take a look at the eigenvalues, we can already see that 2 eigenvalues are close to 0 and conclude that the eigenpairs are less informative than the other two. Let's express the "explained variance" as percentage:
End of explanation
W = np.hstack((eig_pairs[0][1].reshape(4,1), eig_pairs[1][1].reshape(4,1)))
print('Matrix W:\n', W.real)
Explanation: <br>
<br>
The first eigenpair is by far the most informative one, and we won't loose much information if we would form a 1D-feature spaced based on this eigenpair.
<br>
<br>
4.2. Choosing k eigenvectors with the largest eigenvalues
[back to top]
After sorting the eigenpairs by decreasing eigenvalues, it is now time to construct our $k \times d$-dimensional eigenvector matrix $\pmb W$ (here $4 \times 2$: based on the 2 most informative eigenpairs) and thereby reducing the initial 4-dimensional feature space into a 2-dimensional feature subspace.
End of explanation
X_lda = X.dot(W)
assert X_lda.shape == (150,2), "The matrix is not 2x150 dimensional."
from matplotlib import pyplot as plt
def plot_step_lda():
    
    ax = plt.subplot(111)
    for label,marker,color in zip(
        range(1,4),('^', 's', 'o'),('blue', 'red', 'green')):
        plt.scatter(x=X_lda[:,0].real[y == label],
                y=X_lda[:,1].real[y == label],
                marker=marker,
                color=color,
                alpha=0.5,
                label=label_dict[label]
                )
    plt.xlabel('LD1')
    plt.ylabel('LD2')
    leg = plt.legend(loc='upper right', fancybox=True)
    leg.get_frame().set_alpha(0.5)
    plt.title('LDA: Iris projection onto the first 2 linear discriminants')
    
    # hide axis ticks
    plt.tick_params(axis="both", which="both", bottom="off", top="off",  
            labelbottom="on", left="off", right="off", labelleft="on")
    # remove axis spines
    ax.spines["top"].set_visible(False)  
    ax.spines["right"].set_visible(False) 
    ax.spines["bottom"].set_visible(False) 
    ax.spines["left"].set_visible(False)    
 
    plt.grid()
    plt.tight_layout
    plt.show()
    
plot_step_lda()
Explanation: <br>
<br>
Step 5: Transforming the samples onto the new subspace
[back to top]
In the last step, we use the $4 \times 2$-dimensional matrix $\pmb W$ that we just computed to transform our samples onto the new subspace via the equation   
$\pmb Y = \pmb X \times \pmb W $.
(where $\pmb X$ is a $n \times d$-dimensional matrix representing the $n$ samples, and $\pmb Y$ are the transformed $n \times k$-dimensional samples in the new subspace).
End of explanation
from sklearn.decomposition import PCA as sklearnPCA
sklearn_pca = sklearnPCA(n_components=2)
X_pca = sklearn_pca.fit_transform(X)
def plot_pca():
    ax = plt.subplot(111)
    
    for label,marker,color in zip(
        range(1,4),('^', 's', 'o'),('blue', 'red', 'green')):
        plt.scatter(x=X_pca[:,0][y == label],
                y=X_pca[:,1][y == label],
                marker=marker,
                color=color,
                alpha=0.5,
                label=label_dict[label]
                )
    plt.xlabel('PC1')
    plt.ylabel('PC2')
    leg = plt.legend(loc='upper right', fancybox=True)
    leg.get_frame().set_alpha(0.5)
    plt.title('PCA: Iris projection onto the first 2 principal components')
    # hide axis ticks
    plt.tick_params(axis="both", which="both", bottom="off", top="off",  
            labelbottom="on", left="off", right="off", labelleft="on")
    # remove axis spines
    ax.spines["top"].set_visible(False)  
    ax.spines["right"].set_visible(False) 
    ax.spines["bottom"].set_visible(False) 
    ax.spines["left"].set_visible(False)    
    
    plt.tight_layout
    plt.grid()
    
    plt.show()
plot_pca()
plot_step_lda()
Explanation: The scatter plot above represents our new feature subspace that we constructed via LDA. We can see that the first linear discriminant "LD1" separates the classes quite nicely. However, the second discriminant, "LD2", does not add much valuable information, which we've already concluded when we looked at the ranked eigenvalues is step 4.
<br>
<br>
A comparison of PCA and LDA
[back to top]
In order to compare the feature subspace that we obtained via the Linear Discriminant Analysis, we will use the PCA class from the scikit-learn machine-learning library. The documentation can be found here:
http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html.  
For our convenience, we can directly specify to how many components we want to retain in our input dataset via the n_components parameter. 
n_components : int, None or string
Number of components to keep. if n_components is not set all components are kept:
    n_components == min(n_samples, n_features)
    if n_components == ‘mle’, Minka’s MLE is used to guess the dimension if 0 &lt; n_components &lt; 1, 
    select the number of components such that the amount of variance that needs to be explained 
    is greater than the percentage specified by n_components
<br>
<br>
But before we skip to the results of the respective linear transformations, let us quickly recapitulate the purposes of PCA and LDA: PCA finds the axes with maximum variance for the whole data set where LDA tries to find the axes for best class seperability. In practice, often a LDA is done followed by a PCA for dimensionality reduction.
<br>
<br>
<br>
<br>
End of explanation
from sklearn.lda import LDA
# LDA
sklearn_lda = LDA(n_components=2)
X_lda_sklearn = sklearn_lda.fit_transform(X, y)
def plot_scikit_lda(X, title, mirror=1):
    
    ax = plt.subplot(111)
    for label,marker,color in zip(
        range(1,4),('^', 's', 'o'),('blue', 'red', 'green')):
        
        plt.scatter(x=X[:,0][y == label]*mirror,
                y=X[:,1][y == label],
                marker=marker,
                color=color,
                alpha=0.5,
                label=label_dict[label]
                )
    plt.xlabel('LD1')
    plt.ylabel('LD2')
    leg = plt.legend(loc='upper right', fancybox=True)
    leg.get_frame().set_alpha(0.5)
    plt.title(title)
    
    # hide axis ticks
    plt.tick_params(axis="both", which="both", bottom="off", top="off",  
            labelbottom="on", left="off", right="off", labelleft="on")
    # remove axis spines
    ax.spines["top"].set_visible(False)  
    ax.spines["right"].set_visible(False) 
    ax.spines["bottom"].set_visible(False) 
    ax.spines["left"].set_visible(False)    
 
    plt.grid()
    plt.tight_layout
    plt.show()
plot_step_lda()
plot_scikit_lda(X_lda_sklearn, title='Default LDA via scikit-learn', mirror=(-1))
Explanation: <br>
<br>
The two plots above nicely confirm what we have discussed before: Where the PCA accounts for the most variance in the whole dataset, the LDA gives us the axes that account for the most variance between the individual classes.
<br>
<br>
LDA via scikit-learn
[back to top]
Now, after we have seen how an Linear Discriminant Analysis works using a step-by-step approach, there is also a more convenient way to achive the same via the LDA class implemented in the scikit-learn machine learning library.
End of explanation
<END_TASK>