PCA - “na piechotę”
Jednym z zadań eksploracji danych jest redukcja wymiaru, czyli określenie, które ze składowych wektora obserwacji są nieistotne lub też jakie inne kombinacje składowych mogą się okazać przydatne do dalszej analizy. Standardową metodą redukcji wymiaru jest analiza składowych głównych (PCA - principal component analysis). Polega ona na znalezieniu nowego kierunku, który maksymalizuje wariancję zrzutowanych na niego obserwacji. Następnie szukamy kolejnego kierunku, również o jak największej wariancji, tyle, że ortogonalnego do poprzedniego etc. Okazuje się, że takie cechy odpowiadają wektorom własnym związanym z kolejnymi wartościami własnymi (począwszy od największej) macierzy kowarinacji.
Zaczynamy od przykładu ilustrującego podstawowe cechy PCA dla sztucznie wykreowanych danych postaci \(y + \eta = x + \xi\), gdzie \(\eta\) i \(\xi\) są szumem z rozkładu jednorodnego \(U[-1/2, 1/2]\).
import numpy as np
import pandas as pd
x = np.arange(-5, 5, 0.1)
y = x + np.random.uniform(-0.5, 0.5, len(x))
x += np.random.uniform(-0.5, 0.5, len(x))
df = pd.DataFrame({'x':x, 'y':y})import matplotlib.pyplot as plt
plt.axes().set_aspect('equal')
plt.scatter(df.x, df.y, s = 10)
plt.axhline(0, color = "gray")
plt.axvline(0, color = "gray")
plt.axline((0,0), slope = 1, color = "green")
plt.axline((0,0), slope = -1, color = "red")
plt.text(4, -0.5, "$x$", fontsize = 14)
plt.text(-1, 5, "$y$", fontsize = 14)
plt.text(4, 3, "$x_1$", fontsize = 14, color = "green")
plt.text(-4, 3, "$y_1$", fontsize = 14, color = "red")
plt.show()
Otrzymalismy spodziewany wykres, gdzie zależność \(y=x\) jest lekko zaburzona ale dalej bardzo wyraźna. Nie uciekając się do generycznej funkji realizującej PCA, sami możemy policzyć skłądowe główne poprzez wyznaczenie wektorów i wartości własnych macierzy kowariancji \(\mathbf{S}\):
S = np.cov(df, rowvar=False)
print(S)## [[8.77998917 8.56780868]
## [8.56780868 8.52687869]]lam, u = np.linalg.eig(S.T)
idx = lam.argsort()[::-1]
lam = lam[idx]
u = u[:,idx]
print(lam)## [17.22217723 0.08469063]print(u)## [[ 0.71230941 -0.70186558]
## [ 0.70186558 0.71230941]]lub macierzy korelacji \(\mathbf{S}_c\):
S_c = df.corr()
print(S_c)## x y
## x 1.000000 0.990211
## y 0.990211 1.000000lam_c, u_c = np.linalg.eig(S_c)
idx = lam_c.argsort()[::-1]
lam_c = lam_c[idx]
u_c = u_c[:,idx]
print(lam_c)## [1.99021096 0.00978904]print(u_c)## [[ 0.70710678 -0.70710678]
## [ 0.70710678 0.70710678]]Jak widać, wektory własne są bardzo podobne, jeśli chodzi wartości własne, to są one po prostu przeskalowane (tzn. \(\lambda_2 / lambda_1\) w obu przypadkach powinno być z grubsza takie samo). Rzecz jasna, jeśli znacząco zmienmi kowariancję (np. \(s_{xx}\) będzie znacząco różne od \(s_{yy}\)), to wektory własne powinny być bardziej niepodobne do siebie.
Mozemy teraz nanieść otrzymane dwa wektory na nasz rysunek, przy czym żeby wysunek był OK, musimy zadbać, aby skale były takie same.
import matplotlib.pyplot as plt
plt.axes().set_aspect('equal')
plt.scatter(df.x, df.y, s = 10)
plt.axhline(0, color = "gray")
plt.axvline(0, color = "gray")
plt.text(4, -0.5, "$x$", fontsize = 14)
plt.text(-1, 5, "$y$", fontsize = 14)
plt.arrow(0, 0, 2*u[0,0], 2*u[1,0], color = "orange", width = 0.1)
plt.arrow(0, 0, 2*u[0,1], 2*u[1,1], color = "orange", width = 0.1)
plt.show()
Zgodnie z naszym wykładem, nowe współrzędne powstają poprzez wykonanie mnożenia macierzowego \(\mathbf{\Gamma}^{T}\left(\mathbf{X} - \mathbf{m} \right)\), gdzie \(\mathbf{\Gamma}\) to macierz wektorów własnych, a \(\mathbf{X}\) to macierz naszych danych.
df1 = (df-df.mean()) @ u
print(df1)## 0 1
## 0 -6.661247 0.288764
## 1 -7.024944 0.078740
## 2 -6.466700 -0.015800
## 3 -6.153106 -0.142002
## 4 -6.836449 0.192935
## .. ... ...
## 95 6.627487 -0.478734
## 96 6.494207 -0.372547
## 97 7.235440 -0.027233
## 98 6.786699 -0.020024
## 99 6.959327 0.048557
##
## [100 rows x 2 columns]df1.columns = ['x1', 'y1']
plt.axes().set_aspect('equal')
plt.scatter(df1.x1, df1.y1, s = 10)
plt.axhline(0, color = "gray")
plt.axvline(0, color = "gray")
plt.text(7, -1, "$x_1$", fontsize = 14)
plt.text(-1, 5, "$y_1$", fontsize = 14)
plt.ylim((-6,6))## (-6.0, 6.0)plt.show()
Warto przy tym porównać otrzymane wcześniej wartości własne dla macierzy kowariancji oraz kowariancje nowych zmiennych:
print(lam)## [17.22217723 0.08469063]print(df1.var())## x1 17.222177
## y1 0.084691
## dtype: float64PCA - sklearn
Skorzystajmy teraz już z pakietu sklearn, a dokładniej z
modułu decomposition i funkcji PCA().
from sklearn.decomposition import PCA
pca = PCA(n_components=None)
pca.fit(df)PCA()In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
Parameters
Zacznimy od wypisania tego, co już zrobiliśmy, czyli wariancji
pca.explained_variance_## array([17.22217723, 0.08469063])Bardzo często będziemy potrzebować nie tyle samej wariancji, ale
informacji o tym, jak dużą część całkowitej wariancji (nowych zmiennych)
wyjaśnia dany kierunek. Do tego możemy wykorzystać bezposrednio
explained_variance_ratio_
pca.explained_variance_ratio_## array([0.99510653, 0.00489347])lub jeszcze bardziej wartość skumulowaną:
plt.bar(["PC " + str(x) for x in range(1, 3)], np.cumsum(pca.explained_variance_ratio_));
plt.show()
Powyższy przypadek nie jest może specjalnie interesujący, ale posłużmy się zbiorem irysów:
from sklearn import datasets
iris = datasets.load_iris()
iris_pca = pca.fit(iris.data)
plt.bar(["PC " + str(x) for x in range(1, 5)], iris_pca.explained_variance_ratio_);
plt.show()
Wyraźnie widać, że pierwsza składowa główna praktycznie wystarcza do
opisu danych. Możemy też wykorzystać pole components_ aby
przekonać się jaki wkład mają oryginalne składowe do skłądowych
głównych:
iris_pca.components_## array([[ 0.36138659, -0.08452251, 0.85667061, 0.3582892 ],
## [ 0.65658877, 0.73016143, -0.17337266, -0.07548102],
## [-0.58202985, 0.59791083, 0.07623608, 0.54583143],
## [ 0.31548719, -0.3197231 , -0.47983899, 0.75365743]])
import seaborn as sns
sns.heatmap(iris_pca.components_, cmap = "bwr")
kPCA - nieliniowa PCA
W wielu przypadkach zwykła analiza składowych głównych nie daje pożądanego rezultatu. Wylosujemy dość perfidny zbiór danych, który zdecydowanie nie jest liniowy:
from sklearn.datasets import make_circles
data, target = make_circles(n_samples=1_000, factor = 0.3, noise = 0.05)
plt.scatter(data[:,0], data[:,1], c = target, cmap=plt.cm.coolwarm)
Jeśli wykonamy PCA, to efekt nie bedzie specjalnie różny od oryginału:
data_pca = pca.fit_transform(data)
plt.scatter(data_pca[:,0], data_pca[:,1], c = target, cmap=plt.cm.coolwarm)
W takiej sytuacji warto skorzystać z nieliniowej metody PCA - w tym przypadku z jądrem radialnym sterowanym parametrem \(\gamma\)
from sklearn.decomposition import KernelPCA
kernel_pca = KernelPCA(
n_components=None, kernel="rbf", gamma=2)
data_kpca = kernel_pca.fit_transform(data)
plt.scatter(data_kpca[:, 0], data_kpca[:, 1], c = target, cmap=plt.cm.coolwarm)
Skalowanie wielowymiarowe (MDS)
Bardzo ciekawą i prostą w wizualizacji metodą jest skalowanie wielowymiarowe (ang. multidimensional scaling - MDS). W ogólnym ujęciu, za pomocą tego podejścia możliwe jest przedstawienie odległości pomiędzy poszczególnymi punktami (obserwacjami) w niskowymiarowej (np 2D) przestrzeni. Jest to szczególnie proste w przypadku danych geograficznych, dla których posiadamy macierz wzajemnych odległości:
euro = pd.read_csv("https://jsienkiewicz.pl/pliki/SDE/eurodist.csv", index_col = 0)
euro## Athens Barcelona Brussels ... Rome Stockholm Vienna
## Athens 0 3313 2963 ... 817 3927 1991
## Barcelona 3313 0 1318 ... 1460 2868 1802
## Brussels 2963 1318 0 ... 1511 1616 1175
## Calais 3175 1326 204 ... 1662 1786 1381
## Cherbourg 3339 1294 583 ... 1794 2196 1588
## Cologne 2762 1498 206 ... 1497 1403 937
## Copenhagen 3276 2218 966 ... 2050 650 1455
## Geneva 2610 803 677 ... 995 2068 1019
## Gibraltar 4485 1172 2256 ... 2631 3886 2974
## Hamburg 2977 2018 597 ... 1751 949 1155
## Hook of Holland 3030 1490 172 ... 1683 1500 1205
## Lisbon 4532 1305 2084 ... 2700 3231 2937
## Lyons 2753 645 690 ... 1048 2108 1157
## Madrid 3949 636 1558 ... 2097 3188 2409
## Marseilles 2865 521 1011 ... 1011 2428 1363
## Milan 2282 1014 925 ... 586 2187 898
## Munich 2179 1365 747 ... 946 1754 428
## Paris 3000 1033 285 ... 1476 1827 1249
## Rome 817 1460 1511 ... 0 2707 1209
## Stockholm 3927 2868 1616 ... 2707 0 2105
## Vienna 1991 1802 1175 ... 1209 2105 0
##
## [21 rows x 21 columns]W tym momencie możemy już wykonać skalowanie:
from sklearn.manifold import MDS
embedding = MDS(n_components=2, dissimilarity ="euclidean")
euro_transformed = embedding.fit_transform(euro)## C:\Users\julas\MINICO~1\envs\sde\Lib\site-packages\sklearn\manifold\_mds.py:744: FutureWarning: The default value of `n_init` will change from 4 to 1 in 1.9. To suppress this warning, provide some value of `n_init`.
## warnings.warn(
## C:\Users\julas\MINICO~1\envs\sde\Lib\site-packages\sklearn\manifold\_mds.py:754: FutureWarning: The default value of `init` will change from 'random' to 'classical_mds' in 1.10. To suppress this warning, provide some value of `init`.
## warnings.warn(
## C:\Users\julas\MINICO~1\envs\sde\Lib\site-packages\sklearn\manifold\_mds.py:771: FutureWarning: The `dissimilarity` parameter is deprecated and will be removed in 1.10. Use `metric` instead.
## warnings.warn(
## C:\Users\julas\MINICO~1\envs\sde\Lib\site-packages\sklearn\manifold\_mds.py:794: UserWarning: The provided input is a square matrix. Note that ``fit`` constructs a dissimilarity matrix from data and will treat rows as samples and columns as features. To use a pre-computed dissimilarity matrix, set ``metric='precomputed'``.
## warnings.warn(plt.scatter(euro_transformed[:,0], euro_transformed[:,1])
for i, s in enumerate(euro.index):
plt.text(euro_transformed[i,0], euro_transformed[i,1], s = s)
plt.show()
W podobny sposób da się opracować dane, które nie mają żadnych cech faktycznych odległości - możemy się tu posłużyć zbiorem dotyczącym głosowań z poprzednich zajęć:
df = pd.read_csv("https://jsienkiewicz.pl/pliki/SDE/votes.csv", index_col=0)
votes = df.iloc[:,11:]
votes = votes.drop(votes[votes.apply(lambda x: x.isnull().values.any(), axis = 1)].index)
votes_transformed = embedding.fit_transform(votes)## C:\Users\julas\MINICO~1\envs\sde\Lib\site-packages\sklearn\manifold\_mds.py:744: FutureWarning: The default value of `n_init` will change from 4 to 1 in 1.9. To suppress this warning, provide some value of `n_init`.
## warnings.warn(
## C:\Users\julas\MINICO~1\envs\sde\Lib\site-packages\sklearn\manifold\_mds.py:754: FutureWarning: The default value of `init` will change from 'random' to 'classical_mds' in 1.10. To suppress this warning, provide some value of `init`.
## warnings.warn(
## C:\Users\julas\MINICO~1\envs\sde\Lib\site-packages\sklearn\manifold\_mds.py:771: FutureWarning: The `dissimilarity` parameter is deprecated and will be removed in 1.10. Use `metric` instead.
## warnings.warn(
plt.scatter(votes_transformed[:,0], votes_transformed[:,1])
for i, s in enumerate(votes.index):
plt.text(votes_transformed[i,0], votes_transformed[i,1], s = s)