In [1]:
import numpy as np
import pandas as pd
import matplotlib
from matplotlib import pyplot as plt
import matplotlib.patches as mpatches
matplotlib.style.use('ggplot')
%matplotlib inline
Любой набор данных представляет собой матрицу $X$.
Метод главных компонент последовательно находит следующие линейные комбинации признаков (компоненты) из $X$:
Предположения, в рамках которых данных подход будет работать хорошо:
Как это выглядит математически?
Обозначим следующим образом выборочную матрицу ковариации данных: $\hat{C} \propto Q = X^TX$. ($Q$ отличается от $\hat{C}$ нормировкой на число объектов).
Сингулярное разложение матрицы $Q$ выглядит следующим образом:
$$Q = X^TX = W \Lambda W^T$$Можно строго показать, что столбцы матрицы $W$ являются главными компонентами матрицы $X$, т.е. комбинациями признаков, удовлетворяющих двум условиям, указанным в начале. При этом дисперсия данных вдоль направления, заданного каждой компонентой, равна соответствующему значению диагональной матрицы $\Lambda$.
Как же на основании этого преобразования производить уменьшение размерности? Мы можем отранжировать компоненты, используя значения дисперсий данных вдоль них.
Сделаем это: $\lambda_{(1)} > \lambda_{(2)} > \dots > \lambda_{(D)}$.
Тогда, если мы выберем компоненты, соответствующие первым $d$ дисперсиям из этого списка, мы получим набор из $d$ новых признаков, которые наилучшим образом описывают дисперсию изначального набора данных среди всех других возможных линейных комбинаций исходных признаков матрицы $X$.
Получается, что метод главных компонент позволяет нам ранжировать полученные компоненты по "значимости", а также запустить процесс их отбора.
In [2]:
from sklearn.decomposition import PCA
mu = np.zeros(2)
C = np.array([[3,1],[1,2]])
data = np.random.multivariate_normal(mu, C, size=50)
plt.scatter(data[:,0], data[:,1])
plt.show()
Путём диагонализации истинной матрицы ковариаций $C$, мы можем найти преобразование исходного набора данных, компоненты которого наилучшим образом будут описывать дисперсию, с учётом их ортогональности друг другу:
In [3]:
v, W_true = np.linalg.eig(C)
plt.scatter(data[:,0], data[:,1])
# построим истинные компоненты, вдоль которых максимальна дисперсия данных
plt.plot(data[:,0], (W_true[0,0]/W_true[0,1])*data[:,0], color="g")
plt.plot(data[:,0], (W_true[1,0]/W_true[1,1])*data[:,0], color="g")
g_patch = mpatches.Patch(color='g', label='True components')
plt.legend(handles=[g_patch])
plt.axis('equal')
limits = [np.minimum(np.amin(data[:,0]), np.amin(data[:,1])),
np.maximum(np.amax(data[:,0]), np.amax(data[:,1]))]
plt.xlim(limits[0],limits[1])
plt.ylim(limits[0],limits[1])
plt.draw()
А теперь сравним эти направления с направлениями, которые выбирает метод главных компонент:
In [4]:
def plot_principal_components(data, model, scatter=True, legend=True):
W_pca = model.components_
if scatter:
plt.scatter(data[:,0], data[:,1])
plt.plot(data[:,0], -(W_pca[0,0]/W_pca[0,1])*data[:,0], color="c")
plt.plot(data[:,0], -(W_pca[1,0]/W_pca[1,1])*data[:,0], color="c")
if legend:
c_patch = mpatches.Patch(color='c', label='Principal components')
plt.legend(handles=[c_patch], loc='lower right')
# сделаем графики красивыми:
plt.axis('equal')
limits = [np.minimum(np.amin(data[:,0]), np.amin(data[:,1]))-0.5,
np.maximum(np.amax(data[:,0]), np.amax(data[:,1]))+0.5]
plt.xlim(limits[0],limits[1])
plt.ylim(limits[0],limits[1])
plt.draw()
In [5]:
model = PCA(n_components=2)
model.fit(data)
plt.scatter(data[:,0], data[:,1])
# построим истинные компоненты, вдоль которых максимальна дисперсия данных
plt.plot(data[:,0], (W_true[0,0]/W_true[0,1])*data[:,0], color="g")
plt.plot(data[:,0], (W_true[1,0]/W_true[1,1])*data[:,0], color="g")
# построим компоненты, полученные с использованием метода PCA:
plot_principal_components(data, model, scatter=False, legend=False)
c_patch = mpatches.Patch(color='c', label='Principal components')
plt.legend(handles=[g_patch, c_patch])
plt.draw()
Видно, что уже при небольшом количестве данных они отличаются незначительно. Увеличим размер выборки:
In [6]:
data_large = np.random.multivariate_normal(mu, C, size=5000)
model = PCA(n_components=2)
model.fit(data_large)
plt.scatter(data_large[:,0], data_large[:,1], alpha=0.1)
# построим истинные компоненты, вдоль которых максимальна дисперсия данных
plt.plot(data_large[:,0], (W_true[0,0]/W_true[0,1])*data_large[:,0], color="g")
plt.plot(data_large[:,0], (W_true[1,0]/W_true[1,1])*data_large[:,0], color="g")
# построим компоненты, полученные с использованием метода PCA:
plot_principal_components(data_large, model, scatter=False, legend=False)
c_patch = mpatches.Patch(color='c', label='Principal components')
plt.legend(handles=[g_patch, c_patch])
plt.draw()
В этом случае главные компоненты значительно точнее приближают истинные направления данных, вдоль которых наблюдается наибольшая дисперсия.
Как формализовать предположения метода, указанные выше? При помощи вероятностной модели!
Задача, стоящая за любым методом уменьшения размерности: получить из набора зашумлённых признаков $X$ истинные значения $Y$, которые на самом деле определяют набор данных (т.е. сведение датасета с большим количеством признаков к данным, имеющим т.н. "эффективную размерность").
В случае метода главных компонент мы хотим найти направления, вдоль которых максимальна дисперсия, с учётом описанных выше предположений о структуре данных и компонент.
Материал, описанный ниже в данной секции, не обязателен для ознакомления для выполнения следующего задания, т.к. требует некоторых знаний статистики.
Для тех, кто собирается его пропустить: в конце раздела мы получим метрику качества, которая должна определять, насколько данные хорошо описываются построенной моделью при заданном числе компонент. Отбор признаков при этом сводится к тому, что мы выбираем то количеством компонент, при котором используемая метрика (логарифм правдоподобия) является максимальной.
С учётом предположений задача метода главных компонент выглядит следующим образом:
$$ x = Wy + \mu + \epsilon$$где:
Исходя из распределения шума, выпишем распределение на $x$:
$$p(x \mid y) = \mathcal{N}(Wx + \mu, \sigma^2I) $$Введём априорное распределение на $y$:
$$p(y) = \mathcal{N}(0, 1)$$Выведем из этого при помощи формулы Байеса маргинальное распределение на $p(x)$:
$$p(x) = \mathcal{N}(\mu, \sigma^2I + WW^T)$$Тогда правдоподобие набора данных при условии используемой модели выглядит следующим образом:
$$\mathcal{L} = \sum_{i=1}^N \log p(x_i) = -N/2 \Big( d\log(2\pi) + \log |C| + \text{tr}(C^{-1}S) \Big)$$где:
Значение $\mathcal{L}$ имеет смысл логарифма вероятности пронаблюдать данные $X$ при условии, что они удовлетворяют предположениям модели метода главных компонент. Чем оно больше -- тем лучше модель описывает наблюдаемые данные.
Рассмотрим набор данных размерности $D$, чья реальная размерность значительно меньше наблюдаемой (назовём её $d$). От вас требуется:
Для оценки логарифма правдоподобия модели для заданного числа главных компонент при помощи метода кросс-валидации используйте следующие функции:
model = PCA(n_components=n)
scores = cv_score(model, data)
Обратите внимание, что scores -- это вектор, длина которого равна числу фолдов. Для получения оценки на правдоподобие модели его значения требуется усреднить.
Для визуализации оценок можете использовать следующую функцию:
plot_scores(d_scores)
которой на вход передаётся вектор полученных оценок логарифма правдоподобия данных для каждого $\hat{d}$.
Для интересующихся: данные для заданий 1 и 2 были сгенерированны в соответствии с предполагаемой PCA моделью. То есть: данные $Y$ с эффективной размерностью $d$, полученные из независимых равномерных распределений, линейно траснформированны случайной матрицей $W$ в пространство размерностью $D$, после чего ко всем признакам был добавлен независимый нормальный шум с дисперсией $\sigma$.
In [14]:
from sklearn.decomposition import PCA
from sklearn.cross_validation import cross_val_score as cv_score
def plot_scores(d_scores):
n_components = np.arange(1,len(d_scores)+1)
plt.plot(n_components, d_scores, 'b', label='PCA scores')
plt.xlim(n_components[0], n_components[-1])
plt.xlabel('n components')
plt.ylabel('cv scores')
plt.legend(loc='lower right')
plt.show()
def write_answer_1(optimal_d):
with open("pca_answer1.txt", "w") as fout:
fout.write(str(optimal_d))
data = pd.read_csv('data_task1.csv')
d_scores = []
# place your code here
D = len(data.columns) + 1
for compon in range(1, D):
model = PCA(n_components=compon)
scores = cv_score(model, data)
d_scores.append(scores.mean())
plot_scores(np.array(d_scores))
Мы знаем, что каждой главной компоненте соответствует описываемай ей дисперсия данных (дисперсия данных при проекции на эту компоненту). Она численно равна значению диагональных элементов матрицы $\Lambda$, получаемой из спектрального разложения матрицы ковариации данных (смотри теорию выше).
Исходя из этого, мы можем отсортировать дисперсию данных вдоль этих компонент по убыванию, и уменьшить размерность данных, отбросив $q$ итоговых главных компонент, имеющих наименьшую дисперсию.
Делать это разными двумя способами. Например, если вы вдальнейшем обучаете на данных с уменьшенной размерностью модель классификации или регрессии, то можно запустить итерационный процесс: удалять компоненты с наименьшей дисперсией по одной, пока качество итоговой модели не станет значительно хуже.
Более общий способ отбора признаков заключается в том, что вы можете посмотреть на разности в дисперсиях в отсортированном ряде $\lambda_{(1)} > \lambda_{(2)} > \dots > \lambda_{(D)}$: $\lambda_{(1)}-\lambda_{(2)}, \dots, \lambda_{(D-1)} - \lambda_{(D)}$, и удалить те компоненты, на которых разность будет наибольшей. Именно этим методом вам и предлагается воспользоваться для тестового набора данных.
Рассмотрим ещё один набор данных размерности $D$, чья реальная размерность значительно меньше наблюдаемой (назовём её также $d$). От вас требуется:
Для построения модели PCA используйте функцию:
model.fit(data)
Для трансформации данных используйте метод:
model.transform(data)
Оценку дисперсий на трансформированных данных от вас потребуется реализовать вручную. Для построения графиков можно воспользоваться функцией
plot_variances(d_variances)
которой следует передать на вход отсортированный по убыванию вектор дисперсий вдоль компонент.
In [21]:
from sklearn.decomposition import PCA
from sklearn.cross_validation import cross_val_score as cv_score
def plot_variances(d_variances):
n_components = np.arange(1,d_variances.size+1)
plt.plot(n_components, d_variances, 'b', label='Component variances')
plt.xlim(n_components[0], n_components[-1])
plt.xlabel('n components')
plt.ylabel('variance')
plt.legend(loc='upper right')
plt.show()
def write_answer_2(optimal_d):
with open("pca_answer2.txt", "w") as fout:
fout.write(str(optimal_d))
data = pd.read_csv('data_task2.csv')
# place your code here
model = PCA(n_components=D)
model.fit(data)
X = model.transform(data)
disp = np.var(X, axis=0)
disp = np.sort(disp)[::-1]
print disp.shape, disp
# data.shape
# disp.sort()
# disp.sort()
dispdiff = np.append(disp[1:], [0])
print dispdiff
# model.explained_variance_ratio_
diff = (disp - dispdiff)[:-1]
# np.array([1,2,3,4])[:-1]
np.argmax(diff)
plot_variances(diff)
# diff
ans = np.argmax(diff) + 1
print ans
В качестве главных компонент мы получаем линейные комбинации исходных призанков, поэтому резонно возникает вопрос об их интерпретации.
Для этого существует несколько подходов, мы рассмотрим два:
Первый способ подходит в том случае, когда все объекты из набора данных не несут для нас никакой семантической информации, которая уже не запечатлена в наборе признаков.
Второй способ подходит для случая, когда данные имеют более сложную структуру. Например, лица для человека несут больший семантический смысл, чем вектор значений пикселей, которые анализирует PCA.
Рассмотрим подробнее способ 1: он заключается в подсчёте коэффициентов корреляций между исходными признаками и набором главных компонент.
Так как метод главных компонент является линейным, то предлагается для анализа использовать корреляцию Пирсона, выборочный аналог которой имеет следующую формулу:
$$r_{jk} = \frac{\sum_{i=1}^N (x_{ij} - \bar{x}_j) (y_{ik} - \bar{y}_k)}{\sum_{i=1}^N (x_{ij} - \bar{x}_j)^2 \sum_{i=1}^N (y_{ik} - \bar{y}_k)^2} $$где:
Корреляция Пирсона является мерой линейной зависимости. Она равна 0 в случае, когда величины независимы, и $\pm 1$, если они линейно зависимы. Исходя из степени корреляции новой компоненты с исходными признаками, можно строить её семантическую интерпретацию, т.к. смысл исходных признаков мы знаем.
Набор данных состоит из 4 признаков, посчитанных для 150 ирисов. Каждый из них принадлежит одному из трёх видов. Визуализацию проекции данного датасета на две компоненты, которые описывают наибольшую дисперсию данных, можно получить при помощи функции
plot_iris(transformed_data, target, target_names)
на вход которой требуется передать данные, преобразованные при помощи PCA, а также информацию о классах. Цвет точек отвечает одному из трёх видов ириса.
Для того чтобы получить имена исходных признаков, используйте следующий список:
iris.feature_names
При подсчёте корреляций не забудьте центрировать признаки и проекции на главные компоненты (вычитать из них среднее).
In [32]:
from sklearn import datasets
def plot_iris(transformed_data, target, target_names):
plt.figure()
for c, i, target_name in zip("rgb", [0, 1, 2], target_names):
plt.scatter(transformed_data[target == i, 0],
transformed_data[target == i, 1], c=c, label=target_name)
plt.legend()
plt.show()
def write_answer_3(list_pc1, list_pc2):
with open("pca_answer3.txt", "w") as fout:
fout.write(" ".join([str(num) for num in list_pc1]))
fout.write(" ")
fout.write(" ".join([str(num) for num in list_pc2]))
# загрузим датасет iris
iris = datasets.load_iris()
data = iris.data
target = iris.target
target_names = iris.target_names
# place your code here
# загрузим датасет iris
iris = datasets.load_iris()
data = iris.data
target = iris.target
target_names = iris.target_names
# place your code here
model = PCA(n_components=4)
X = model.fit_transform(data)
Xtwo = X[:,:2]
from scipy.stats import pearsonr
print pearsonr(data[:,3],Xtwo[:,0])
corr1 = [1, 3, 4]
corr2 = [2]
print corr1, corr2
write_answer_3(corr1, corr2)
plot_iris(Xtwo, target, target_names)
Рассмотрим теперь величину, которую можно проинтерпретировать, как квадрат косинуса угла между объектом выборки и главной компонентой:
$$ cos^2_{ik} = \frac{(f_{ik} - \bar{f}_{k})^2}{\sum_{k=1}^d (f_{ik} - \bar{f}_{k})^2} $$где
Очевидно, что
$$ \sum_{k=1}^d cos^2_{ik} = 1 $$Это значит, что для каждого объекта мы в виде данной величины получили веса, пропорциональные вкладу, которую вносит данный объект в дисперсию каждой компоненты. Чем больше вклад, тем более значим объект для описания конкретной главной компоненты.
Для каждой компоненты найдите и визуализируйте лицо, которое вносит наибольший относительный вклад в неё. Для визуализации используйте функцию
plt.imshow(image.reshape(image_shape))Передайте в функцию write_answer_4 список номеров лиц с наибольшим относительным вкладом в дисперсию каждой из компонент, список начинается с 0.
In [36]:
from sklearn.datasets import fetch_olivetti_faces
from sklearn.decomposition import RandomizedPCA
def write_answer_4(list_pc):
with open("pca_answer4.txt", "w") as fout:
fout.write(" ".join([str(num) for num in list_pc]))
data = fetch_olivetti_faces(shuffle=True, random_state=0).data
image_shape = (64, 64)
In [38]:
model = RandomizedPCA(n_components=10)
X = model.fit_transform(data)
# fmean = X[:,0].mean(axis=0)
fmean = X.mean(axis=0)
# fmean
# plt.imshow(data[3].reshape(image_shape))
disp = []
for i in range(len(X)):
divis = np.sum((X[i] - fmean)**2)
cosin = (X[i] - fmean)**2 / divis
disp.append(cosin)
disp = np.array(disp)
ans = []
for i in range(10):
ans.append(np.argmax(disp[:,i]))
print ans
write_answer_4(ans)
Рассмотренные выше задачи являются, безусловно, модельными, потому что данные для них были сгенерированы в соответствии с предположениями метода главных компонент. На практике эти предположения, естественно, выполняются далеко не всегда. Рассмотрим типичные ошибки PCA, которые следует иметь в виду перед тем, как его применять.
In [30]:
C1 = np.array([[10,0],[0,0.5]])
phi = np.pi/3
C2 = np.dot(C1, np.array([[np.cos(phi), np.sin(phi)],
[-np.sin(phi),np.cos(phi)]]))
data = np.vstack([np.random.multivariate_normal(mu, C1, size=50),
np.random.multivariate_normal(mu, C2, size=50)])
plt.scatter(data[:,0], data[:,1])
# построим истинные интересующие нас компоненты
plt.plot(data[:,0], np.zeros(data[:,0].size), color="g")
plt.plot(data[:,0], 3**0.5*data[:,0], color="g")
# обучим модель pca и построим главные компоненты
model = PCA(n_components=2)
model.fit(data)
plot_principal_components(data, model, scatter=False, legend=False)
c_patch = mpatches.Patch(color='c', label='Principal components')
plt.legend(handles=[g_patch, c_patch])
plt.draw()
В чём проблема, почему pca здесь работает плохо? Ответ прост: интересующие нас компоненты в данных коррелированны между собой (или неортогональны, в зависимости от того, какой терминологией пользоваться). Для поиска подобных преобразований требуются более сложные методы, которые уже выходят за рамки метода главных компонент.
Для интересующихся: то, что можно применить непосредственно к выходу метода главных компонент, для получения подобных неортогональных преобразований, называется методами ротации. Почитать о них можно в связи с другим методом уменьшения размерности, который называется Factor Analysis (FA), но ничего не мешает их применять и к главным компонентам.
In [31]:
C = np.array([[0.5,0],[0,10]])
mu1 = np.array([-2,0])
mu2 = np.array([2,0])
data = np.vstack([np.random.multivariate_normal(mu1, C, size=50),
np.random.multivariate_normal(mu2, C, size=50)])
plt.scatter(data[:,0], data[:,1])
# обучим модель pca и построим главные компоненты
model = PCA(n_components=2)
model.fit(data)
plot_principal_components(data, model)
plt.draw()
Очевидно, что в данном случае метод главных компонент будет считать вертикальную компоненту более значимой для описания набора данных, чем горизонтальную.
Но, например, в случае, когда данные из левого и правого кластера относятся к разным классам, для их линейной разделимости вертикальная компонента является шумовой. Несмотря на это, её метод главных компонент никогда шумовой не признает, и есть вероятность, что отбор признаков с его помощью выкинет из ваших данных значимые для решаемой вами задачи компоненты просто потому, что вдоль них значения имеют низкую дисперсию.
Справляться с такими ситуациями могут некоторые другие методы уменьшения размерности данных, например, метод независимых компонент (Independent Component Analysis, ICA).