PCA

pca.py

#!/usr/bin/env python3

import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

# 1. calculate mu
# X = np.array([
#   [-1, 1],
#   [ 0, 0],
#   [ 1, 1]
# ])
X = np.array([
    [1, 0],
    [2, 0],
    [3, 0],
    [5, 6],
    [6, 6],
    [7, 6]
])
mean = np.mean(X, axis=0)
mu = np.repeat(mean[np.newaxis, :], X.shape[0], 0)
print("1. Matrix mu is: \n{0}".format(mu))

# 2. calculate `Sigma_hat` i.e. `COV(X)`
X_wave = X - mu
Sigma_hat = (1 / X.shape[0]) * np.dot(np.transpose(X_wave), X_wave)
print("2. Matrix Sigma_hat is: \n{0}".format(Sigma_hat))

# 3. calculate eigen-pairs
eigen_value, eigen_vectors = LA.eig(Sigma_hat)
print("3. Eigen pairs are:     \
      \n  eigen value: {0},    \
      \n  eigen vectors:\n{1}"
      .format(eigen_value, eigen_vectors))

# 4. reconstruct original data matrix
U_T = np.array([eigen_vectors[0]])
U = np.transpose(U_T)
X_hat = mu + np.dot(X_wave, np.dot(U, U_T))
print("4. Result is: \n {0}".format(X_hat))

# 5. visualization
f, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True)
ax1.scatter(X[:, 0], X[:, 1])
ax1.set_title('Original Data')
ax1.set(adjustable='box-forced', aspect='equal')
ax2.scatter(X_hat[:, 0], X_hat[:, 1])
ax2.set_title('PCA Reconstruction')
ax2.set(adjustable='box-forced', aspect='equal')
plt.show()