Laten we in een concreet voorbeeld de principale componenten analyse in Python uitvoeren.
We beginnen met een dataset te creëren voor twee grootheden en de gegevens te visualiseren.

>>> import numpy as np
>>> import numpy.linalg as la
>>> import matplotlib.pyplot as plt
>>> X1 = np.random.normal(0, 1, 25)
>>> X1 = np.sort(X1)
>>> X2 = 2 + 1.5*X1 + 0.25*np.random.normal(0, 1, 25)
>>> plt.plot(X1, X2, 'bo', mfc='none')
>>> plt.gca().set_aspect('equal')
>>> plt.show()

Om de analyse goed te kunnen uitvoeren moeten we de gegevens eerst bewerken zodanig dat de gemiddelde waarden van de twee grootheden gelijk aan nul zijn.

>>> X1 = X1 - np.mean(X1)
>>> X2 = X2 - np.mean(X2)
>>> X = np.vstack((X1,X2)))

We kunnen voor de datamatrix \(X\) met 2 grootheden (2 rijen) en \(N=25\) gegevens (25 kolommen) dan de covariantiematrix \(C\) berekenen als \[C = \frac{1}{N-1}X\,X^{\top}\] Het numpy pakket heeft hiervoor de functie cov aan boord.

>>> C = np.cov(X); print(C)
[[1.48630631 2.31901205]
[2.31901205 3.66047268]]

We willen nu nieuwe grootheden definiëren wier waarden ongecorreleerd zijn. Hiervoor diagonaliseren we de covariantiematrix door het eigenwaardenprobleem voor de matrix \(C\) numeriek op te lossen.

>>> [vals, vecs] = la.eig(C)
>>> print(vals)
[0.01222438 5.13455461]
>>> print(vecs)
[[-0.84393385 -0.53644725]
[ 0.53644725 -0.84393385]]

We sorteren nu de eigenwaarde in dalende volgorde en sorteren de bijpassende eigenvectoren mee.

>>> sorted_index = np.argsort(vals)[::-1]
>>> sorted_vals = vals[sorted_index]
>>> sorted_vecs = vecs[:, sorted_index]

De eigenvectoren definiëren nieuwe grootheden (als lineaire combinaties van de oude grootheden) die ongecorreleerd zijn omdat correlatie niets ander is dan inproduct gelijk aan nul voor de kolomvectoren. In onderstaande figuur staan de twee vectoren afgebeeld als stukje van een nieuw coordinatensysteem.

>>> plt.plot(X1, X2, 'bo', mfc='none')
>>> plt.arrow(0,0,-sorted_vecs[0,0],-sorted_vecs[1,0])
>>> plt.arrow(0,0,sorted_vecs[0,1],sorted_vecs[1,1])
>>> plt.gca().set_aspect('equal')
>>> plt.show()

We kunnen de gegevens transformeren voor het nieuwe coördinatensysteem en visualiseren. Omdat de eigenwaarden in opklimmende volgorde staan moeten we tussendoor kolommen verwisselen.

>>> Y = np.transpose(sorted_vecs) @ X
>>> Y1 = -Y[0,:]  # minteken vanwege tegengestelde eerste gesorteerde eigenvector 
>>> Y2 = Y[1,:]
>>> plt.plot(Y1, Y2, 'bo', mfc='none')
[<matplotlib.lines.Line2D object at 0x000001D11480EA90>]
>>> plt.gca().set_aspect('equal')
>>> plt.show()

Met de functie var kunnen we de variatie in de grootheden bepalen. De ene grootheid betreft dan \(99.8\%\) van de variatie van de data, en de andere grootheid slechts voor \(0.2\%\).

>>> varY1 = np.var(Y1); round(varY1, 5)
4.92917
>>> varY2 = np.var(Y2); round(varY2, 5)
0.01174
>>> sumvar = varY1 + varY2
>>> varsY = [varY1/sumvar, varY2/sumvar]; 
>>> print(np.round(varsY, 3))
[0.998 0.002]

\(\phantom{x}\)

We kunnen natuurlijk ook ingebouwde functies voor PCA gebruiken; we passen hier pca uit het sklearn pakket toe en laten zien dat we dan hetzelfde resultaat krijgen.

>>> from sklearn.decomposition import PCA
>>> pca = PCA(2)
>>> pca.fit(X)
PCA(n_components=2)
>>> pca.fit(np.transpose(X))
PCA(n_components=2)
>>> print(pca.components_) # coëfficiënten van de principale componenten (factorladingen)
[[ 0.53644725 0.84393385]
[-0.84393385 0.53644725]]
>>> print(pca.explained_variance_)
[5.13455461 0.01222438]
>>> print(np.round(pca.explained_variance_ /sum(pca.explained_variance_), 3))
[0.998 0.002]
>>> YY = pca.transform(np.transpose(X))  # getransformeerde data
>>> round(sum(abs(Y1 - YY[:,0]))) # Y1 = eerste kolom van YY
0
>>> round(sum(abs(Y2 - YY[:,1]))) # Y2 = tweede kolom van YY
0