В данной заметке подйет речь о преобразования над самими признаками, т.е. над матрицей X. Значения каждого из прецендента нам не важно. Постараемся из X выжать как можно больше данных.
Это предварительная версия! Любые замечания приветсвуются.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
Метод Главных Компонент
Исходя из предыдущих заметок матрица X признаков не должна содержать избыточные признаки (столбцы матрицы X). Избыточные признаки это те, что могут быть восстановлены по некой формуле из других. В частности, если некий столбец может быть восстановлен как линейная комбинация других столбцов. В последнем случае столбцы будут линейно зависимоми и матрицы будет вырождена, даже в псевдообратном случае. Поэтому, необходимо такие столбцы выявить и устранить.
Синтезированные примеры
Формирование данных на приямой, которая под углом.
Для начала рассмотрим синтезированные данные. Причем сначала вообще идеальные "одномерные", т.е. по сути одномерные, то лежащие в двумерной плоскости.
In [2]:
# Сгенерируем данные
x = np.random.uniform(-10, 10, 30) # Раномерная выборка на отрезке [-10, 10] в количестве 30 элементов.
x.shape
Out[2]:
In [3]:
# Две строчки -- к каждому числу добавлена вторая координата (равная 0).
X = np.row_stack( (x, np.zeros_like(x) ) )
X # Выводим матрицу.
Out[3]:
In [4]:
plt.plot(X[0,:], X[1,:], ".") # Нарисуем их
# Они находятся на прямой.
Out[4]:
In [5]:
# Теперь подкрутим данные. Разместим их под уголом.
theta = np.radians( 30 ) # Градусы в радианы.
c, s = np.cos(theta), np.sin(theta) # Поэлементвное присвоение,
# Фактически эквивалентно:
# c = np.cos(theta)
# s = np.sin(theta)
In [6]:
R = np.array([[c,-s], [s, c]]) # Формируем матрицу поворота,
T = np.array([[2.5], [-1.3]]) # а также вектор смещения.
In [7]:
XX = R.dot(X) + T # Применяем к нашим данным вращение и трансляцию.
XXt = XX.transpose()
XXt
Out[7]:
In [8]:
plt.plot(XXt[:,0], XXt[:,1], "*", 0, 0, "r.")
# Данные теперь под углом. Ордината (красная точка) находится сбоку.
Out[8]:
Метод главных компонент (PCA), данные на прямой.
Как восстановить их? Как устранить факт того, что они под уголом? Для этого есть метод главных компонент.
Он позволяет снизить размерность данных, т.е. найти подпространство, которое содержит данные. Самое пространство и преобразование линейны.
In [9]:
# Загружаем метод главных компонент (PCA -- Principal component analysis):
from sklearn.decomposition import PCA
In [10]:
# PCA это конструктор объектов. Для выполнения вычислений нужно объект создать.
# Самое главное задать требуемое размер пространства.
pca = PCA( n_components = 1 ) # Хотим преобзовать данные к прямой, т.е. сделать их одномерными (как были исходные).
In [11]:
pca.fit( XXt ) # Запускаем метод pca на данных XXt, т.е. ищим подплоскость (смещение и поворот) содержащую данные.
# Выводятся используемые значения параметров метода.
Out[11]:
In [12]:
XXtd = pca.transform( XXt ) # Применяем найденное (сохраненное в объекте pca) преобразование к данным.
XXtd
Out[12]:
In [13]:
pca.explained_variance_ratio_ # Хорошо-ли данные коррелируют с вычисленным преобзованием?
# Фактически говорит о том насколько хорошо удалось снизить размерность,
# много ли данных оказалось искаженными (спроецированными с сильным смещением)
# 1 -- 100% по данных по первой оси (она же единственная).
Out[13]:
In [14]:
# Теперь нужно сравнить восстановленные данные с исходными. Для этого их нужно привести к нормальному виду:
from sklearn import preprocessing # preprocessing -- первичные преобразования данных.
# Будем применять preprocessing.scale -- делает выборочное среднее равно нулю, а дисперсию равной 1.
# Иногда, для корретной работы примера, требуется:
XXtd = -XXtd # см ниже.
XXtds = preprocessing.scale( XXtd[:, 0] ) # Нормируем. Используем функию scale из модуля sklearn.preprocessing.
XXtds
Out[14]:
In [15]:
xs = preprocessing.scale( x ) # x -- это исходные даные. Их тоже нормируем.
xs
Out[15]:
In [16]:
# Вычислим расстояние между векторами, т.е. отнормированными данными.
np.linalg.norm((xs - XXtds)) # Вызываем функцию из модуля, без from...
# Длина фактически равна нулю. Значит данные были восстановлены полностью.
# Если норма не равна нулю, то требуется массив отразить. см выше: XXtd = -XXtd. Нужно убрать комментарий.
Out[16]:
In [17]:
# Что будет если на две оси спроецировать?
pca = PCA( n_components = 2 ) # Настраиваем объект на двумерную плоскость.
XXtd = pca.fit_transform( XXt ) # Настроим параметры модели и сразу преобразуем данные к двум компонентам.
XXtd # Видим, что даные состоят из двумерных точек. Вторая компонента фактически равна нулю (машинное эпсилон).
Out[17]:
In [18]:
# Построим преобращованные (с целью снижения размерность) точки.
plt.plot(XXtd[:,0], XXtd[:,1], "*" ) # Но график будет обманчевым. Ось по y отмасштабирована (разные масштабы)!
Out[18]:
In [19]:
plt.plot(XXtd[:,0], XXtd[:,1], "*" )
plt.axis('equal') # Деламем оси одинакового масштаба.
# Правильнее так. Теперь видно, что данные расположены вдоль оси x.
Out[19]:
In [20]:
pca.explained_variance_ratio_ # Теперь тоже проверяем насколько хорошо данные коррелируют с осями.
# Видно, что первая ось все содержит. Для второй оси коеффициент равен 0.
Out[20]:
Видно, что увеличиение количества осей мы не теряем информацию о том, что данные на самом деле одномерные.
Метод главных компонент (PCA), данные на прямой плюс шум
In [21]:
# Добавим шума (к y, раньше он был раен 0)
X[1,:] = np.random.randn( (30) )/4 # вторая компонента содержит номальный шум, раньше там был 0.
XX = R.dot(X) + T # Как и раньше применяем вращение и трансляцию.
XXt = XX.transpose() # Для удобства. Например, при вызове plot, fit_transform.
In [22]:
plt.plot(XXt[:,0], XXt[:,1], "*", 0, 0, "r.")
Out[22]:
In [23]:
pca = PCA( n_components = 1 ) # Настраиваем объект на одну компоненту.
XXtd = pca.fit_transform( XXt ) # Ищим параметры преобразования и преобразуем.
In [24]:
pca.explained_variance_ratio_ # Теперь важность первой оси немного упала. Но совсем немного!
Out[24]:
In [25]:
# Теперь спроецируем на две оси.
pca = PCA( n_components = 2 )
XXtd = pca.fit_transform( XXt )
XXtd # Видно, что вторая компонента фактически равна шуму.
Out[25]:
In [26]:
plt.plot(XXtd[:,0], XXtd[:,1], "*" ) # Но график как и раньше будет обманный. Ось по y отмасштабирована!
Out[26]:
In [27]:
plt.plot(XXtd[:,0], XXtd[:,1], "*" ) # Правильнее так. Теперь видно, что данные расположены вдоль оси x.
plt.axis('equal') # Деламем оси одинакового масштаба. Вдоль y расположен шум.
Out[27]:
In [28]:
pca.explained_variance_ratio_ # Теперь тоже проверяем насколько хорошо данные коррелируют с осями.
# Видно, что первая ось содержит почти все. Для второй оси коеффициент почти равен 0, но не 0.
Out[28]:
In [29]:
XXtd[:,1] # Это шум.
Out[29]:
In [30]:
import scipy.stats as models
models.norm.fit( XXtd[:,1] )[1], models.norm.fit( X[1,:] )[1]
# Почти как исходное 1/4 = 0.25. У исходного и преобразованного оценка почти совпадает.
Out[30]:
Аналогичный пример в 3d
данные находятся в плоскости.
In [31]:
from sklearn.decomposition import PCA
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
In [32]:
e = np.exp(1)
np.random.seed(4)
def pdf(x):
return 0.5 * (stats.norm(scale=0.25 / e).pdf(x)
+ stats.norm(scale=4 / e).pdf(x))
y = np.random.normal(scale=0.5, size=(30000))
x = np.random.normal(scale=0.5, size=(30000))
z = np.random.normal(scale=0.1, size=len(x))
density = pdf(x) * pdf(y)
pdf_z = pdf(5 * z)
density *= pdf_z
a = x + y
b = 2 * y
c = a - b + z
norm = np.sqrt(a.var() + b.var())
a /= norm
b /= norm
In [33]:
density.shape
Out[33]:
In [34]:
def plot_figs(fig_num, elev, azim):
fig = plt.figure(fig_num, figsize=(8, 6))
plt.clf()
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=elev, azim=azim)
ax.scatter(a[::10], b[::10], c[::10], c=density[::10], marker='+', alpha=.4)
Y = np.c_[a, b, c]
# Using SciPy's SVD, this would be:
# _, pca_score, V = scipy.linalg.svd(Y, full_matrices=False)
pca = PCA(n_components=3)
pca.fit(Y)
pca_score = pca.explained_variance_ratio_
V = pca.components_
print( V.shape )
x_pca_axis, y_pca_axis, z_pca_axis = V.T * pca_score / pca_score.min()
x_pca_axis, y_pca_axis, z_pca_axis = 3 * V.T
x_pca_plane = np.r_[x_pca_axis[:2], - x_pca_axis[1::-1]]
y_pca_plane = np.r_[y_pca_axis[:2], - y_pca_axis[1::-1]]
z_pca_plane = np.r_[z_pca_axis[:2], - z_pca_axis[1::-1]]
x_pca_plane.shape = (2, 2)
y_pca_plane.shape = (2, 2)
z_pca_plane.shape = (2, 2)
print(x_pca_plane)
ax.plot_surface(x_pca_plane, y_pca_plane, z_pca_plane)
ax.w_xaxis.set_ticklabels([])
ax.w_yaxis.set_ticklabels([])
ax.w_zaxis.set_ticklabels([])
In [35]:
elev = -40
azim = -80
plot_figs(1, elev, azim)
In [36]:
elev = 30 # Так можно сменить точки из котрой смотрим на 3д график.
azim = 20
plot_figs(2, elev, azim)
Другие примеры с PCA
In [37]:
X = np.array([[-1, -1], [-2, -1], [-3, -2], [1, 1], [2, 1], [3, 2]])
X
Out[37]:
In [38]:
plt.plot(X[:,0], X[:,1], "*")
plt.axis('equal')
Out[38]:
In [39]:
pca = PCA(n_components=2)
In [40]:
pca.fit(X)
Out[40]:
In [41]:
print(pca.explained_variance_ratio_)
In [42]:
print(pca.singular_values_)
In [43]:
X2 = pca.transform(X)
X2
Out[43]:
In [44]:
plt.plot(X2[:,0], X2[:,1], "*")
#plt.axis('equal')
Out[44]:
In [45]:
pca_full = PCA(n_components=2, svd_solver='full')
pca_full.fit( X )
Out[45]:
In [46]:
pca_arpack = PCA(n_components=1, svd_solver='arpack') #долюно быть n_components < n_features
pca_arpack.fit( X )
Out[46]:
In [47]:
X_arpack=pca_arpack.transform(X)
X_arpack
Out[47]:
In [48]:
plt.plot(np.sort(X_arpack[:,0]), "*")
Out[48]:
Ядровый Метод Главных Компонент (Kernel PCA)
Выделяемые компоненты совсем не обязательно должны быть "прямыми". Они могут быть и кривыми засчет подмены понятия скалярного произведения. Оно заменяется на так называемое ядро. В частности, данные могут находится на кривой поверхности (например, сфере) в признаковом пространстве. Усовершенствование за счет ядер позволит и в таком случае найти компоненты.In [49]:
from sklearn.decomposition import KernelPCA
from sklearn.datasets import make_circles
In [50]:
# Для обеспечения повтаряемости эксперимента
np.random.seed(0)
Синтезированные данные
In [51]:
# Создадим "кривые" данные, две концентрические окружности с центром в точке (0, 0)
X, y = make_circles(n_samples=400, factor=.3, noise=.05)
X.shape, y.shape
Out[51]:
In [52]:
np.unique(y)
Out[52]:
In [53]:
plt.title("Исходное пространство")
reds = y == 0
blues = y == 1
plt.scatter(X[reds, 0], X[reds, 1], c="red",
s=20, edgecolor='k')
plt.scatter(X[blues, 0], X[blues, 1], c="blue",
s=20, edgecolor='k')
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
Out[53]:
In [54]:
pca = PCA(n_components=2)
X_pca = pca.fit_transform(X)
In [55]:
plt.scatter(X_pca[reds, 0], X_pca[reds, 1], c="red",
s=20, edgecolor='k')
plt.scatter(X_pca[blues, 0], X_pca[blues, 1], c="blue",
s=20, edgecolor='k')
plt.title("После МГК преобразования")
plt.xlabel("1ая главная компонента")
plt.ylabel("2ая главная компонента")
Out[55]:
В общем ничего не поменялось. Теперь пробуем ядровый метод главных компонент.
In [56]:
# Ядро rbf, radial basis function
kpca = KernelPCA(n_components=2, kernel="rbf", fit_inverse_transform=True, gamma=10)
X_kpca = kpca.fit_transform(X)
X_back = kpca.inverse_transform(X_kpca)
In [57]:
plt.scatter(X_kpca[reds, 0], X_kpca[reds, 1], c="red",
s=20, edgecolor='orange')
plt.scatter(X_kpca[blues, 0], X_kpca[blues, 1], c="blue",
s=20, edgecolor='c')
plt.title("После Ядровый МГК преобразования")
plt.xlabel("1ая главная компонента")
plt.ylabel("2ая главная компонента")
Out[57]:
Теперь победа. Исходные два множества четко разделились.
Начертим кривые координаты.
In [58]:
l = np.linspace(-1.5, 1.5, 5)
l
Out[58]:
In [59]:
grid = np.meshgrid( l, l )
grid
Out[59]:
In [60]:
grid[0][1,3],grid[1][1,3]
Out[60]:
In [61]:
type(grid), len(grid) #! Это список, а не массив !
Out[61]:
In [62]:
grid[0]
Out[62]:
In [63]:
grid[1]
Out[63]:
In [64]:
np.ravel([[1, 2],[3,4]]) # выпрямляет массив в одномерный, схоже с flatten
Out[64]:
In [65]:
plt.title("Кривые координаты в исходном пространстве")
reds = y == 0
blues = y == 1
plt.scatter(X[reds, 0], X[reds, 1], c="red",
s=20, edgecolor='k')
plt.scatter(X[blues, 0], X[blues, 1], c="blue",
s=20, edgecolor='k')
X1, X2 = np.meshgrid(np.linspace(-1.0, 1.0, 50), np.linspace(-1.0, 1.0, 50))
X_grid = np.array([np.ravel(X1), np.ravel(X2)]).T
# projection on the first principal component (in the phi space)
Z_grid = kpca.transform(X_grid)[:, 0].reshape(X1.shape)
plt.contour(X1, X2, Z_grid, colors='grey', linewidths=1, origin='lower')
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
Out[65]:
In [66]:
plt.scatter(X_back[reds, 0], X_back[reds, 1], c="red",
s=20, edgecolor='k')
plt.scatter(X_back[blues, 0], X_back[blues, 1], c="blue",
s=20, edgecolor='k')
plt.title("Обратно в исходном пространстве")
plt.xlabel("$x_1$")
plt.ylabel("$x_2$")
Out[66]:
Метод главных компонент на примере Iris dataset
In [67]:
# Загрузка данных из самой бибилотеки.
from sklearn import datasets # В dataset есть нексолько традиционых наборов.
In [68]:
iris = datasets.load_iris()
X = iris.data # Сами данные
y = iris.target # Класс, в данной заметке он нужен нам только для визуализации.
target_names = iris.target_names # Имена классов.
X.shape, y.shape, target_names.shape # Покажем чему раны размеры.
Out[68]:
Инача говоря данные содержат 150 примеров. Для каждого примера заданы 4 признака, а также разметку (чем пример является). Для числовых обозначений классов даны соответствующие текстовые названия.
In [69]:
X # Данные, четыре признака. 150 примера (прецендента).
Out[69]:
In [70]:
y # Класса. Фактически сначала идут преценденты относящиеся к первому классу, потом второму, а уже потом третьему.
Out[70]:
In [71]:
target_names # Имена классов, т.е. 0 классу соответствует имя setosa, ..., 2 -- virginica
# Да, тип элементов может быть текстом!
Out[71]:
In [72]:
np.unique(y) # Найдем номера различающихся классов. Функция ищет уникальные элементы массива.
Out[72]:
In [73]:
# Напомню.
y * 2 + 1 # Выполняет данное вычисление поэлементно.
Out[73]:
In [74]:
# По аналогии.
y==1 # Применили проверку на равенство 1 для каждого из элементов массива. Истина -- True, ложь -- False.
Out[74]:
In [75]:
# Напечатам размер каждого из классов, т.е. количество прецендентов в каждом классе.
# Далее идет так называемый цикл, конкретно цикл for. Он выполняет некий код для каждого из элементов по одному.
for i in np.unique(y): # i пробегает список элементов [0,1,2], т.е. принимает каждое из данных значений по отдельности.
# Тело начинается с отступа. Замечу, в предыдущей строчке в конце был символ :!
print( "размер класса ", i, sum( y==i ) ) # Для каждого из жлементов выполянется данный код. Идет печать.
# В арифметических операциях истина (True) равна 1, а лож (False).
# Видим, что классы имеют один и тот же размер, равный 50.
In [76]:
nums = np.unique(y) # Фактически это [0, 1, 2]
list( zip( nums, target_names ) ) # Списки сшиваются по-элементно вместе. Как молния с бегунком у куртки.
# nums[0] и target_names[0] в первый элемент нового списка... (0, 'setosa') ... и так далее
# по аналогии все далается при большем количестве списоков. Рзамеры их должны сопадать!
Out[76]:
In [77]:
X_1 = X[ y == 1 ] # Выбираем только те элементы для которых условие истинно, т.е. относящиеся к классу 1.
X_1
Out[77]:
In [78]:
X_1.shape, sum(y==1) # Перепроверяем размеры.
Out[78]:
In [79]:
X[y == 1, 0] # Берем массив нулевых признаков отсносящихся к классу 1. Фактически 0ой вектор-столбец из X_1.
# Иначе говоря, идет проекция всего признакового пространства (там было 4 признака) на 0 ось.
Out[79]:
In [80]:
colors = ['navy', 'turquoise', 'darkorange']
lw = 2
for color, i, target_name in zip(colors, [0, 1, 2], target_names): # Хитрым образом записаный цикл for.
plt.scatter(X[y == i, 0], X[y == i, 1], color=color, alpha=.8, lw=lw, # Отрисовываем каждый класс по отдельности.
label=target_name) # Названия наборов данных.
plt.legend(loc='best', shadow=False, scatterpoints = 1 ) # Отображение легенды.
plt.title('Исходные IRIS данные') # Название графика.
Out[80]:
Синий класс можно почти отделить, а вот оранжевый и берюзовый скорее нет. Как быть? Приобразуем систему координат так, что бы было видно. Обращаю внимание: исходные даные 4х мерные, а картинка 2х мерная.
Упр. Изучить другие проекции.
In [81]:
pca = PCA(n_components=2, svd_solver='full')
X_r = pca.fit_transform(X) # Сразу строим модель и преобразуем.
In [82]:
for color, i, target_name in zip(colors, [0, 1, 2], target_names):
plt.scatter(X_r[y == i, 0], X_r[y == i, 1], color=color, alpha=.8, lw=lw,
label=target_name)
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('PCA of IRIS dataset')
Out[82]:
Теперь данные гораздо проще отделить данные друг отдруга. По крайней мере можно провести четкую границу/ломанную между классами.
Посмотрим на все это дело в 3d
In [83]:
from mpl_toolkits.mplot3d import Axes3D
from sklearn import decomposition
from sklearn import datasets
np.random.seed(5)
In [84]:
iris = datasets.load_iris()
X = iris.data
y = iris.target
In [85]:
np.choose([0, 2, 1, 2], ['a','b','c']) # Один список используется для индексации элементов из второго списка
Out[85]:
In [86]:
fig = plt.figure(figsize=(8, 6))
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=20, azim=60)
for name, label in [('Setosa', 0), ('Versicolour', 1), ('Virginica', 2)]:
ax.text3D(X[y == label, 0].mean(),
X[y == label, 1].mean() + 1.0,
X[y == label, 2].mean(), name,
horizontalalignment='center')
# Отрисовываем данные.
yy = np.choose(y, [1, 2, 0]).astype(np.float)
ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=yy, cmap=plt.cm.spring, edgecolor='k')
Out[86]:
In [87]:
pca = decomposition.PCA(n_components=3)
pca.fit(X)
Xpca = pca.transform(X)
fig = plt.figure(figsize=(8, 6))
ax = Axes3D(fig, rect=[0, 0, .95, 1], elev=20, azim=90)
for name, label in [('Setosa', 0), ('Versicolour', 1), ('Virginica', 2)]:
ax.text3D(Xpca[y == label, 0].mean(),
Xpca[y == label, 1].mean() + 1.0,
Xpca[y == label, 2].mean(), name,
horizontalalignment='center')
# Отрисовываем данные.
yy = np.choose(y, [1, 2, 0]).astype(np.float)
ax.scatter(Xpca[:, 0], Xpca[:, 1], Xpca[:, 2], c=yy, cmap=plt.cm.spring, edgecolor='k')
Out[87]:
Ослабим часть условий на метод главных компонент
In [88]:
# В Методе главных компонент считалось, что разброс энергии (шума) постоянен по всем осям.
# В данном методе каждая из ортогональных осей независима.
from sklearn.decomposition import FactorAnalysis
In [89]:
fa = FactorAnalysis( n_components=2 )
X_fa = fa.fit_transform(X)
In [90]:
for color, i, target_name in zip(colors, [0, 1, 2], target_names):
plt.scatter(X_fa[y == i, 0], X_fa[y == i, 1], color=color, alpha=.8, lw=lw,
label=target_name)
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('FA of IRIS dataset')
Out[90]:
Когда будет разница видна?
В методе главных компонентов оси ортонормированы. Это жесткое требование. от него можно отказатся применив следующий подход.
Еще более общий случай метода главных компонент
In [91]:
# Теперь оси необязательно ортогональны
from sklearn.decomposition import FastICA, FactorAnalysis
In [92]:
# ?rng
In [93]:
# Создадим данные
rng = np.random.RandomState(42)
S = rng.standard_t(1.5, size=(20000, 2)) #standard_t 1.5, size=(20000, 2) randn
S[:, 0] *= 2.
S, S.shape
Out[93]:
In [94]:
SS = S / S.std()
plt.scatter(SS[:, 0], SS[:, 1], s=2, marker='o', zorder=10,
color='steelblue', alpha=0.5)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
Out[94]:
In [95]:
# Смешиваем данные распределенные по разным осям
A = np.array([[1, 1], [0, 2]]) # Матрица смешивания
X = np.dot(S, A.T) # Формируем наблюдения
XX = X / X.std()
plt.scatter(XX[:, 0], XX[:, 1], s=2, marker='o', zorder=10,
color='steelblue', alpha=0.5)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
Out[95]:
In [96]:
pca = FactorAnalysis() #PCA()
S_pca_ = pca.fit(X).transform(X)
ica = FastICA(random_state=rng)
S_ica_ = ica.fit(X).transform(X) # Estimate the sources
S_ica_ /= S_ica_.std(axis=0)
In [97]:
def plot_samples(S, axis_list=None):
plt.scatter(S[:, 0], S[:, 1], s=2, marker='o', zorder=10,
color='steelblue', alpha=0.5)
if axis_list is not None:
colors = ['orange', 'red']
for color, axis in zip(colors, axis_list):
axis /= axis.std()
x_axis, y_axis = axis
# Trick to get legend to work
plt.plot(0.1 * x_axis, 0.1 * y_axis, linewidth=2, color=color)
plt.quiver(0, 0, x_axis, y_axis, zorder=11, width=0.01, scale=6,
color=color)
plt.hlines(0, -3, 3)
plt.vlines(0, -3, 3)
plt.xlim(-3, 3)
plt.ylim(-3, 3)
plt.xlabel('x')
plt.ylabel('y')
plt.figure(figsize=(8, 6))
plt.subplot(2, 2, 1)
plot_samples(S / S.std())
plt.title('True Independent Sources')
axis_list = [pca.components_.T, ica.mixing_]
plt.subplot(2, 2, 2)
plot_samples(X / np.std(X), axis_list=axis_list)
legend = plt.legend(['PCA', 'ICA'], loc='upper right')
legend.set_zorder(100)
plt.title('Observations')
plt.subplot(2, 2, 3)
plot_samples(S_pca_ / np.std(S_pca_, axis=0))
plt.title('PCA recovered signals')
plt.subplot(2, 2, 4)
plot_samples(S_ica_ / np.std(S_ica_))
plt.title('ICA recovered signals')
Out[97]:
In [ ]: