mhjensen
844 lines
29 KiB
Plaintext
844 lines
29 KiB
Plaintext
TITLE: Data Analysis and Machine Learning: Preprocessing and Dimensionality Reduction
|
||
AUTHOR: Morten Hjorth-Jensen {copyright, 1999-present|CC BY-NC} at Department of Physics, University of Oslo & Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University
|
||
DATE: today
|
||
|
||
|
||
!split
|
||
===== Reducing the number of degrees of freedom, overarching view =====
|
||
!bblock
|
||
|
||
Many Machine Learning problems involve thousands or even millions of
|
||
features for each training instance. Not only does this make training
|
||
extremely slow, it can also make it much harder to find a good
|
||
solution, as we will see. This problem is often referred to as the
|
||
curse of dimensionality. Fortunately, in real-world problems, it is
|
||
often possible to reduce the number of features considerably, turning
|
||
an intractable problem into a tractable one.
|
||
|
||
Here we will discuss some of the most popular dimensionality reduction
|
||
techniques: the principal component analysis PCA, Kernel PCA, and
|
||
Locally Linear Embedding (LLE). Furthermore, we will start by looking
|
||
at some simple preprocessing of the data which allow us to rescale the
|
||
data.
|
||
|
||
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Preprocessing our data =====
|
||
!bblock
|
||
|
||
Before we proceed however, we will discuss how to preprocess our
|
||
data. Till now and in connection with our previous examples we have
|
||
not met so many cases where we are too sensitive to the scaling of our
|
||
data. Normally the data may need a rescaling and/or may be sensitive
|
||
to extreme values. Scaling the data renders our inputs much more
|
||
suitable for the algorithms we want to employ.
|
||
|
||
_Scikit-Learn_ has several functions which allow us to rescale the
|
||
data, normally resulting in much better results in terms of various
|
||
accuracy scores. The _StandardScaler_ function in _Scikit-Learn_
|
||
ensures that for each feature/predictor we study the mean value is
|
||
zero and the variance is one (every column in the design/feature
|
||
matrix). This scaling has the drawback that it does not ensure that
|
||
we have a particular maximum or minimum in our data set. Another
|
||
function included in _Scikit-Learn_ is the _MinMaxScaler_ which
|
||
ensures that all features are exactly between $0$ and $1$. The
|
||
|
||
!split
|
||
===== More preprocessing =====
|
||
|
||
|
||
The _Normalizer_ scales each data
|
||
point such that the feature vector has a euclidean length of one. In other words, it
|
||
projects a data point on the circle (or sphere in the case of higher dimensions) with a
|
||
radius of 1. This means every data point is scaled by a different number (by the
|
||
inverse of it’s length).
|
||
This normalization is often used when only the direction (or angle) of the data matters,
|
||
not the length of the feature vector.
|
||
|
||
The _RobustScaler_ works similarly to the StandardScaler in that it
|
||
ensures statistical properties for each feature that guarantee that
|
||
they are on the same scale. However, the RobustScaler uses the median
|
||
and quartiles, instead of mean and variance. This makes the
|
||
RobustScaler ignore data points that are very different from the rest
|
||
(like measurement errors). These odd data points are also called
|
||
outliers, and might often lead to trouble for other scaling
|
||
techniques.
|
||
|
||
!eblock
|
||
|
||
!split
|
||
===== Simple preprocessing examples, Franke function and regression =====
|
||
|
||
!bc pycod
|
||
# Common imports
|
||
import os
|
||
import numpy as np
|
||
import pandas as pd
|
||
import matplotlib.pyplot as plt
|
||
import sklearn.linear_model as skl
|
||
from sklearn.metrics import mean_squared_error
|
||
from sklearn.model_selection import train_test_split
|
||
from sklearn.preprocessing import MinMaxScaler, StandardScaler, Normalizer
|
||
from sklearn.svm import SVR
|
||
|
||
# Where to save the figures and data files
|
||
PROJECT_ROOT_DIR = "Results"
|
||
FIGURE_ID = "Results/FigureFiles"
|
||
DATA_ID = "DataFiles/"
|
||
|
||
if not os.path.exists(PROJECT_ROOT_DIR):
|
||
os.mkdir(PROJECT_ROOT_DIR)
|
||
|
||
if not os.path.exists(FIGURE_ID):
|
||
os.makedirs(FIGURE_ID)
|
||
|
||
if not os.path.exists(DATA_ID):
|
||
os.makedirs(DATA_ID)
|
||
|
||
def image_path(fig_id):
|
||
return os.path.join(FIGURE_ID, fig_id)
|
||
|
||
def data_path(dat_id):
|
||
return os.path.join(DATA_ID, dat_id)
|
||
|
||
def save_fig(fig_id):
|
||
plt.savefig(image_path(fig_id) + ".png", format='png')
|
||
|
||
|
||
def FrankeFunction(x,y):
|
||
term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))
|
||
term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))
|
||
term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))
|
||
term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)
|
||
return term1 + term2 + term3 + term4
|
||
|
||
|
||
def create_X(x, y, n ):
|
||
if len(x.shape) > 1:
|
||
x = np.ravel(x)
|
||
y = np.ravel(y)
|
||
|
||
N = len(x)
|
||
l = int((n+1)*(n+2)/2) # Number of elements in beta
|
||
X = np.ones((N,l))
|
||
|
||
for i in range(1,n+1):
|
||
q = int((i)*(i+1)/2)
|
||
for k in range(i+1):
|
||
X[:,q+k] = (x**(i-k))*(y**k)
|
||
|
||
return X
|
||
|
||
|
||
# Making meshgrid of datapoints and compute Franke's function
|
||
n = 5
|
||
N = 1000
|
||
x = np.sort(np.random.uniform(0, 1, N))
|
||
y = np.sort(np.random.uniform(0, 1, N))
|
||
z = FrankeFunction(x, y)
|
||
X = create_X(x, y, n=n)
|
||
# split in training and test data
|
||
X_train, X_test, y_train, y_test = train_test_split(X,z,test_size=0.2)
|
||
|
||
|
||
svm = SVR(gamma='auto',C=10.0)
|
||
svm.fit(X_train, y_train)
|
||
|
||
# The mean squared error and R2 score
|
||
print("MSE before scaling: {:.2f}".format(mean_squared_error(svm.predict(X_test), y_test)))
|
||
print("R2 score before scaling {:.2f}".format(svm.score(X_test,y_test)))
|
||
|
||
scaler = StandardScaler()
|
||
scaler.fit(X_train)
|
||
X_train_scaled = scaler.transform(X_train)
|
||
X_test_scaled = scaler.transform(X_test)
|
||
|
||
print("Feature min values before scaling:\n {}".format(X_train.min(axis=0)))
|
||
print("Feature max values before scaling:\n {}".format(X_train.max(axis=0)))
|
||
|
||
print("Feature min values after scaling:\n {}".format(X_train_scaled.min(axis=0)))
|
||
print("Feature max values after scaling:\n {}".format(X_train_scaled.max(axis=0)))
|
||
|
||
svm = SVR(gamma='auto',C=10.0)
|
||
svm.fit(X_train_scaled, y_train)
|
||
|
||
print("MSE after scaling: {:.2f}".format(mean_squared_error(svm.predict(X_test_scaled), y_test)))
|
||
print("R2 score for scaled data: {:.2f}".format(svm.score(X_test_scaled,y_test)))
|
||
|
||
!ec
|
||
|
||
|
||
|
||
!split
|
||
===== Simple preprocessing examples, breast cancer data and classification, Support Vector Machines =====
|
||
|
||
We show here how we can use a simple regression case on the breast
|
||
cancer data using support vector machines (SVM) as algorithm for
|
||
classification.
|
||
|
||
|
||
!bc pycod
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
from sklearn.model_selection import train_test_split
|
||
from sklearn.datasets import load_breast_cancer
|
||
from sklearn.svm import SVC
|
||
cancer = load_breast_cancer()
|
||
|
||
X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
|
||
print(X_train.shape)
|
||
print(X_test.shape)
|
||
|
||
svm = SVC(C=100)
|
||
svm.fit(X_train, y_train)
|
||
print("Test set accuracy: {:.2f}".format(svm.score(X_test,y_test)))
|
||
|
||
from sklearn.preprocessing import MinMaxScaler, StandardScaler
|
||
scaler = MinMaxScaler()
|
||
scaler.fit(X_train)
|
||
X_train_scaled = scaler.transform(X_train)
|
||
X_test_scaled = scaler.transform(X_test)
|
||
|
||
print("Feature min values before scaling:\n {}".format(X_train.min(axis=0)))
|
||
print("Feature max values before scaling:\n {}".format(X_train.max(axis=0)))
|
||
|
||
print("Feature min values before scaling:\n {}".format(X_train_scaled.min(axis=0)))
|
||
print("Feature max values before scaling:\n {}".format(X_train_scaled.max(axis=0)))
|
||
|
||
|
||
svm.fit(X_train_scaled, y_train)
|
||
print("Test set accuracy scaled data with Min-Max scaling: {:.2f}".format(svm.score(X_test_scaled,y_test)))
|
||
|
||
scaler = StandardScaler()
|
||
scaler.fit(X_train)
|
||
X_train_scaled = scaler.transform(X_train)
|
||
X_test_scaled = scaler.transform(X_test)
|
||
|
||
svm.fit(X_train_scaled, y_train)
|
||
print("Test set accuracy scaled data with Standar Scaler: {:.2f}".format(svm.score(X_test_scaled,y_test)))
|
||
|
||
!ec
|
||
|
||
!split
|
||
===== More on Cancer Data, now with Logistic Regression =====
|
||
|
||
|
||
!bc pycod
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
from sklearn.model_selection import train_test_split
|
||
from sklearn.datasets import load_breast_cancer
|
||
from sklearn.linear_model import LogisticRegression
|
||
cancer = load_breast_cancer()
|
||
|
||
# Set up training data
|
||
X_train, X_test, y_train, y_test = train_test_split(cancer.data,cancer.target,random_state=0)
|
||
logreg = LogisticRegression()
|
||
logreg.fit(X_train, y_train)
|
||
print("Test set accuracy: {:.2f}".format(logreg.score(X_test,y_test)))
|
||
|
||
# Scale data
|
||
from sklearn.preprocessing import StandardScaler
|
||
scaler = StandardScaler()
|
||
scaler.fit(X_train)
|
||
X_train_scaled = scaler.transform(X_train)
|
||
X_test_scaled = scaler.transform(X_test)
|
||
logreg.fit(X_train_scaled, y_train)
|
||
print("Test set accuracy scaled data: {:.2f}".format(logreg.score(X_test_scaled,y_test)))
|
||
|
||
!ec
|
||
|
||
|
||
|
||
|
||
!split
|
||
===== Why should we think of reducing the dimensionality =====
|
||
|
||
In addition to the plot of the features, we study now also the covariance (or rather the correlation matrix).
|
||
We use also _Pandas_ to compute the correlation matrix.
|
||
!bc pycod
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
from sklearn.model_selection import train_test_split
|
||
from sklearn.datasets import load_breast_cancer
|
||
from sklearn.linear_model import LogisticRegression
|
||
cancer = load_breast_cancer()
|
||
import pandas as pd
|
||
# Making a data frame
|
||
cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names)
|
||
|
||
fig, axes = plt.subplots(15,2,figsize=(10,20))
|
||
malignant = cancer.data[cancer.target == 0]
|
||
benign = cancer.data[cancer.target == 1]
|
||
ax = axes.ravel()
|
||
|
||
for i in range(30):
|
||
_, bins = np.histogram(cancer.data[:,i], bins =50)
|
||
ax[i].hist(malignant[:,i], bins = bins, alpha = 0.5)
|
||
ax[i].hist(benign[:,i], bins = bins, alpha = 0.5)
|
||
ax[i].set_title(cancer.feature_names[i])
|
||
ax[i].set_yticks(())
|
||
ax[0].set_xlabel("Feature magnitude")
|
||
ax[0].set_ylabel("Frequency")
|
||
ax[0].legend(["Malignant", "Benign"], loc ="best")
|
||
fig.tight_layout()
|
||
plt.show()
|
||
|
||
import seaborn as sns
|
||
correlation_matrix = cancerpd.corr().round(1)
|
||
# use the heatmap function from seaborn to plot the correlation matrix
|
||
# annot = True to print the values inside the square
|
||
sns.heatmap(data=correlation_matrix, annot=True)
|
||
plt.show()
|
||
|
||
#print eigvalues of correlation matrix
|
||
EigValues, EigVectors = np.linalg.eig(correlation_matrix)
|
||
print(EigValues)
|
||
!ec
|
||
|
||
In the above example we note two things. In the first plot we display
|
||
the overlap of benign and malignant tumors as functions of the various
|
||
features in the Wisconsing breast cancer data set. We see that for
|
||
some of the features we can distinguish clearly the benign and
|
||
malignant cases while for other features we cannot. This can point to
|
||
us which features may be of greater interest when we wish to classify
|
||
a benign or not benign tumour.
|
||
|
||
In the second figure we have computed the so-called correlation
|
||
matrix, which in our case with thirty features becomes a $30\times 30$
|
||
matrix.
|
||
|
||
We constructed this matrix using _pandas_ via the statements
|
||
!bc pycod
|
||
cancerpd = pd.DataFrame(cancer.data, columns=cancer.feature_names)
|
||
!ec
|
||
and then
|
||
!bc pycod
|
||
correlation_matrix = cancerpd.corr().round(1)
|
||
!ec
|
||
|
||
Diagonalizing this matrix we can in turn say something about which
|
||
features are of relevance and which are not. But before we proceed we
|
||
need to define covariance and correlation matrices. This leads us to
|
||
the classical Principal Component Analysis (PCA) theorem with
|
||
applications.
|
||
|
||
|
||
|
||
!split
|
||
===== Basic ideas of the Principal Component Analysis (PCA) =====
|
||
|
||
We have a data set defined by a design/feature matrix $\bm{X}$ (see below for its definition)
|
||
* Each data point is determined by $p$ extrinsic (measurement) variables
|
||
* We may want to ask the following question: Are there fewer intrinsic variables (say $d << p$) that still approximately describe the data?
|
||
* If so, these intrinsic variables may tell us something important and finding these intrinsic variables is what dimension reduction methods do.
|
||
|
||
|
||
!split
|
||
===== Introducing the Covariance and Correlation functions =====
|
||
|
||
Before we discuss the PCA theorem, we need to remind ourselves about the definition of the covariance and the correlation function.
|
||
|
||
Suppose we have defined two vectors
|
||
$\hat{x}$ and $\hat{y}$ with $n$ elements each. The covariance matrix $\bm{C}$ is defined as
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x},\bm{y}] = \begin{bmatrix} \mathrm{cov}[\bm{x},\bm{x}] & \mathrm{cov}[\bm{x},\bm{y}] \\
|
||
\mathrm{cov}[\bm{y},\bm{x}] & \mathrm{cov}[\bm{y},\bm{y}] \\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
where for example
|
||
!bt
|
||
\[
|
||
\mathrm{cov}[\bm{x},\bm{y}] =\frac{1}{n} \sum_{i=0}^{n-1}(x_i- \overline{x})(y_i- \overline{y}).
|
||
\]
|
||
!et
|
||
With this definition and recalling that the variance is defined as
|
||
!bt
|
||
\[
|
||
\mathrm{var}[\bm{x}]=\frac{1}{n} \sum_{i=0}^{n-1}(x_i- \overline{x})^2,
|
||
\]
|
||
!et
|
||
we can rewrite the covariance matrix as
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x},\bm{y}] = \begin{bmatrix} \mathrm{var}[\bm{x}] & \mathrm{cov}[\bm{x},\bm{y}] \\
|
||
\mathrm{cov}[\bm{x},\bm{y}] & \mathrm{var}[\bm{y}] \\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
|
||
The covariance takes values between zero and infinity and may thus
|
||
lead to problems with loss of numerical precision for particularly
|
||
large values. It is common to scale the covariance matrix by
|
||
introducing instead the correlation matrix defined via the so-called
|
||
correlation function
|
||
|
||
!bt
|
||
\[
|
||
\mathrm{corr}[\bm{x},\bm{y}]=\frac{\mathrm{cov}[\bm{x},\bm{y}]}{\sqrt{\mathrm{var}[\bm{x}] \mathrm{var}[\bm{y}]}}.
|
||
\]
|
||
!et
|
||
|
||
The correlation function is then given by values $\mathrm{corr}[\bm{x},\bm{y}]
|
||
\in [-1,1]$. This avoids eventual problems with too large values. We
|
||
can then define the correlation matrix for the two vectors $\bm{x}$
|
||
and $\bm{y}$ as
|
||
|
||
!bt
|
||
\[
|
||
\bm{K}[\bm{x},\bm{y}] = \begin{bmatrix} 1 & \mathrm{corr}[\bm{x},\bm{y}] \\
|
||
\mathrm{corr}[\bm{y},\bm{x}] & 1 \\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
|
||
In the above example this is the function we constructed using _pandas_.
|
||
|
||
!split
|
||
===== Correlation Function and Design/Feature Matrix =====
|
||
|
||
In our derivation of the various regression algorithms like Ordinary Least Squares or Ridge regression we defined the design/feature matrix $\bm{X}$ as
|
||
!bt
|
||
\[
|
||
\bm{X}=\begin{bmatrix}
|
||
x_{0,0} & x_{0,1} & x_{0,2}& \dots & \dots x_{0,p-1}\\
|
||
x_{1,0} & x_{1,1} & x_{1,2}& \dots & \dots x_{1,p-1}\\
|
||
x_{2,0} & x_{2,1} & x_{2,2}& \dots & \dots x_{2,p-1}\\
|
||
\dots & \dots & \dots & \dots \dots & \dots \\
|
||
x_{n-2,0} & x_{n-2,1} & x_{n-2,2}& \dots & \dots x_{n-2,p-1}\\
|
||
x_{n-1,0} & x_{n-1,1} & x_{n-1,2}& \dots & \dots x_{n-1,p-1}\\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
with $\bm{X}\in {\mathbb{R}}^{n\times p}$, with the predictors/features $p$ refering to the column numbers and the
|
||
entries $n$ being the row elements.
|
||
We can rewrite the design/feature matrix in terms of its column vectors as
|
||
!bt
|
||
\[
|
||
\bm{X}=\begin{bmatrix} \bm{x}_0 & \bm{x}_0 & \bm{x}_0 & \dots & \dots & \bm{x}_{p-1}\end{bmatrix},
|
||
\]
|
||
!et
|
||
with a given vector
|
||
!bt
|
||
\[
|
||
\bm{x}_i^T = \begin{bmatrix}x_{0,i} & x_{1,i} & x_{2,i}& \dots & \dots x_{n-1,i}\end{bmatrix}.
|
||
\]
|
||
!et
|
||
|
||
With these definitions, we can now rewrite our $2\times 2$ correaltion/covariance matrix in terms of a moe general design/feature matrix $\bm{X}\in {\mathbb{R}}^{n\times p}$. This leads to a $p\times p$ covariance matrix for the vectors $\bm{x}_i$ with $i =0,1,\dots,p-1$
|
||
!bt
|
||
\[
|
||
\bm{C}[\bm{x}] = \begin{bmatrix}
|
||
\mathrm{var}[\bm{x}_0] & \mathrm{cov}[\bm{x}_0,\bm{x}_1] & \mathrm{cov}[\bm{x}_0,\bm{x}_2] & \dots & \dots & \mathrm{cov}[\bm{x}_0,\bm{x}_{p-1}]\\
|
||
\mathrm{cov}[\bm{x}_1,\bm{x}_0] & \mathrm{var}[\bm{x}_1] & \mathrm{cov}[\bm{x}_1,\bm{x}_2] & \dots & \dots & \mathrm{cov}[\bm{x}_1,\bm{x}_{p-1}]\\
|
||
\mathrm{cov}[\bm{x}_2,\bm{x}_0] & \mathrm{cov}[\bm{x}_2,\bm{x}_1] & \mathrm{var}[\bm{x}_2] & \dots & \dots & \mathrm{cov}[\bm{x}_2,\bm{x}_{p-1}]\\
|
||
\dots & \dots & \dots & \dots & \dots & \dots \\
|
||
\dots & \dots & \dots & \dots & \dots & \dots \\
|
||
\mathrm{cov}[\bm{x}_{p-1},\bm{x}_0] & \mathrm{cov}[\bm{x}_{p-1},\bm{x}_1] & \mathrm{cov}[\bm{x}_{p-1},\bm{x}_{2}] & \dots & \dots & \mathrm{var}[\bm{x}_{p-1}]\\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
and the correlation matrix
|
||
!bt
|
||
\[
|
||
\bm{K}[\bm{x}] = \begin{bmatrix}
|
||
1 & \mathrm{corr}[\bm{x}_0,\bm{x}_1] & \mathrm{corr}[\bm{x}_0,\bm{x}_2] & \dots & \dots & \mathrm{corr}[\bm{x}_0,\bm{x}_{p-1}]\\
|
||
\mathrm{corr}[\bm{x}_1,\bm{x}_0] & 1 & \mathrm{corr}[\bm{x}_1,\bm{x}_2] & \dots & \dots & \mathrm{corr}[\bm{x}_1,\bm{x}_{p-1}]\\
|
||
\mathrm{corr}[\bm{x}_2,\bm{x}_0] & \mathrm{corr}[\bm{x}_2,\bm{x}_1] & 1 & \dots & \dots & \mathrm{corr}[\bm{x}_2,\bm{x}_{p-1}]\\
|
||
\dots & \dots & \dots & \dots & \dots & \dots \\
|
||
\dots & \dots & \dots & \dots & \dots & \dots \\
|
||
\mathrm{corr}[\bm{x}_{p-1},\bm{x}_0] & \mathrm{corr}[\bm{x}_{p-1},\bm{x}_1] & \mathrm{corr}[\bm{x}_{p-1},\bm{x}_{2}] & \dots & \dots & 1\\
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
|
||
|
||
!split
|
||
===== Covariance Matrix Examples =====
|
||
|
||
|
||
The Numpy function _np.cov_ calculates the covariance elements using
|
||
the factor $1/(n-1)$ instead of $1/n$ since it assumes we do not have
|
||
the exact mean values. The following simple function uses the
|
||
_np.vstack_ function which takes each vector of dimension $1\times n$
|
||
and produces a $2\times n$ matrix $\bm{W}$
|
||
|
||
|
||
!bt
|
||
\[
|
||
\bm{W} = \begin{bmatrix} x_0 & y_0 \\
|
||
x_1 & y_1 \\
|
||
x_2 & y_2\\
|
||
\dots & \dots \\
|
||
x_{n-2} & y_{n-2}\\
|
||
x_{n-1} & y_{n-1} &
|
||
\end{bmatrix},
|
||
\]
|
||
!et
|
||
|
||
which in turn is converted into into the $2\times 2$ covariance matrix
|
||
$\bm{C}$ via the Numpy function _np.cov()_. We note that we can also calculate
|
||
the mean value of each set of samples $\bm{x}$ etc using the Numpy
|
||
function _np.mean(x)_. We can also extract the eigenvalues of the
|
||
covariance matrix through the _np.linalg.eig()_ function.
|
||
|
||
!bc pycod
|
||
# Importing various packages
|
||
import numpy as np
|
||
n = 100
|
||
x = np.random.normal(size=n)
|
||
print(np.mean(x))
|
||
y = 4+3*x+np.random.normal(size=n)
|
||
print(np.mean(y))
|
||
W = np.vstack((x, y))
|
||
C = np.cov(W)
|
||
print(C)
|
||
!ec
|
||
|
||
!split
|
||
===== Correlation Matrix =====
|
||
|
||
The previous example can be converted into the correlation matrix by
|
||
simply scaling the matrix elements with the variances. We should also
|
||
subtract the mean values for each column. This leads to the following
|
||
code which sets up the correlations matrix for the previous example in
|
||
a more brute force way. Here we scale the mean values for each column of the design matrix, calculate the relevant mean values and variances and then finally set up the $2\times 2$ correlation matrix (since we have only two vectors).
|
||
|
||
!bc pycod
|
||
import numpy as np
|
||
n = 100
|
||
# define two vectors
|
||
x = np.random.random(size=n)
|
||
y = 4+3*x+np.random.normal(size=n)
|
||
#scaling the x and y vectors
|
||
x = x - np.mean(x)
|
||
y = y - np.mean(y)
|
||
variance_x = np.sum(x@x)/n
|
||
variance_y = np.sum(y@y)/n
|
||
print(variance_x)
|
||
print(variance_y)
|
||
cov_xy = np.sum(x@y)/n
|
||
cov_xx = np.sum(x@x)/n
|
||
cov_yy = np.sum(y@y)/n
|
||
C = np.zeros((2,2))
|
||
C[0,0]= cov_xx/variance_x
|
||
C[1,1]= cov_yy/variance_y
|
||
C[0,1]= cov_xy/np.sqrt(variance_y*variance_x)
|
||
C[1,0]= C[0,1]
|
||
print(C)
|
||
!ec
|
||
|
||
We see that the matrix elements along the diagonal are one as they
|
||
should be and that the matrix is symmetric. Furthermore, diagonalizing
|
||
this matrix we easily see that it is a positive definite matrix.
|
||
|
||
The above procedure with _numpy_ can be made more compact if we use _pandas_.
|
||
|
||
!split
|
||
===== Correlation Matrix with Pandas =====
|
||
|
||
We whow here how we can set up the correlation matrix using _pandas_, as done in this simple code
|
||
!bc pycod
|
||
import numpy as np
|
||
import pandas as pd
|
||
n = 10
|
||
x = np.random.normal(size=n)
|
||
x = x - np.mean(x)
|
||
y = 4+3*x+np.random.normal(size=n)
|
||
y = y - np.mean(y)
|
||
X = (np.vstack((x, y))).T
|
||
print(X)
|
||
Xpd = pd.DataFrame(X)
|
||
print(Xpd)
|
||
correlation_matrix = Xpd.corr()
|
||
print(correlation_matrix)
|
||
!ec
|
||
|
||
|
||
We expand this model to the Franke function discussed above.
|
||
|
||
!split
|
||
===== Correlation Matrix with Pandas and the Franke function =====
|
||
|
||
!bc pycod
|
||
# Common imports
|
||
import numpy as np
|
||
import pandas as pd
|
||
|
||
|
||
def FrankeFunction(x,y):
|
||
term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))
|
||
term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))
|
||
term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))
|
||
term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)
|
||
return term1 + term2 + term3 + term4
|
||
|
||
|
||
def create_X(x, y, n ):
|
||
if len(x.shape) > 1:
|
||
x = np.ravel(x)
|
||
y = np.ravel(y)
|
||
|
||
N = len(x)
|
||
l = int((n+1)*(n+2)/2) # Number of elements in beta
|
||
X = np.ones((N,l))
|
||
|
||
for i in range(1,n+1):
|
||
q = int((i)*(i+1)/2)
|
||
for k in range(i+1):
|
||
X[:,q+k] = (x**(i-k))*(y**k)
|
||
|
||
return X
|
||
|
||
|
||
# Making meshgrid of datapoints and compute Franke's function
|
||
n = 4
|
||
N = 100
|
||
x = np.sort(np.random.uniform(0, 1, N))
|
||
y = np.sort(np.random.uniform(0, 1, N))
|
||
z = FrankeFunction(x, y)
|
||
X = create_X(x, y, n=n)
|
||
|
||
Xpd = pd.DataFrame(X)
|
||
# subtract the mean values and set up the covariance matrix
|
||
Xpd = Xpd - Xpd.mean()
|
||
covariance_matrix = Xpd.cov()
|
||
print(covariance_matrix)
|
||
!ec
|
||
|
||
We note here that the covariance is zero for the first rows and
|
||
columns since all matrix elements in the design matrix were set to one
|
||
(we are fitting the function in terms of a polynomial of degree $n$).
|
||
|
||
This means that the variance for these elements will be zero and will
|
||
cause problems when we set up the correlation matrix. We can simply
|
||
drop these elements as follows and then construct the correlation
|
||
matrix.
|
||
|
||
|
||
|
||
!split
|
||
===== Classical PCA Theorem =====
|
||
|
||
|
||
!split
|
||
===== Prof of the PCA Theorem =====
|
||
|
||
|
||
|
||
|
||
|
||
|
||
!split
|
||
===== Getting started with PCA =====
|
||
|
||
|
||
!bc pycod
|
||
# Now add PCA
|
||
from sklearn.decomposition import PCA
|
||
pca = PCA(n_components = 2)
|
||
pca.fit(X_train_scaled)
|
||
|
||
X_pca = pca.transform(X_train_scaled)
|
||
!ec
|
||
|
||
|
||
|
||
!split
|
||
===== Principal Component Analysis =====
|
||
!bblock
|
||
Principal Component Analysis (PCA) is by far the most popular dimensionality reduction algorithm.
|
||
First it identifies the hyperplane that lies closest to the data, and then it projects the data onto it.
|
||
|
||
The following Python code uses NumPy’s _svd()_ function to obtain all the principal components of the
|
||
training set, then extracts the first two principal components
|
||
!bc pycod
|
||
X_centered = X - X.mean(axis=0)
|
||
U, s, V = np.linalg.svd(X_centered)
|
||
c1 = V.T[:, 0]
|
||
c2 = V.T[:, 1]
|
||
!ec
|
||
|
||
PCA assumes that the dataset is centered around the origin. Scikit-Learn’s PCA classes take care of centering
|
||
the data for you. However, if you implement PCA yourself (as in the preceding example), or if you use other libraries, don’t
|
||
forget to center the data first.
|
||
|
||
Once you have identified all the principal components, you can reduce the dimensionality of the dataset
|
||
down to $d$ dimensions by projecting it onto the hyperplane defined by the first $d$ principal components.
|
||
Selecting this hyperplane ensures that the projection will preserve as much variance as possible.
|
||
!bc pycod
|
||
W2 = V.T[:, :2]
|
||
X2D = X_centered.dot(W2)
|
||
!ec
|
||
|
||
!split
|
||
===== PCA and scikit-learn =====
|
||
|
||
Scikit-Learn’s PCA class implements PCA using SVD decomposition just like we did before. The
|
||
following code applies PCA to reduce the dimensionality of the dataset down to two dimensions (note
|
||
that it automatically takes care of centering the data):
|
||
!bc pycod
|
||
from sklearn.decomposition import PCA
|
||
pca = PCA(n_components = 2)
|
||
X2D = pca.fit_transform(X)
|
||
!ec
|
||
After fitting the PCA transformer to the dataset, you can access the principal components using the
|
||
components variable (note that it contains the PCs as horizontal vectors, so, for example, the first
|
||
principal component is equal to
|
||
!bc pycod
|
||
pca.components_.T[:, 0]).
|
||
!ec
|
||
Another very useful piece of information is the explained variance ratio of each principal component,
|
||
available via the $explained\_variance\_ratio$ variable. It indicates the proportion of the dataset’s
|
||
variance that lies along the axis of each principal component.
|
||
More material to come here.
|
||
|
||
!split
|
||
===== More on the PCA =====
|
||
|
||
Instead of arbitrarily choosing the number of dimensions to reduce down to, it is generally preferable to
|
||
choose the number of dimensions that add up to a sufficiently large portion of the variance (e.g., 95%).
|
||
Unless, of course, you are reducing dimensionality for data visualization — in that case you will
|
||
generally want to reduce the dimensionality down to 2 or 3.
|
||
The following code computes PCA without reducing dimensionality, then computes the minimum number
|
||
of dimensions required to preserve 95% of the training set’s variance:
|
||
!bc pycod
|
||
pca = PCA()
|
||
pca.fit(X)
|
||
cumsum = np.cumsum(pca.explained_variance_ratio_)
|
||
d = np.argmax(cumsum >= 0.95) + 1
|
||
!ec
|
||
You could then set $n\_components=d$ and run PCA again. However, there is a much better option: instead
|
||
of specifying the number of principal components you want to preserve, you can set $n\_components$ to be
|
||
a float between 0.0 and 1.0, indicating the ratio of variance you wish to preserve:
|
||
!bc pycod
|
||
pca = PCA(n_components=0.95)
|
||
X_reduced = pca.fit_transform(X)
|
||
!ec
|
||
|
||
!split
|
||
===== Incremental PCA =====
|
||
|
||
One problem with the preceding implementation of PCA is that it requires the whole training set to fit in
|
||
memory in order for the SVD algorithm to run. Fortunately, Incremental PCA (IPCA) algorithms have
|
||
been developed: you can split the training set into mini-batches and feed an IPCA algorithm one minibatch
|
||
at a time. This is useful for large training sets, and also to apply PCA online (i.e., on the fly, as new
|
||
instances arrive).
|
||
|
||
!split
|
||
===== Randomized PCA =====
|
||
|
||
Scikit-Learn offers yet another option to perform PCA, called Randomized PCA. This is a stochastic
|
||
algorithm that quickly finds an approximation of the first d principal components. Its computational
|
||
complexity is $O(m \times d^2)+O(d^3)$, instead of $O(m \times n^2) + O(n^3)$, so it is dramatically faster than the
|
||
previous algorithms when $d$ is much smaller than $n$.
|
||
|
||
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== Kernel PCA =====
|
||
!bblock
|
||
|
||
The kernel trick is a mathematical technique that implicitly maps instances into a
|
||
very high-dimensional space (called the feature space), enabling nonlinear classification and regression
|
||
with Support Vector Machines. Recall that a linear decision boundary in the high-dimensional feature
|
||
space corresponds to a complex nonlinear decision boundary in the original space.
|
||
It turns out that the same trick can be applied to PCA, making it possible to perform complex nonlinear
|
||
projections for dimensionality reduction. This is called Kernel PCA (kPCA). It is often good at
|
||
preserving clusters of instances after projection, or sometimes even unrolling datasets that lie close to a
|
||
twisted manifold.
|
||
For example, the following code uses Scikit-Learn’s KernelPCA class to perform kPCA with an
|
||
!bc pycod
|
||
from sklearn.decomposition import KernelPCA
|
||
rbf_pca = KernelPCA(n_components = 2, kernel="rbf", gamma=0.04)
|
||
X_reduced = rbf_pca.fit_transform(X)
|
||
!ec
|
||
|
||
!eblock
|
||
|
||
|
||
!split
|
||
===== LLE =====
|
||
|
||
Locally Linear Embedding (LLE) is another very powerful nonlinear dimensionality reduction
|
||
(NLDR) technique. It is a Manifold Learning technique that does not rely on projections like the previous
|
||
algorithms. In a nutshell, LLE works by first measuring how each training instance linearly relates to its
|
||
closest neighbors (c.n.), and then looking for a low-dimensional representation of the training set where
|
||
these local relationships are best preserved (more details shortly).
|
||
|
||
|
||
|
||
!split
|
||
===== Other techniques =====
|
||
|
||
|
||
There are many other dimensionality reduction techniques, several of which are available in Scikit-Learn.
|
||
|
||
Here are some of the most popular:
|
||
* _Multidimensional Scaling (MDS)_ reduces dimensionality while trying to preserve the distances between the instances.
|
||
* _Isomap_ creates a graph by connecting each instance to its nearest neighbors, then reduces dimensionality while trying to preserve the geodesic distances between the instances.
|
||
* _t-Distributed Stochastic Neighbor Embedding_ (t-SNE) reduces dimensionality while trying to keep similar instances close and dissimilar instances apart. It is mostly used for visualization, in particular to visualize clusters of instances in high-dimensional space (e.g., to visualize the MNIST images in 2D).
|
||
* Linear Discriminant Analysis (LDA) is actually a classification algorithm, but during training it learns the most discriminative axes between the classes, and these axes can then be used to define a hyperplane onto which to project the data. The benefit is that the projection will keep classes as far apart as possible, so LDA is a good technique to reduce dimensionality before running another classification algorithm such as a Support Vector Machine (SVM) classifier discussed in the SVM lectures.
|
||
|
||
|
||
|
||
|
||
|
||
|
||
|
||
Here are other examples where we use the _DataFrame_ functionality to handle arrays, now with more interesting features for us, namely numbers. We set up a matrix
|
||
of dimensionality $10\times 5$ and compute the mean value and standard deviation of each column. Similarly, we can perform mathematial operations like squaring the matrix elements and many other operations.
|
||
!bc pycod
|
||
import numpy as np
|
||
import pandas as pd
|
||
from IPython.display import display
|
||
np.random.seed(100)
|
||
# setting up a 10 x 5 matrix
|
||
rows = 10
|
||
cols = 5
|
||
a = np.random.randn(rows,cols)
|
||
df = pd.DataFrame(a)
|
||
display(df)
|
||
print(df.mean())
|
||
print(df.std())
|
||
display(df**2)
|
||
!ec
|
||
|
||
Thereafter we can select specific columns only and plot final results
|
||
!bc pycod
|
||
df.columns = ['First', 'Second', 'Third', 'Fourth', 'Fifth']
|
||
df.index = np.arange(10)
|
||
|
||
display(df)
|
||
print(df['Second'].mean() )
|
||
print(df.info())
|
||
print(df.describe())
|
||
|
||
from pylab import plt, mpl
|
||
plt.style.use('seaborn')
|
||
mpl.rcParams['font.family'] = 'serif'
|
||
|
||
df.cumsum().plot(lw=2.0, figsize=(10,6))
|
||
plt.show()
|
||
|
||
|
||
df.plot.bar(figsize=(10,6), rot=15)
|
||
plt.show()
|
||
!ec
|
||
We can produce a $4\times 4$ matrix
|
||
!bc pycod
|
||
b = np.arange(16).reshape((4,4))
|
||
print(b)
|
||
df1 = pd.DataFrame(b)
|
||
print(df1)
|
||
!ec
|
||
and many other operations.
|
||
|