Principal Component Analysis¶
Intro¶
Is a procedure that uses orthogonal transformation to convert correlated variables into linearly uncorrelated variables, which it calls the principal components. It does this by repeatedly determining which variables account/explain the greatest amount of variability in the data and is orthogonal to the preceding components. (This constraint does not apply to the first component). The result is a set of vectors that forms an uncorrelated orthogonal basis set for the data
For those needing a refresher (like i did)
Correlation := a standard measure of the interdependence between 2 variables
Orthogonal := statistical context implies independence which is the same as two vectors being at 90-degrees to each other. The orthogonal basis set is a potential axis set, aka the principal axes. I like to think of these as the x,y,z axes in a 3-dimensional diagram
In Machine learning application it is used to reduce the number of dimensions. Higher dimensions can be easier to work with mathematically but considerably more difficult to apply to the real world.
Basically it's a great way to take a lot of data and reduce it to the important stuff.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()
rng = np.random.RandomState(1)
X = np.dot(rng.rand(2, 3), rng.randn(3, 250)).T
plt.scatter(X[:, 0], X[:, 1])
plt.axis('equal');
# Obviously this dataset is highly correlated.
# Instead of trying to predict y from x we shall try to learn about their relationship.
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
pca.fit(X)
print("Principal Components: ")
print(pca.components_)
print("Explained Variance: ")
print(pca.explained_variance_)
# The above appears rather cryptic so let's try to visualize it
def draw_vector(v0, v1, ax=None):
ax = ax or plt.gca()
arrowprops=dict(arrowstyle='->',
linewidth=2,
shrinkA=0, shrinkB=0)
ax.annotate('', v1, v0, arrowprops=arrowprops)
# plot data
plt.scatter(X[:, 0], X[:, 1], alpha=0.2)
for length, vector in zip(pca.explained_variance_, pca.components_):
v = vector * 3 * np.sqrt(length)
draw_vector(pca.mean_, pca.mean_ + v)
plt.axis('equal');
# The axes above are the principal axes for the dataset
# the length of the axis/vector is a reflection of it's importance in explaining the variance
# and thus describing the distribution
#Now let's see what happens when we reduce the dimensionality
# recreate the original diagram and draw the vectors
pca = PCA(n_components=2)
pca.fit(X)
pca2 = PCA(n_components=1)
pca2.fit(X)
X_pca = pca2.transform(X)
#print("original shape: ", X.shape)
#print("transformed shape:", X_pca.shape)
X_new = pca2.inverse_transform(X_pca)
plt.scatter(X[:, 0], X[:, 1], alpha=0.2)
plt.scatter(X_new[:, 0], X_new[:, 1], c='red', alpha=0.8)
for length, vector in zip(pca.explained_variance_, pca.components_):
v = vector * 3 * np.sqrt(length)
draw_vector(pca.mean_, pca.mean_ + v)
plt.axis('equal');
What we've done is create 2 models.
Model 1: Retains both dimensions and computes the principal component vectors
Model 2:
Sets the numbers of components=1 thereby reducing the number of dimensions(2).
Transforms the data (zero's out the less valuable component) and stores in X_pca
Apply inverse transform (which will extrapolate the Ys from X_pca)
The end result here is that the final extrapolated data (the red points) contains the majority of the information in the original but came at a much lower price point.
RW Example 1¶
It can be difficult to understand the application of PCA from the concocted data above so let's take a look at something in higher dimensions.
The following example can be found at:
https://www.kaggle.com/nirajvermafcb/principal-component-analysis-explained/notebook
import pandas as pd
raw = pd.read_csv('./data/hr_comma_sep.csv')
print('Shape: ', raw.shape)
columns_names = raw.columns.tolist()
print("Columns names:")
print(columns_names)
raw.head()
correlation = raw.corr()
plt.figure(figsize=(10,10))
sns.heatmap(correlation, vmax=1, square=True,annot=True,cmap='cubehelix')
raw2 = raw.drop(labels=['sales','salary'],axis=1) # drop non numeric data
# Move target column to the front
cols = raw2.columns.tolist()
cols.insert(0, cols.pop(cols.index('left')))
raw2 = raw2.reindex(columns= cols)
#raw2.columns.tolist()
#Seperate features and labels
X = raw2.iloc[:,1:8].values #Cols 1 thru 8 inclusive
y = raw2.iloc[:,0].values #Col 0 = left
X
from sklearn.preprocessing import StandardScaler
#Standardize features by removing the mean and scaling to unit variance
X_std = StandardScaler().fit_transform(X)
#Get mean vector
mean_vec = np.mean(X_std, axis=0)
#Get feature var-covar matrix
cov_mat = (X_std - mean_vec).T.dot((X_std - mean_vec)) / (X_std.shape[0]-1)
print('Covariance matrix \n%s' %cov_mat)
plt.figure(figsize=(8,8))
sns.heatmap(cov_mat, vmax=1, square=True,annot=True,cmap='cubehelix')
plt.title('Correlation between features')
# Eigenvector/Value decomposition
eig_vals, eig_vecs = np.linalg.eig(cov_mat)
print('Eigenvectors \n%s' %eig_vecs)
print('\nEigenvalues \n%s' %eig_vals)
# Selecting Principal Components
# 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.sort(key=lambda x: x[0], reverse=True)
# Visually confirm that the list is correctly sorted by decreasing eigenvalues
#print('Eigenvalues in descending order:')
#for i in eig_pairs:
# print(i[0])
tot = sum(eig_vals)
var_exp = [(i / tot)*100 for i in sorted(eig_vals, reverse=True)]
with plt.style.context('dark_background'):
plt.figure(figsize=(6, 4))
plt.bar(range(7), var_exp, alpha=0.5, align='center',
label='individual explained variance')
plt.ylabel('Explained variance ratio')
plt.xlabel('Principal components')
plt.legend(loc='best')
plt.tight_layout()
matrix_w = np.hstack((eig_pairs[0][1].reshape(7,1),
eig_pairs[1][1].reshape(7,1)
))
print('Matrix W:\n', matrix_w)
Y = X_std.dot(matrix_w)
Y
from sklearn.decomposition import PCA
pca = PCA().fit(X_std)
plt.plot(np.cumsum(pca.explained_variance_ratio_))
plt.xlim(0,7,1)
plt.xlabel('Number of components')
plt.ylabel('Cumulative explained variance')
#above plot shows almost 90% variance by the first 6 components. Therfore we can drop 7th component.
from sklearn.decomposition import PCA
sklearn_pca = PCA(n_components=6)
Y_sklearn = sklearn_pca.fit_transform(X_std)
print(Y_sklearn.shape)
print(Y_sklearn)