Метод главных компонент или компонентный анализ (principal component analysis, PCA) - один из важнейших методов в арсенале зоолога или эколога. К сожалению, в тех случаях, когда вполне уместным является применение компонентного анализа, сплошь и рядом применяют кластерный анализ.
Типичная задача, для которой полезен компонентный анализ, такова: есть некое множество объектов, каждый из которых охарактеризован по определенному (достаточно большому) количеству признаков. Исследователя интересуют закономерности, отраженные в разнообразии этих объектов. В том случае, когда есть основания предполагать, что объекты распределены по иерархически соподчиненным группам, можно использовать кластерный анализ - метод классификации (распределения по группам). Если нет оснований ожидать, что в разнообразии объектов отражена какая-то иерархия, логично использовать ординацию (упорядоченное расположение). Если каждый объект охарактеризован по достаточно большому количеству признаков (по крайней мере - такому количеству признаков, какое не получается адекватно отразить на одном графике), оптимально начинать исследование данных с анализа главных компонент. Дело в том, что этот метод является одновременно методом понижения размерности (количества измерений) данных.
Если группа рассматриваемых объектов охарактеризована значениями одного признака, для характеристики их разнообразия можно использовать гистограмму (для непрерывных признаков) или столбчатую диаграмму (для характеристики частот дискретного признака). Если объекты охарактеризованы двумя признаками, можно использовать двумерный график рассеяния, если тремя - трехмерный. А если признаков много? Можно попытаться на двумерном графике отразить взаимное расположение объектов друг относительно друга в многомерном пространстве. Обычно такое понижение размерности связано с потерей информации. Из разных возможных способов такого отображения надо выбрать тот, при котором потеря информации будет минимальной.
Поясним сказанное на самом простом примере: переходе от двумерного пространства к одномерному. Минимальное количество точек, которое задает двумерное пространство (плоскость) - 3. На рис. 9.1.1 показано расположение трех точек на плоскости. Координаты этих точек легко читаются по самому рисунку. Как выбрать прямую, которая будет нести максимальную информацию о взаиморасположении точек?
Рис. 9.1.1. Три точки на плоскости, заданной двумя признаками. На какую прямую будет проецироваться максимальная дисперсия этих точек?
Рассмотрим проекции точек на прямую A (показанную синим цветом). Координаты проекций этих точек на прямую A таковы: 2, 8, 10. Среднее значение - 6 2 / 3 . Дисперсия (2-6 2 / 3)+ (8-6 2 / 3)+ (10-6 2 / 3)=34 2 / 3 .
Теперь рассмотрим прямую B (показанную зеленым цветом). Координаты точек - 2, 3, 7; среднее значение - 4, дисперсия - 14. Таким образом, на прямую B отражается меньшая доля дисперсии, чем на прямую A.
Какова эта доля? Поскольку прямые A и B ортогональны (перпендикулярны), доли общей дисперсии, проецирующиеся на A и B, не пересекаются. Значит, общую дисперсию расположения интересующих нас точек можно вычислить как сумму этих двух слагаемых: 34 2 / 3 +14=48 2 / 3 . При этом на прямую A проецируется 71,2% общей дисперсии, а на прямую B - 28,8%.
А как определить, на какую прямую отразится максимальная доля дисперсии? Эта прямая будет соответствовать линии регрессии для интересующих нас точек, которая обозначена как C (красный цвет). На эту прямую отразится 77,2% общей дисперсии, и это - максимально возможное значение при данном расположении точек. Такую прямую, на которую проецируется максимальная доля общей дисперсии, называют первой главной компонентой .
А на какую прямую отразить оставшиеся 22,8% общей дисперсии? На прямую, перпендикулярную первой главной компоненте. Эта прямая тоже будет являться главной компонентой, ведь на нее отразится максимально возможная доля дисперсии (естественно, без учета той, которая отразилась на первую главную компоненту). Таким образом, это - вторая главная компонента .
Вычислив эти главные компоненты с помощью Statistica (диалог мы опишем чуть позже), мы получим картину, показанную на рис. 9.1.2. Координаты точек на главных компонентах показываются в стандартных отклонениях.
Рис. 9.1.2. Расположение трех точек, показанных на рис. 9.1.1, на плоскости двух главных компонент. Почему эти точки располагаются друг относительно друга иначе, чем на рис. 9.1.1?
На рис. 9.1.2 взаиморасположение точек оказывается измененным. Чтобы в дальнейшем правильно интерпретировать подобные картинки, следует рассмотреть причины отличий в расположении точек на рис. 9.1.1 и 9.1.2 подробнее. Точка 1 в обоих случаях находится правее (имеет большую координату по первому признаку и первой главной компоненте), чем точка 2. Но, почему-то, точка 3 на исходном расположении находится ниже двух других точек (имеет наименьшее значение признака 2), и выше двух других точек на плоскости главных компонент (имеет большую координату по второй компоненте). Это связано с тем, что метод главных компонент оптимизирует именно дисперсию исходных данных, проецирующихся на выбираемые им оси. Если главная компонента коррелирована с какой-то исходной осью, компонента и ось могут быть направлены в одну сторону (иметь положительную корреляцию) или в противоположные стороны (иметь отрицательные корреляции). Оба эти варианта равнозначны. Алгоритм метода главных компонент может «перевернуть» или не «перевернуть» любую плоскость; никаких выводов на основании этого делать не следует.
Однако точки на рис. 9.1.2 не просто «перевернуты» по сравнению с их взаиморасположением на рис. 9.1.1; определенным образом изменилось и их взаиморасположения. Отличия между точками по второй главной компоненте кажутся усиленными. 22,76% общей дисперсии, приходящиеся на вторую компоненту, «раздвинули» точки на такую же дистанцию, как и 77,24% дисперсии, приходящихся на первую главную компоненту.
Чтобы расположение точек на плоскости главных компонент соответствовало их действительному расположению, эту плоскость следовало бы исказить. На рис. 9.1.3. показаны два концентрических круга; их радиусы соотносятся как доли дисперсий, отражаемых первой и второй главными компонентами. Картинка, соответствующая рис. 9.1.2, искажена так, чтобы среднеквадратичное отклонение по первой главной компоненте соответствовало большему кругу, а по второй - меньшему.
Рис. 9.1.3. Мы учли, что на первую главную компоненту приходится бо льшая доля дисперсии, чем на вторую. Для этого мы исказили рис. 9.1.2, подогнав его под два концентрических круга, радиусы которых соотносятся, как доли дисперсий, приходящихся на главные компоненты. Но расположение точек все равно не соответствует исходному, показанному на рис. 9.1.1!
А почему взаимное расположение точек на рис. 9.1.3 не соответствует таковому на рис. 9.1.1? На исходном рисунке, рис. 9.1 точки расположены в соответствии со своими координатами, а не в соответствии с долями дисперсии, приходящимися на каждую ось. Расстоянию в 1 единицу по первому признаку (по оси абсцисс) на рис. 9.1.1 приходятся меньшая доля дисперсии точек по этой оси, чем расстоянию в 1 единицу по второму признаку (по оси ординат). А на рис 9.1.1 расстояния между точками определяются именно теми единицами, в которых измеряются признаки, по которым они описаны.
Несколько усложним задачу. В табл. 9.1.1 показаны координаты 10 точек в 10-мерном пространстве. Первые три точки и первые два измерения - это тот пример, который мы только что рассматривали.
Таблица 9.1.1. Координаты точек для дальнейшего анализа
|
Координаты |
||||||||||
В учебных целях вначале рассмотрим только часть данных из табл. 9.1.1. На рис. 9.1.4 мы видим положение десяти точек на плоскости первых двух признаков. Обратите внимание, что первая главная компонента (прямая C) прошла несколько иначе, чем в предыдущем случае. Ничего удивительного: на ее положение влияют все рассматриваемые точки.
Рис. 9.1.4. Мы увеличили количество точек. Первая главная компонента проходит уже несколько иначе, ведь на нее оказали влияние добавленные точки
На рис. 9.1.5 показано положение рассмотренных нами 10 точек на плоскости двух первых компонент. Обратите внимание: все изменилось, не только доля дисперсии, приходящейся на каждую главную компоненту, но даже положение первых трех точек!
Рис. 9.1.5. Ординация в плоскости первых главных компонент 10 точек, охарактеризованных в табл. 9.1.1. Рассматривались только значения двух первых признаков, последние 8 столбцов табл. 9.1.1 не использовались
В общем, это естественно: раз главные компоненты расположены иначе, то изменилось и взаиморасположение точек.
Трудности в сопоставлении расположения точек на плоскости главных компонент и на исходной плоскости значений их признаков могут вызвать недоумение: зачем использовать такой трудноинтерпретируемый метод? Ответ прост. В том случае, если сравниваемые объекты описаны всего по двум признакам, вполне можно использовать их ординацию по этим, исходным признакам. Все преимущества метода главных компонент проявляются в случае многомерных данных. Метод главных компонент в таком случае оказывается эффективным способом снижения размерности данных.
9.2. Переход к начальным данным с большим количеством измерений
Рассмотрим более сложный случай: проанализируем данные, представленные в табл. 9.1.1 по всем десяти признакам. На рис. 9.2.1 показано, как вызывается окно интересующего нас метода.
Рис. 9.2.1. Запуск метода главных компонент
Нас будет интересовать только выбор признаков для анализа, хотя диалог Statistica позмоляет намного более тонкую настройку (рис. 9.2.2).
Рис. 9.2.2. Выбор переменных для анализа
После выполнения анализа появляется окно его результатов с несколькими вкладками (рис. 9.2.3). Все основные окна доступны уже из первой вкладки.
Рис. 9.2.3. Первая вкладка диалога результатов анализа главных компонент
Можно увидеть, что анализ выделил 9 главных компонент, причем описал с их помощью 100% дисперсии, отраженной в 10 начальных признаках. Это означает, что один признак был лишним, избыточным.
Начнем просматривать результаты с кнопки «Plot case factor voordinates, 2D»: она покажет расположение точек на плоскости, заданной двумя главными компонентами. Нажав эту кнопку, мы попадем в диалог, где надо будет указать, какие мы будем использовать компоненты; естественно начинать анализ с первой и второй компонент. Результат - на рис. 9.2.4.
Рис. 9.2.4. Ординация рассматриваемых объектов на плоскости двух первых главных компонент
Положение точек изменилось, и это естественно: в анализ вовлечены новые признаки. На рис. 9.2.4 отражено более 65% всего разнообразия в положении точек друг относительно друга, и это уже нетривиальный результат. К примеру, вернувшись к табл. 9.1.1, можно убедиться в том, что точки 4 и 7, а также 8 и 10 действительно достаточно близки друг к другу. Впрочем, отличия между ними могут касаться других главных компонент, не показанных на рисунке: на них, все-таки, тоже приходится треть оставшейся изменчивости.
Кстати, при анализе размещения точек на плоскости главных компонент может возникнуть необходимость проанализировать расстояния между ними. Проще всего получить матрицу дистанций между точками с использованием модуля для кластерного анализа.
А как выделенные главные компоненты связаны с исходными признаками? Это можно узнать, нажав кнопку (рис. 9.2.3) Plot var. factor coordinates, 2D. Результат - на рис. 9.2.5.
Рис. 9.2.5. Проекции исходных признаков на плоскость двух первых главных компонент
Мы смотрим на плоскость двух главных компонент «сверху». Исходные признаки, которые никак не связаны с главными компонентами, будет перпендикулярны (или почти перпендикулярны) им и отразятся короткими отрезками, заканчивающимися вблизи начала координат. Так, меньше всего с двумя первыми главными компонентами связан признак № 6 (хотя он демонстрирует определенную положительную корреляцию с первой компонентой). Отрезки, соответствующие тем признакам, которые полностью отразятся на плоскости главных компонент, будут заканчиваться на охватывающей центр рисунка окружности единичного радиуса.
Например, можно увидеть, что на первую главную компоненту сильнее всего повлияли признаки 10 (связан положительной корреляцией), а также 7 и 8 (связаны отрицательной корреляцией). Чтобы рассмотреть структуру таких корреляций подробнее, можно нажать кнопку Factor coordinates of variables, и получить таблицу, показанную на рис. 9.2.6.
Рис. 9.2.6. Корреляции между исходными признаками и выделенными главными компонентами (Factors)
Кнопка Eigenvalues выводит величины, которые называются собственными значениями главных компонент . В верхней части окна, показанного на рис. 9.2.3, выведены такие значения для нескольких первых компонент; кнопка Scree plot показывает их в удобной для восприятия форме (рис. 9.2.7).
Рис. 9.2.7. Собственные значения выделенных главных компонент и доли отраженной ими общей дисперсии
Для начала надо понять, что именно показывает значение eigenvalue. Это - мера дисперсии, отразившейся на главную компоненту, измеренная в количестве дисперсии, приходившейся на каждый признак в начальных данных. Если eigenvalue первой главной компоненты равен 3,4, это означает, что на нее отражается больше дисперсии, чем на три признака из начального набора. Собственные величины линейно связаны с долей дисперсии, приходящейся на главную компоненту, единое что, сумма собственных значений равна количеству исходных признаков, а сумма долей дисперсии равна 100%.
А что означает, что информацию об изменчивости по 10 признакам удалось отразить в 9 главных компонентах? Что один из начальных признаков был избыточным, не добавлял никакой новой информации. Так и было; на рис. 9.2.8 показано, как был сгенерирован набор точек, отраженный в табл. 9.1.1.
В этой статье я бы хотел рассказать о том, как именно работает метод анализа главных компонент (PCA – principal component analysis) с точки зрения интуиции, стоящей за ее математическим аппаратом. Максимально просто, но подробно.
Математика вообще очень красивая и изящная наука, но порой ее красота скрывается за кучей слоев абстракции. Показать эту красоту лучше всего на простых примерах, которые, так сказать, можно покрутить, поиграть и пощупать, потому что в конце концов все оказывается гораздо проще, чем кажется на первый взгляд – самое главное понять и представить.
В анализе данных, как и в любом другом анализе, порой бывает нелишним создать упрощенную модель, максимально точно описывающую реальное положение дел. Часто бывает так, что признаки довольно сильно зависят друг от друга и их одновременное наличие избыточно.
К примеру, расход топлива у нас меряется в литрах на 100 км, а в США в милях на галлон. На первый взгляд, величиные разные, но на самом деле они строго зависят друг от друга. В миле 1600км, а в галлоне 3.8л. Один признак строго зависит от другого, зная один, знаем и другой.
Но гораздо чаще бывает так, что признаки зависят друг от друга не так строго и (что важно!) не так явно. Объем двигателя в целом положительно влияет на разгон до 100 км/ч, но это верно не всегда. А еще может оказаться, что с учетом не видимых на первый взгляд факторов (типа улучшения качества топлива, использования более легких материалов и прочих современных достижений), год автомобиля не сильно, но тоже влияет на его разгон.
Зная зависимости и их силу, мы можем выразить несколько признаков через один, слить воедино, так сказать, и работать уже с более простой моделью. Конечно, избежать потерь информации, скорее всего не удастся, но минимизировать ее нам поможет как раз метод PCA.
Выражаясь более строго, данный метод аппроксимирует n-размерное облако наблюдений до эллипсоида (тоже n-мерного), полуоси которого и будут являться будущими главными компонентами. И при проекции на такие оси (снижении размерности) сохраняется наибольшее количество информации.
Шаг 1. Подготовка данных
Здесь для простоты примера я не буду брать реальные обучающие датасеты на десятки признаков и сотни наблюдений, а сделаю свой, максимально простой игрушечный пример. 2 признака и 10 наблюдений будет вполне достаточно для описания того, что, а главное – зачем, происходит в недрах алгоритма.
Сгенерируем выборку:
X = np.arange(1,11) y = 2 * x + np.random.randn(10)*2 X = np.vstack((x,y)) print X OUT: [[ 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. ] [ 2.73446908 4.35122722 7.21132988 11.24872601 9.58103444 12.09865079 13.78706794 13.85301221 15.29003911 18.0998018 ]]
В данной выборке у нас имеются два признака, сильно коррелирующие друг с другом. С помощью алгоритма PCA мы сможем легко найти признак-комбинацию и, ценой части информации, выразить оба этих признака одним новым. Итак, давайте разбираться!
Для начала немного статистики. Вспомним, что для описания случайной величины используются моменты. Нужные нам – мат. ожидание и дисперсия. Можно сказать, что мат. ожидание – это «центр тяжести» величины, а дисперсия – это ее «размеры». Грубо говоря, мат. ожидание задает положение случайной величины, а дисперсия – ее размер.
Сам процесс проецирования на вектор никак не влияет на значения средних, так как для минимизации потерь информации наш вектор должен проходить через центр нашей выборки. Поэтому нет ничего страшного, если мы отцентрируем нашу выборку – линейно сдвинем ее так, чтобы средние значения признаков были равны 0. Это очень сильно упростит наши дальнейшие вычисления (хотя, стоит отметить, что можно обойтись и без центрирования).
Оператор, обратный сдвигу будет равен вектору изначальных средних значений – он понадобится для восстановления выборки в исходной размерности.
Xcentered = (X - x.mean(), X - y.mean()) m = (x.mean(), y.mean()) print Xcentered print "Mean vector: ", m OUT: (array([-4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5]), array([-8.44644233, -8.32845585, -4.93314426, -2.56723136, 1.01013247, 0.58413394, 1.86599939, 7.00558491, 4.21440647, 9.59501658])) Mean vector: (5.5, 10.314393916)
Дисперсия же сильно зависит от порядков значений случайной величины, т.е. чувствительна к масштабированию. Поэтому если единицы измерения признаков сильно различаются своими порядками, крайне рекомендуется стандартизировать их. В нашем случае значения не сильно разнятся в порядках, так что для простоты примера мы не будем выполнять эту операцию.
Шаг 2. Ковариационная матрица
В случае с многомерной случайной величиной (случайным вектором) положение центра все так же будет являться мат. ожиданиями ее проекций на оси. А вот для описания ее формы уже недостаточно толькое ее дисперсий по осям. Посмотрите на эти графики, у всех трех случайных величин одинаковые мат.ожидания и дисперсии, а их проекции на оси в целом окажутся одинаковы!
Для описания формы случайного вектора необходима ковариационная матрица.
Это матрица, у которой (i,j) -элемент является корреляцией признаков (X i , X j). Вспомним формулу ковариации:
В нашем случае она упрощается, так как E(X i) = E(X j) = 0:
Заметим, что когда X i = X j:
и это справедливо для любых случайных величин.
Таким образом, в нашей матрице по диагонали будут дисперсии признаков (т.к. i = j), а в остальных ячейках – ковариации соответствующих пар признаков. А в силу симметричности ковариации матрица тоже будет симметрична.
Замечание: Ковариационная матрица является обобщением дисперсии на случай многомерных случайных величин – она так же описывает форму (разброс) случайной величины, как и дисперсия.
И действительно, дисперсия одномерной случайной величины – это ковариационная матрица размера 1x1, в которой ее единственный член задан формулой Cov(X,X) = Var(X).
Итак, сформируем ковариационную матрицу Σ
для нашей выборки. Для этого посчитаем дисперсии X i и X j , а также их ковариацию. Можно воспользоваться вышенаписанной формулой, но раз уж мы вооружились Python’ом, то грех не воспользоваться функцией numpy.cov(X)
. Она принимает на вход список всех признаков случайной величины и возвращает ее ковариационную матрицу и где X – n-мерный случайный вектор (n-количество строк). Функция отлично подходит и для расчета несмещенной дисперсии, и для ковариации двух величин, и для составления ковариационной матрицы.
(Напомню, что в Python матрица представляется массивом-столбцом массивов-строк.)
Covmat = np.cov(Xcentered) print covmat, "n" print "Variance of X: ", np.cov(Xcentered) print "Variance of Y: ", np.cov(Xcentered) print "Covariance X and Y: ", np.cov(Xcentered) OUT: [[ 9.16666667 17.93002811] [ 17.93002811 37.26438587]] Variance of X: 9.16666666667 Variance of Y: 37.2643858743 Covariance X and Y: 17.9300281124
Шаг 3. Собственные вектора и значения (айгенпары)
О"кей, мы получили матрицу, описывающую форму нашей случайной величины, из которой мы можем получить ее размеры по x и y (т.е. X 1 и X 2), а также примерную форму на плоскости. Теперь надо надо найти такой вектор (в нашем случае только один), при котором максимизировался бы размер (дисперсия) проекции нашей выборки на него.
Замечание: Обобщение дисперсии на высшие размерности - ковариационная матрица, и эти два понятия эквивалентны. При проекции на вектор максимизируется дисперсия проекции, при проекции на пространства больших порядков – вся ее ковариационная матрица.
Итак, возьмем единичный вектор на который будем проецировать наш случайный вектор X. Тогда проекция на него будет равна v T X. Дисперсия проекции на вектор будет соответственно равна Var(v T X). В общем виде в векторной форме (для центрированных величин) дисперсия выражается так:
Соответственно, дисперсия проекции:
Легко заметить, что дисперсия максимизируется при максимальном значении v T Σv. Здесь нам поможет отношение Рэлея. Не вдаваясь слишком глубоко в математику, просто скажу, что у отношения Рэлея есть специальный случай для ковариационных матриц:
Последняя формула должна быть знакома по теме разложения матрицы на собственные вектора и значения. x является собственным вектором, а λ – собственным значением. Количество собственных векторов и значений равны размеру матрицы (и значения могут повторяться).
Кстати, в английском языке собственные значения и векторы именуются eigenvalues и eigenvectors соответственно.
Мне кажется, это звучит намного более красиво (и кратко), чем наши термины.
Таким образом, направление максимальной дисперсии у проекции всегда совпадает с айгенвектором имеющим максимальное собственное значение, равное величине этой дисперсии .
И это справедливо также для проекций на большее количество измерений – дисперсия (ковариационная матрица) проекции на m-мерное пространство будет максимальна в направлении m айгенвекторов, имеющих максимальные собственные значения.
Размерность нашей выборки равна двум и количество айгенвекторов у нее, соответственно, 2. Найдем их.
В библиотеке numpy реализована функция numpy.linalg.eig(X) , где X – квадратная матрица. Она возвращает 2 массива – массив айгензначений и массив айгенвекторов (векторы-столбцы). И векторы нормированы - их длина равна 1. Как раз то, что надо. Эти 2 вектора задают новый базис для выборки, такой что его оси совпадают с полуосями аппроксимирующего эллипса нашей выборки.
На этом графике мы апроксимировали нашу выборку эллипсом с радиусами в 2 сигмы (т.е. он должен содержать в себе 95% всех наблюдений – что в принципе мы здесь и наблюдаем). Я инвертировал больший вектор (функция eig(X) направляла его в обратную сторону) – нам важно направление, а не ориентация вектора.
Шаг 4. Снижение размерности (проекция)
Наибольший вектор имеет направление, схожее с линией регрессии и спроецировав на него нашу выборку мы потеряем информацию, сравнимую с суммой остаточных членов регрессии (только расстояние теперь евклидово, а не дельта по Y). В нашем случае зависимость между признаками очень сильная, так что потеря информации будет минимальна. «Цена» проекции - дисперсия по меньшему айгенвектору - как видно из предыдущего графика, очень невелика.
Замечание: диагональные элементы ковариационной матрицы показывают дисперсии по изначальному базису, а ее собственные значения – по новому (по главным компонентам).
Часто требуется оценить объем потерянной (и сохраненной) информации. Удобнее всего представить в процентах. Мы берем дисперсии по каждой из осей и делим на общую сумму дисперсий по осям (т.е. сумму всех собственных чисел ковариационной матрицы).
Таким образом, наш больший вектор описывает 45.994 / 46.431 * 100% = 99.06%, а меньший, соответственно, примерно 0.94%. Отбросив меньший вектор и спроецировав данные на больший, мы потеряем меньше 1% информации! Отличный результат!
Замечание: На практике, в большинстве случаев, если суммарная потеря информации составляет не более 10-20%, то можно спокойно снижать размерность.
Для проведения проекции, как уже упоминалось ранее на шаге 3, надо провести операцию v T X (вектор должен быть длины 1). Или, если у нас не один вектор, а гиперплоскость, то вместо вектора v T берем матрицу базисных векторов V T . Полученный вектор (или матрица) будет являться массивом проекций наших наблюдений.
V = (-vecs, -vecs) Xnew = dot(v,Xcentered) print Xnew OUT: [ -9.56404107 -9.02021624 -5.52974822 -2.96481262 0.68933859 0.74406645 2.33433492 7.39307974 5.3212742 10.59672425]
dot(X,Y) - почленное произведение (так мы перемножаем векторы и матрицы в Python)
Нетрудно заметить, что значения проекций соответствуют картине на предыдущем графике.
Шаг 5. Восстановление данных
С проекцией удобно работать, строить на ее основе гипотезы и разрабатывать модели. Но не всегда полученные главные компоненты будут иметь явный, понятный постороннему человеку, смысл. Иногда полезно раскодировать, к примеру, обнаруженные выбросы, чтобы посмотреть, что за наблюдения за ними стоят.
Это очень просто. У нас есть вся необходимая информация, а именно координаты базисных векторов в исходном базисе (векторы, на которые мы проецировали) и вектор средних (для отмены центровки). Возьмем, к примеру, наибольшее значение: 10.596… и раскодируем его. Для этого умножим его справа на транспонированный вектор и прибавим вектор средних, или в общем виде для всей выбоки: X T v T +m
Xrestored = dot(Xnew,v) + m print "Restored: ", Xrestored print "Original: ", X[:,9] OUT: Restored: [ 10.13864361 19.84190935] Original: [ 10. 19.9094105]
Разница небольшая, но она есть. Ведь потерянная информация не восстанавливается. Тем не менее, если простота важнее точности, восстановленное значение отлично аппроксимирует исходное.
Вместо заключения – проверка алгоритма
Итак, мы разобрали алгоритм, показали как он работает на игрушечном примере, теперь осталось только сравнить его с PCA, реализованным в sklearn – ведь пользоваться будем именно им.
From sklearn.decomposition import PCA pca = PCA(n_components = 1) XPCAreduced = pca.fit_transform(transpose(X))
Параметр n_components указывает на количество измерений, на которые будет производиться проекция, то есть до скольки измерений мы хотим снизить наш датасет. Другими словами – это n айгенвекторов с самыми большими собственными числами. Проверим результат снижения размерности:
Print "Our reduced X: n", Xnew print "Sklearn reduced X: n", XPCAreduced OUT: Our reduced X: [ -9.56404106 -9.02021625 -5.52974822 -2.96481262 0.68933859 0.74406645 2.33433492 7.39307974 5.3212742 10.59672425] Sklearn reduced X: [[ -9.56404106] [ -9.02021625] [ -5.52974822] [ -2.96481262] [ 0.68933859] [ 0.74406645] [ 2.33433492] [ 7.39307974] [ 5.3212742 ] [ 10.59672425]]
Мы возвращали результат как матрицу вектор-столбцов наблюдений (это более канонический вид с точки зрения линейной алгебры), PCA в sklearn же возвращает вертикальный массив.
В принципе, это не критично, просто стоит отметить, что в линейной алгебре канонично записывать матрицы через вектор-столбцы, а в анализе данных (и прочих связанных с БД областях) наблюдения (транзакции, записи) обычно записываются строками.
Проверим и прочие параметры модели – функция имеет ряд атрибутов, позволяющих получить доступ к промежуточным переменным:
Вектор средних: mean_
- Вектор(матрица) проекции: components_
- Дисперсии осей проекции (выборочная): explained_variance_
- Доля информации (доля от общей дисперсии): explained_variance_ratio_
Замечание: explained_variance_ показывает выборочную дисперсию, тогда как функция cov() для построения ковариационной матрицы рассчитывает несмещенные дисперсии!
Сравним полученные нами значения со значениями библиотечной функции.
Print "Mean vector: ", pca.mean_, m print "Projection: ", pca.components_, v print "Explained variance ratio: ", pca.explained_variance_ratio_, l/sum(l) OUT: Mean vector: [ 5.5 10.31439392] (5.5, 10.314393916) Projection: [[ 0.43774316 0.89910006]] (0.43774316434772387, 0.89910006232167594) Explained variance: [ 41.39455058] 45.9939450918 Explained variance ratio: [ 0.99058588] 0.990585881238
Единственное различие – в дисперсиях, но как уже упоминалось, мы использовали функцию cov(), которая использует несмещенную дисперсию, тогда как атрибут explained_variance_ возвращает выборочную. Они отличаются только тем, что первая для получения мат.ожидания делит на (n-1), а вторая – на n. Легко проверить, что 45.99 ∙ (10 - 1) / 10 = 41.39.
Все остальные значения совпадают, что означает, что наши алгоритмы эквивалентны. И напоследок замечу, что атрибуты библиотечного алгоритма имеют меньшую точность, поскольку он наверняка оптимизирован под быстродействие, либо просто для удобства округляет значения (либо у меня какие-то глюки).
Замечание: библиотечный метод автоматически проецирует на оси, максимизирующие дисперсию. Это не всегда рационально. К примеру, на данном рисунке неаккуратное снижение размерности приведет к тому, что классификация станет невозможна. Тем не менее, проекция на меньший вектор успешно снизит размерность и сохранит классификатор.
Итак, мы рассмотрели принципы работы алгоритма PCA и его реализации в sklearn. Я надеюсь, эта статья была достаточно понятна тем, кто только начинает знакомство с анализом данных, а также хоть немного информативна для тех, кто хорошо знает данный алгоритм. Интуитивное представление крайне полезно для понимания того как работает метод, а понимание очень важно для правильной настройки выбранной модели. Спасибо за внимание!
P.S.: Просьба не ругать автора за возможные неточности. Автор сам в процессе знакомства с дата-анализом и хочет помочь таким же как он в процессе освоения этой удивительной области знаний! Но конструктивная критика и разнообразный опыт всячески приветствуются!
В стремлении предельно точно описать исследуемую область аналитики часто отбирают большое число независимых переменных (p). В этом случае может возникнуть серьезная ошибка: несколько описывающих переменных могут характеризовать одну и ту же сторону зависимой переменной и, как следствие, высоко коррелировать между собой. Мультиколлинеарность независимых переменных серьезно искажает результаты исследования, поэтому от нее следует избавляться.
Метод главных компонент (как упрощенная модель факторного анализа, поскольку при этом методе не используются индивидуальные факторы, описывающие только одну переменную x i) позволяет объединить влияние высоко коррелированных переменных в один фактор, характеризующий зависимую переменную с одной единственной стороны. В результате анализа, осуществленного по методу главных компонент, мы добьемся сжатия информации до необходимых размеров, описания зависимой переменной m (m Для начала необходимо решить, сколько факторов необходимо выделить в данном исследовании. В рамках метода главных компонент первый главный фактор описывает наибольших процент дисперсии независимых переменных, далее – по убывающей. Таким образом, каждая следующая главная компонента, выделенная последовательно, объясняет все меньшую долю изменчивости факторов x i . Задача исследователя состоит в том, чтобы определить, когда изменчивость становится действительно малой и случайной. Другими словами – сколько главных компонент необходимо выбрать для дальнейшего анализа. Существует несколько методов рационального выделения необходимого числа факторов. Наиболее используемый из них – критерий Кайзера. Согласно этому критерию, отбираются только те факторы, собственные значения которых больше 1. Таким образом, фактор, который не объясняет дисперсию, эквивалентную, по крайней мере, дисперсии одной переменной, опускается. Проанализируем Таблицу 19, построенную в SPSS: Таблица 19. Полная объясненная дисперсия Как видно из Таблицы 19, в данном исследовании переменные x i высоко коррелирут между собой (это также выявлено ранее и видно из Таблицы 5 «Парные коэффициенты корреляции»), а следовательно, характеризуют зависимую переменную Y практически с одной стороны: изначально первая главная компонента объясняет 90,7 % дисперсии x i , и только собственное значение, соответствующее первой главной компоненте, больше 1. Конечно, это является недостатком отбора данных, однако в процессе самого отбора этот недостаток не был очевиден. Анализ в пакете SPSS позволяет самостоятельно выбрать число главных компонент. Выберем число 6 – равное количеству независимых переменных. Второй столбец Таблицы 19 показывает суммы квадратов нагрузок вращения, именно по этим результатам и сделаем вывод о числе факторов. Собственные значения, соответствующие первым двум главным компонентам, больше 1 (55,246% и 38,396% соответственно), поэтому, согласно методу Кайзера, выделим 2 наиболее значимые главные компоненты. Второй метод выделения необходимого числа факторов – критерий «каменистой осыпи». Согласно этому методу, собственные значения представляются в виде простого графика, и выбирается такое место на графике, где убывание собственных значений слева направо максимально замедляется: Рисунок 3. Критерий "каменистой осыпи" Как видно на Рисунке 3, убывание собственных значений замедляется уже со второй компоненты, однако постоянная скорость убывания (очень маленькая) начинается лишь с третьей компоненты. Следовательно, для дальнейшего анализа будут отобраны первые две главные компоненты. Это умозаключение согласуется с выводом, полученным при использовании метода Кайзера. Таким образом, окончательно выбираются первые две последовательно полученные главные компоненты. После выделения главных компонент, которые будут использоваться в дальнейшем анализе, необходимо определить корреляцию исходных переменных x i c полученными факторами и, исходя из этого, дать названия компонентам. Для анализа воспользуемся матрицей факторных нагрузок А, элементы которой являются коэффициентами корреляции факторов с исходными независимыми переменными: Таблица 20. Матрица факторных нагрузок В данном случае интерпретация коэффициентов корреляции затруднена, следовательно, довольно сложно дать названия первым двум главным компонентам. Поэтому далее воспользуемся методом ортогонального поворота системы координат Варимакс, целью которого является поворот факторов так, чтобы выбрать простейшую для интерпретации факторную структуру: Таблица 21. Коэффициенты интерпретации Из Таблицы 21 видно, что первая главная компонента больше всего связана с переменными x1, x2, x3; а вторая – с переменными x4, x5, x6. Таким образом, можно сделать вывод, что объем инвестиций в основные средства в регионе (переменная Y)
зависит от двух факторов: - объема собственных и заемных средств, поступивших в предприятия региона за период (первая компонента, z1);
-
а также от интенсивности вложений предприятий региона в финансовые активы и количества иностранного капитала в регионе (вторая компонента, z2).
Рисунок 4. Диаграмма рассеивания Данная диаграмма демонстрирует неутешительные результаты. Еще в самом начале исследования мы старались подобрать данные так, чтобы результирующая переменная Y была распределена нормально, и нам практически это удалось. Законы распределения независимых переменных были достаточно далеки от нормального, однако мы старались максимально приблизить их к нормальному закону (соответствующим образом подобрать данные). Рисунок 4 показывает, что первоначальная гипотеза о близости закона распределения независимых переменных к нормальному закону не подтверждается: форма облака должна напоминать эллипс, в центре объекты должны быть расположены более густо, нежели чем по краям. Стоит заметить, что сделать многомерную выборку, в которой все переменные распределены по нормальному закону – задача, выполнимая с огромным трудом (более того, не всегда имеющая решение). Однако к этой цели нужно стремиться: тогда результаты анализа будут более значимыми и понятными при интерпретации. К сожалению, в нашем случае, когда проделана большая часть работы по анализу собранных данных, менять выборку достаточно затруднительно. Но далее, в последующих работах, стоит более серьезно подходить в выборке независимых переменных и максимально приближать закон их распределения к нормальному. Последним этапом анализа методом главных компонент является построение уравнения регрессии на главные компоненты (в данном случае – на первую и вторую главные компоненты). При помощи SPSS рассчитаем параметры регрессионной модели: Таблица 22. Параметры уравнения регресии на главные компоненты Уравнение регрессии примет вид: y=47 414,184 + 0,916*z1+0,213*z2,
(b0) (b1) (b2)
т. о. b0
=47 414,184
показывает точку пересечения прямой регрессии с осью результирующего показателя; b1= 0,916 –
при увеличении значения фактора z1 на 1 ожидаемое среднее значение суммы объема инвестиций в основные средства увеличится на 0,916; b2= 0,213 –
при увеличении значения фактора z2 на 1 ожидаемое среднее значение суммы объема инвестиций в основные средства увеличится на 0,213. В данном случае значение tкр («альфа»=0,001, «ню»=53) = 3,46 меньше tнабл для всех коэффициентов «бета». Следовательно, все коэффициенты значимы. Таблица 24. Качество регрессионной модели на главные компоненты В Таблице 24 отражены показатели, которые характеризуют качество построенной модели, а именно: R – множественный к-т корреляции – говорит о том, какая доля дисперсии Y объясняется вариацией Z; R^2 – к-т детерминации – показывает долю объяснённой дисперсии отклонений Y от её среднего значения. Стандартная ошибка оценки характеризует ошибку построенной модели. Сравним эти показатели с аналогичными показателями степенной регрессионной модели (ее качество оказалось выше качества линейной модели, поэтому сравниваем именно со степенной): Таблица 25. Качество степенной регрессионной модели Так, множественный к-т корреляции R и к-т детерминации R^2 в степенной модели несколько выше, чем в модели главных компонент. Кроме того, стандартная ошибка модели главных компонент НАМНОГО выше, чем в степенной модели. Поэтому качество степенной регрессионной модели выше, чем регрессионной модели, построенной на главных компонентах. Проведем верификацию регрессионной модели главных компонент, т. е. проанализируем ее значимость. Проверим гипотезу о незначимости модели, рассчитаем F(набл.) = 204,784 (рассчитано в SPSS), F(крит) (0,001; 2; 53)=7,76. F(набл)>F(крит), следовательно, гипотеза о незначимости модели отвергается. Модель значима. Итак, в результате проведения компонентного анализа, было выяснено, что из отобранных независимых переменных x i можно выделить 2 главные компоненты – z1 и z2, причем на z1 в большей степени влияют переменные x1, x2, x3, а на z2 – x4, x5, x6. Уравнение регрессии, построенное на главных компонентах, оказалось значимым, хотя и уступает по качеству степенному уравнению регрессии. Согласно уравнению регрессии на главные компоненты, Y положительно зависит как от Z1, так и от Z2. Однако изначальная мультиколлинеарность переменных xi и то, что они не распределены по нормальному закону распределения, может искажать результаты построенной модели и делать ее менее значимой. Кластерный анализ
Следующим этапом данного исследования является кластерный анализ. Задачей кластерного анализа является разбиение выбранных регионов (n=56) на сравнительно небольшое число групп (кластеров) на основе их естественной близости относительно значений переменных x i . При проведении кластерного анализа мы предполагаем, что геометрическая близость двух или нескольких точек в пространстве означает физическую близость соответствующих объектов, их однородность (в нашем случае - однородность регионов по показателям, влияющим на инвестиции в основные средства). На первой стадии кластерного анализа необходимо определиться с оптимальным числом выделяемых кластеров. Для этого необходимо провести иерархическую кластеризацию – последовательное объединение объектов в кластеры до тех пор, пока не останется два больших кластера, объединяющиеся в один на максимальном расстоянии друг от друга. Результат иерархического анализа (вывод об оптимальном количестве кластеров) зависит от способа расчета расстояния между кластерами. Таким образом, протестируем различные методы и сделаем соответствующие выводы. Метод «ближнего соседа»
Если расстояние между отдельными объектами мы рассчитываем единым способом – как простое евклидово расстояние – расстояние между кластерами вычисляется разными методами. Согласно методу «ближайшего соседа», расстояние между кластерами соответствует минимальному расстоянию между двумя объектами разных кластеров. Анализ в пакете SPSS проходит следующим образом. Сначала рассчитывается матрица расстояний между всеми объектами, а затем, на основе матрицы расстояний, объекты последовательно объединяются в кластеры (для каждого шага матрица составляется заново). Шаги последовательного объединения представлены в таблице: Таблица 26. Шаги агломерации. Метод «ближайшего соседа» Как видно из Таблицы 26, на первом этапе объединились элементы 7 и 8, т. к. расстояние между ними было минимальным – 0,003. Далее расстояние между объединенными объектами увеличивается. По таблице также можно сделать вывод об оптимальном числе кластеров. Для этого нужно посмотреть, после какого шага происходит резкий скачок в величине расстояния, и вычесть номер этой агломерации из числа исследуемых объектов. В нашем случае: (56-53)=3 – оптимальное число кластеров. Рисунок 5. Дендрограмма. Метод "ближайшего соседа" Аналогичный вывод об оптимальном количестве кластеров можно сделать и глядя на дендрограмму (Рис. 5): следует выделить 3 кластера, причем в первый кластер войдут объекты под номерами 1-54 (всего 54 объекта), а во второй и третий кластеры – по одному объекту (под номерами 55 и 56 соответственно). Данный результат говорит о том, что первые 54 региона относительно однородны по показателям, влияющим на инвестиции в основные средства, в то время как объекты под номерами 55 (Республика Дагестан) и 56 (Новосибирская область) значительно выделяются на общем фоне. Стоит заметить, что данные субъекты имеют самые большие объемы инвестиций в основные средства среди всех отобранных регионов. Этот факт еще раз доказывает высокую зависимость результирующей переменной (объема инвестиций) от выбранных независимых переменных. Аналогичные рассуждения проводятся для других методов расчета расстояния между кластерами. Метод «дальнего соседа»
Таблица 27. Шаги агломерации. Метод "дальнего соседа" При методе «дальнего соседа» расстояние между кластерами рассчитывается как максимальное расстояние между двумя объектами в двух разных кластерах. Согласно Таблице 27, оптимальное число кластеров равно (56-53)=3. Рисунок 6. Дендрограмма. Метод "дальнего соседа" Согласно дендрограмме, оптимальным решением также будет выделение 3 кластеров: в первый кластер войдут регионы под номерами 1-50 (50 регионов), во второй – под номерами 51-55 (5 регионов), в третий – последний регион под номером 56. Метод «центра тяжести»
При методе «центра тяжести» за расстояние между кластерами принимается евклидово расстояние между «центрами тяжести» кластеров – средними арифметическими их показателей x i . Рисунок 7. Дендрограмма. Метод "центра тяжести" На Рисунке 7 видно, что оптимальное число кластеров следующее: 1 кластер – 1-47 объекты; 2 кластер – 48-54 объекты (всего 6); 3 кластер – 55 объект; 4 кластер – 56 объект. Принцип «средней связи»
В данном случае расстояние между кластерами равно среднему значению расстояний между всеми возможными парами наблюдений, причем одно наблюдение берется из одного кластера, а второе – соответственно, из другого. Анализ таблицы шагов агломерации показал, что оптимальное количество кластеров равно (56-52)=4. Сравним этот вывод с выводом, полученным при анализе дендрограммы. На Рисунке 8 видно, что в 1 кластер войдут объекты под номерами 1-50, во 2 кластер – объекты 51-54 (4 объекта), в 3 кластер – 55 регион, в 4 кластер – 56 регион. Рисунок 8. Дендрограмма. Метод "средней связи" Компонентный анализ относится к многомерным методам снижения размерности. Он содержит один метод - метод главных компонент. Главные компоненты представляют собой ортогональную систему координат, в которой дисперсии компонент характеризуют их статистические свойства. Учитывая, что объекты исследования в экономике характеризуются большим, но конечным количеством признаков, влияние которых подвергается воздействию большого количества случайных причин. Первой главной компонентой Z1 исследуемой системы признаков Х1, Х2, Х3 , Х4 ,…, Хn называется такая центрировано - нормированная линейная комбинация этих признаков, которая среди прочих центрировано - нормированных линейных комбинаций этих признаков, имеет дисперсию наиболее изменчивую. В качестве второй главной компоненты Z2 мы будем брать такую центрировано - нормированную комбинацию этих признаков, которая: не коррелированна с первой главной компонентой, не коррелированны с первой главной компонентой, эта комбинация имеет наибольшую дисперсию. K-ой главной компонентой Zk (k=1…m) мы будем называть такую центрировано - нормированную комбинацию признаков, которая: не коррелированна с к-1 предыдущими главными компонентами, среди всех возможных комбинаций исходных признаков, которые не не коррелированны с к-1 предыдущими главными компонентами, эта комбинация имеет наибольшую дисперсию. Введём ортогональную матрицу U и перейдём от переменных Х к переменным Z, причём Вектор выбирается т. о., чтобы дисперсия была максимальной. После получения выбирается т. о., чтобы дисперсия была максимальной при условии, что не коррелированно с и т. д. Так как признаки измерены в несопоставимых величинах, то удобнее будет перейти к центрированно-нормированным величинам. Матрицу исходных центрированно-нормированных значений признаков найдем из соотношения: где - несмещенная, состоятельная и эффективная оценка математического ожидания, Несмещенная, состоятельная и эффективная оценка дисперсии. Матрица наблюденных значений исходных признаков приведена в Приложении. Центрирование и нормирование произведено с помощью программы"Stadia". Так как признаки центрированы и нормированы, то оценку корреляционной матрицы можно произвести по формуле: Перед тем как проводить компонентный анализ, проведем анализ независимости исходных признаков. Проверка значимости матрицы парных корреляций с помощью критерия Уилкса. Выдвигаем гипотезу: Н0: незначима Н1: значима 125,7; (0,05;3,3) = 7,8 т.к > , то гипотеза Н0 отвергается и матрица является значимой, следовательно, имеет смысл проводить компонентный анализ. Проверим гипотезу о диагональности ковариационной матрицы Выдвигаем гипотезу: Строим статистику, распределена по закону с степенями свободы. 123,21, (0,05;10) =18,307 т.к >, то гипотеза Н0 отвергается и имеет смысл проводить компонентный анализ. Для построения матрицы факторных нагрузок необходимо найти собственные числа матрицы, решив уравнение. Используем для этой операции функцию eigenvals системы MathCAD, которая возвращает собственные числа матрицы: Т.к. исходные данные представляют собой выборку из генеральной совокупности, то мы получили не собственные числа и собственные вектора матрицы, а их оценки. Нас будет интересовать на сколько “хорошо” со статистической точки зрения выборочные характеристики описывают соответствующие параметры для генеральной совокупности. Доверительный интервал для i-го собственного числа ищется по формуле: Доверительные интервалы для собственных чисел в итоге принимают вид: Оценка значения нескольких собственных чисел попадает в доверительный интервал других собственных чисел. Необходимо проверить гипотезу о кратности собственных чисел. Проверка кратности производится с помощью статистики где r-количество кратных корней. Данная статистика в случае справедливости распределена по закону с числом степеней свободы. Выдвинем гипотезы: Так как, то гипотеза отвергается, то есть собственные числа и не кратны. Так как, то гипотеза отвергается, то есть собственные числа и не кратны. Необходимо выделить главные компоненты на уровне информативности 0,85. Мера информативности показывает какую часть или какую долю дисперсии исходных признаков составляют k-первых главных компонент. Мерой информативности будем называть величину: На заданном уровне информативности выделено три главных компоненты. Запишем матрицу = Для получения нормализованного вектора перехода от исходных признаков к главным компонентам необходимо решить систему уравнений: , где - соответствующее собственное число. После получения решения системы необходимо затем нормировать полученный вектор. Для решения данной задачи воспользуемся функцией eigenvec системы MathCAD, которая возвращает нормированный вектор для соответствующего собственного числа. В нашем случае первых четырех главных компонент достаточно для достижения заданного уровня информативности, поэтому матрица U (матрица перехода от исходного базиса к базису из собственных векторов) Строим матрицу U, столбцами которой являются собственные вектора: Матрица весовых коэффициентов: Коэффициенты матрицы А являются коэффициентами корреляции между центрировано - нормированными исходными признаками и ненормированными главными компонентами, и показывают наличие, силу и направление линейной связи между соответствующими исходными признаками и соответствующими главными компонентами. Метод главных компонент
Метод главных компонент
(англ. Principal component analysis, PCA
) - один из основных способов уменьшить размерность данных, потеряв наименьшее количество информации . Изобретен К. Пирсоном (англ. Karl Pearson
) в г. Применяется во многих областях, таких как распознавание образов , компьютерное зрение , сжатие данных и т. п. Вычисление главных компонент сводится к вычислению собственных векторов и собственных значений ковариационной матрицы исходных данных. Иногда метод главных компонент называют преобразованием Кархунена-Лоэва
(англ. Karhunen-Loeve
) или преобразованием Хотеллинга (англ. Hotelling transform
). Другие способы уменьшения размерности данных - это метод независимых компонент, многомерное шкалирование, а также многочисленные нелинейные обобщения: метод главных кривых и многообразий, метод упругих карт , поиск наилучшей проекции (англ. Projection Pursuit
), нейросетевые методы «узкого горла », и др. Задача анализа главных компонент, имеет, как минимум, четыре базовых версии: Первые три версии оперируют конечными множествами данных. Они эквивалентны и не используют никакой гипотезы о статистическом порождении данных. Четвёртая версия оперирует случайными величинами. Конечные множества появляются здесь как выборки из данного распределения, а решение трёх первых задач - как приближение к «истинному» преобразованию Кархунена-Лоэва. При этом возникает дополнительный и не вполне тривиальный вопрос о точности этого приближения. Иллюстрация к знаменитой работе К. Пирсона (1901): даны точки на плоскости, - расстояние от до прямой . Ищется прямая , минимизирующая сумму Метод главных компонент начинался с задачи наилучшей аппроксимации конечного множества точек прямыми и плоскостями (К. Пирсон, 1901). Дано конечное множество векторов . Для каждого среди всех -мерных линейных многообразий в найти такое , что сумма квадратов уклонений от минимальна: где - евклидово расстояние от точки до линейного многообразия. Всякое -мерное линейное многообразие в может быть задано как множество линейных комбинаций , где параметры пробегают вещественную прямую , а - ортонормированный набор векторов где евклидова норма, - евклидово скалярное произведение, или в координатной форме: Решение задачи аппроксимации для даётся набором вложенных линейных многообразий , . Эти линейные многообразия определяются ортонормированным набором векторов (векторами главных компонент) и вектором . Вектор ищется, как решение задачи минимизации для : Векторы главных компонент могут быть найдены как решения однотипных задач оптимизации : 1) централизуем данные (вычитаем среднее): . Теперь ;
2) находим первую главную компоненту как решение задачи;
.
Если решение не единственно, то выбираем одно из них.
3) Вычитаем из данных проекцию на первую главную компоненту:
;
4) находим вторую главную компоненту как решение задачи
.
Если решение не единственно, то выбираем одно из них.
…
2k-1) Вычитаем проекцию на -ю главную компоненту (напомним, что проекции на предшествующие главные компоненты уже вычтены):
;
2k) находим k-ю главную компоненту как решение задачи:
.
Если решение не единственно, то выбираем одно из них.
…
На каждом подготовительном шаге вычитаем проекцию на предшествующую главную компоненту. Найденные векторы ортонормированы просто в результате решения описанной задачи оптимизации, однако чтобы не дать ошибкам вычисления нарушить взаимную ортогональность векторов главных компонент, можно включать в условия задачи оптимизации. Неединственность в определении помимо тривиального произвола в выборе знака ( и решают ту же задачу) может быть более существенной и происходить, например, из условий симметрии данных. Последняя главная компонента - единичный вектор, ортогональный всем предыдущим . Первая главная компонента максимизирует выборочную дисперсию проекции данных Пусть нам дан центрированный набор векторов данных (среднее арифметическое значение равно нулю). Задача - найти такое ортогональное преобразование в новую систему координат , для которого были бы верны следующие условия: Теория сингулярного разложения была создана Дж. Дж. Сильвестром (англ. James Joseph Sylvester
) в г. и изложена во всех подробных руководствах по теории матриц . Основная процедура - поиск наилучшего приближения произвольной матрицы матрицей вида (где - -мерный вектор, а - -мерный вектор) методом наименьших квадратов: Решение этой задачи дается последовательными итерациями по явным формулам. При фиксированном векторе значения , доставляющие минимум форме , однозначно и явно определяются из равенств : Аналогично, при фиксированном векторе определяются значения : B качестве начального приближения вектора возьмем случайный вектор единичной длины, вычисляем вектор , далее для этого вектора вычисляем вектор и т. д. Каждый шаг уменьшает значение . В качестве критерия остановки используется малость относительного уменьшения значения минимизируемого функционала за шаг итерации () или малость самого значения . В результате для матрицы получили наилучшее приближение матрицей вида (здесь верхним индексом обозначен номер приближения). Далее, из матрицы вычитаем полученную матрицу , и для полученной матрицы уклонений вновь ищем наилучшее приближение этого же вида и т. д., пока, например, норма не станет достаточно малой. В результате получили итерационную процедуру разложения матрицы в виде суммы матриц ранга 1, то есть . Полагаем и нормируем векторы : В результате получена аппроксимация сингулярных чисел и сингулярных векторов (правых - и левых - ). К достоинствам этого алгоритма относится его исключительная простота и возможность почти без изменений перенести его на данные с пробелами , а также взвешенные данные. Существуют различные модификации базового алгоритма, улучшающие точность и устойчивость. Например, векторы главных компонент при разных должны быть ортогональны «по построению», однако при большом числе итерации (большая размерность, много компонент) малые отклонения от ортогональности накапливаются и может потребоваться специальная коррекция на каждом шаге, обеспечивающая его ортогональность ранее найденным главным компонентам. Часто вектор данных имеет дополнительную структуру прямоугольной таблицы (например, плоское изображение) или даже многомерной таблицы - то есть тензора : , . В этом случае также эффективно применять сингулярное разложение. Определение, основные формулы и алгоритмы переносятся практически без изменений: вместо матрицы данных имеем -индексную величину , где первый индекс -номер точки (тензора) данных. Основная процедура - поиск наилучшего приближения тензора тензором вида (где - -мерный вектор ( - число точек данных), - вектор размерности при ) методом наименьших квадратов: Решение этой задачи дается последовательными итерациями по явным формулам. Если заданы все векторы-сомножители кроме одного , то этот оставшийся определяется явно из достаточных условий минимума. B качестве начального приближения векторов () возьмем случайные векторы единичной длины, вычислим вектор , далее для этого вектора и данных векторов вычисляем вектор и т. д. (циклически перебирая индексы) Каждый шаг уменьшает значение . Алгоритм, очевидно, сходится. В качестве критерия остановки используется малость относительного уменьшения значения минимизируемого функционала за цикл или малость самого значения . Далее, из тензора вычитаем полученное приближение и для остатка вновь ищем наилучшее приближение этого же вида и т. д., пока, например, норма очередного остатка не станет достаточно малой. Это многокомпонентное сингулярное разложение (тензорный метод главных компонент) успешно применяется при обработке изображений, видеосигналов, и, шире, любых данных, имеющих табличную или тензорную структуру. Матрица преобразования данных к главным компонентам состоит из векторов главных компонент, расположенных в порядке убывания собственных значений: То есть, матрица является ортогональной . Большая часть вариации данных будет сосредоточена в первых координатах, что позволяет перейти к пространству меньшей размерности. Пусть данные центрированы, . При замене векторов данных на их проекцию на первые главных компонент вносится средний квадрат ошибки в расчете на один вектор данных: где собственные значения эмпирической ковариационной матрицы , расположенные в порядке убывания, с учетом кратности. Эта величина называется остаточной дисперсией
. Величина называется объяснённой дисперсией
. Их сумма равна выборочной дисперсии. Соответствующий квадрат относительной ошибки - это отношение остаточной дисперсии к выборочной дисперсии (то есть доля необъяснённой дисперсии
): По относительной ошибке оценивается применимость метода главных компонент с проецированием на первые компонент. Замечание
: в большинстве вычислительных алгоритмов собственные числа с соответствующими собственными векторами - главными компонентами вычисляются в порядке «от больших - к меньшим». Для вычисления достаточно вычислить первые собственных чисел и след эмпирической ковариационной матрицы , (сумму диагональных элементов , то есть дисперсий по осям). Тогда Целевой подход к оценке числа главных компонент по необходимой доле объяснённой дисперсии формально применим всегда, однако неявно он предполагает, что нет разделения на «сигнал» и «шум», и любая заранее заданная точность имеет смысл. Поэтому часто более продуктивна иная эвристика, основывающаяся на гипотезе о наличии «сигнала» (сравнительно малая размерность, относительно большая амплитуда) и «шума» (большая размерность, относительно малая амплитуда). С этой точки зрения метод главных компонент работает как фильтр: сигнал содержится, в основном, в проекции на первые главные компоненты, а в остальных компонентах пропорция шума намного выше. Вопрос: как оценить число необходимых главных компонент, если отношение «сигнал/шум» заранее неизвестно? Простейший и старейший метод отбора главных компонент даёт правило Кайзера
(англ. Kaiser"s rule
): значимы те главные компоненты, для которых то есть превосходит среднее значение (среднюю выборочную дисперсию координат вектора данных). Правило Кайзера хорошо работает в простых случаях, когда есть несколько главных компонент с , намного превосходящими среднее значение, а остальные собственные числа меньше него. В более сложных случаях оно может давать слишком много значимых главных компонент. Если данные нормированы на единичную выборочную дисперсию по осям, то правило Кайзера приобретает особо простой вид: значимы только те главные компоненты, для которых Пример: оценка числа главных компонент по правилу сломанной трости в размерности 5. Одним из наиболее популярных эвристических подходов к оценке числа необходимых главных компонент является правило сломанной трости
(англ. Broken stick model
) . Набор нормированных на единичную сумму собственных чисел (, ) сравнивается с распределением длин обломков трости единичной длины, сломанной в -й случайно выбранной точке (точки разлома выбираются независимо и равнораспределены по длине трости). Пусть () - длины полученных кусков трости, занумерованные в порядке убывания длины: . Нетрудно найти математическое ожидание : По правилу сломанной трости -й собственный вектор (в порядке убывания собственных чисел ) сохраняется в списке главных компонент, если На Рис. приведён пример для 5-мерного случая: Для примера выбрано По правилу сломанной трости в этом примере следует оставлять 2 главных компоненты: По оценкам пользователей, правило сломанной трости имеет тенденцию занижать количество значимых главных компонент. После
проецирования на первые главных компонент с удобно произвести нормировку на единичную (выборочную) дисперсию по осям. Дисперсия вдоль й главной компоненты равна ), поэтому для нормировки надо разделить соответствующую координату на . Это преобразование не является ортогональным и не сохраняет скалярного произведения. Ковариационная матрица проекции данных после нормировки становится единичной, проекции на любые два ортогональных направления становятся независимыми величинами, а любой ортонормированный базис становится базисом главных компонент (напомним, что нормировка меняет отношение ортогональности векторов). Отображение из пространства исходных данных на первые главных компонент вместе с нормировкой задается матрицей Именно это преобразование чаще всего называется преобразованием Кархунена-Лоэва. Здесь - векторы-столбцы, а верхний индекс означает транспонирование. Предупреждение
: не следует путать нормировку, проводимую после преобразования к главным компонентам, с нормировкой и «обезразмериванием» при предобработке данных
, проводимой до вычисления главных компонент. Предварительная нормировка нужна для обоснованного выбора метрики, в которой будет вычисляться наилучшая аппроксимация данных, или будут искаться направления наибольшего разброса (что эквивалентно). Например, если данные представляют собой трёхмерные векторы из «метров, литров и килограмм», то при использовании стандартного евклидового расстояния разница в 1 метр по первой координате будет вносить тот же вклад, что разница в 1 литр по второй, или в 1 кг по третьей. Обычно системы единиц, в которых представлены исходные данные, недостаточно точно отображают наши представления о естественных масштабах по осям, и проводится «обезразмеривание»: каждая координата делится на некоторый масштаб, определяемый данными, целями их обработки и процессами измерения и сбора данных. Есть три cущественно различных стандартных подхода к такой нормировке: на единичную дисперсию
по осям (масштабы по осям равны средним квадратичным уклонениям - после этого преобразования ковариационная матрица совпадает с матрицей коэффициентов корреляции), на равную точность измерения
(масштаб по оси пропорционален точности измерения данной величины) и на равные требования
в задаче (масштаб по оси определяется требуемой точностью прогноза данной величины или допустимым её искажением - уровнем толерантности). На выбор предобработки влияют содержательная постановка задачи, а также условия сбора данных (например, если коллекция данных принципиально не завершена и данные будут ещё поступать, то нерационально выбирать нормировку строго на единичную дисперсию, даже если это соответствует смыслу задачи, поскольку это предполагает перенормировку всех данных после получения новой порции; разумнее выбрать некоторый масштаб, грубо оценивающий стандартное отклонение, и далее его не менять). Предварительная нормировка на единичную дисперсию по осям разрушается поворотом системы координат, если оси не являются главными компонентами, и нормировка при предобработке данных не заменяет нормировку после приведения к главным компонентам. Если сопоставить каждому вектору данных единичную массу, то эмпирическая ковариационная матрица совпадёт с тензором инерции этой системы точечных масс (делённым на полную массу ), а задача о главных компонентых - с задачей приведения тензора инерции к главным осям. Можно использовать дополнительную свободу в выборе значений масс для учета важности точек данных или надежности их значений (важным данным или данным из более надежных источников приписываются бо́льшие массы). Если вектору данных придаётся масса ,
то вместо эмпирической ковариационной матрицы получим Все дальнейшие операции по приведению к главным компонентам производятся так же, как и в основной версии метода: ищем ортонормированный собственный базис , упорядочиваем его по убыванию собственных значений, оцениваем средневзвешенную ошибку аппроксимации данных первыми компонентами (по суммам собственных чисел ), нормируем и т. п. Более общий способ взвешивания даёт максимизация взвешенной суммы попарных расстояний
между проекциями. Для каждых двух точек данных, вводится вес ; и . Вместо эмпирической ковариационной матрицы используется При симметричная матрица положительно определена, поскольку положительна квадратичная форма: Далее ищем ортонормированный собственный базис , упорядочиваем его по убыванию собственных значений, оцениваем средневзвешенную ошибку аппроксимации данных первыми компонентами и т. д. - в точности так же, как и в основном алгоритме. Этот способ применяется при наличии классов
: для из разных классов вес вес выбирается бо́льшим, чем для точек одного класса. В результате, в проекции на взвешенные главные компоненты различные классы «раздвигаются» на большее расстояние. Другое применение - снижение влияния больших уклонений
(оутлайеров, англ. Outlier
), которые могут искажать картину из-за использования среднеквадратичного расстояния: если выбрать , то влияние больших уклонений будет уменьшено. Таким образом, описанная модификация метода главных компонент является более робастной , чем классическая. В статистике при использовании метода главных компонент используют несколько специальных терминов. Матрица данных
; каждая строка - вектор предобработанных
данных (центрированных
и правильно нормированных
), число строк - (количество векторов данных), число столбцов - (размерность пространства данных); Матрица нагрузок
(Loadings) ; каждый столбец - вектор главных компонент, число строк - (размерность пространства данных), число столбцов - (количество векторов главных компонент, выбранных для проецирования); Матрица счетов
(Scores) ; каждая строка - проекция вектора данных на главных компонент; число строк - (количество векторов данных), число столбцов - (количество векторов главных компонент, выбранных для проецирования); Матрица Z-счетов
(Z-scores) ; каждая строка - проекция вектора данных на главных компонент, нормированная на единичную выборочную дисперсию; число строк - (количество векторов данных), число столбцов - (количество векторов главных компонент, выбранных для проецирования); Матрица ошибок
(или остатков
) (Errors or residuals) . Основная формула: Метод главных компонент применим всегда. Распространённое утверждение о том, что он применим только к нормально распределённым данным (или для распределений, близких к нормальным) неверно: в исходной формулировке К. Пирсона ставится задача об аппроксимации
конечного множества данных и отсутствует даже гипотеза о их статистическом порождении, не говоря уж о распределении. Однако метод не всегда эффективно снижает размерность при заданных ограничениях на точность . Прямые и плоскости не всегда обеспечивают хорошую аппроксимацию. Например, данные могут с хорошей точностью следовать какой-нибудь кривой, а эта кривая может быть сложно расположена в пространстве данных. В этом случае метод главных компонент для приемлемой точности потребует нескольких компонент (вместо одной), или вообще не даст снижения размерности при приемлемой точности. Для работы с такими «кривыми» главными компонентами изобретен метод главных многообразий и различные версии нелинейного метода главных компонент . Больше неприятностей могут доставить данные сложной топологии. Для их аппроксимации также изобретены различные методы, например самоорганизующиеся карты Кохонена , нейронный газ или топологические грамматики . Если данные статистически порождены с распределением, сильно отличающимся от нормального, то для аппроксимации распределения полезно перейти от главных компонент к независимым компонентам
, которые уже не ортогональны в исходном скалярном произведении. Наконец, для изотропного распределения (даже нормального) вместо эллипсоида рассеяния получаем шар, и уменьшить размерность методами аппроксимации невозможно. Визуализация данных - представление в наглядной форме данных эксперимента или результатов теоретического исследования. Первым выбором в визуализации множества данных является ортогональное проецирование на плоскость первых двух главных компонент (или 3-мерное пространство первых трёх главных компонент). Плоскость проектирования является, по сути плоским двумерным «экраном», расположенным таким образом, чтобы обеспечить «картинку» данных с наименьшими искажениями. Такая проекция будет оптимальна (среди всех ортогональных проекций на разные двумерные экраны) в трех отношениях: Визуализация данных является одним из наиболее широко используемых приложений метода главных компонент и его нелинейных обобщений . Для уменьшения пространственной избыточности пикселей при кодировании изображений и видео используется линейные преобразования блоков пикселей. Последующие квантования полученных коэффициентов и кодирование без потерь позволяют получить значительные коэффициенты сжатия. Использование преобразования PCA в качестве линейного преобразования является для некоторых типов данных оптимальным с точки зрения размера полученных данных при одинаковом искажении . На данный момент этот метод активно не используется, в основном из-за большой вычислительной сложности. Также сжатия данных можно достичь отбрасывая последние коэффициенты преобразования. Метод главных компонент - один из основных методов в хемометрике (англ. Chemometrics
). Позволяет разделить матрицу исходных данных X на две части: «содержательную» и «шум». По наиболее популярному определению «Хемометрика - это химическая дисциплина, применяющая математические, статистические и другие методы, основанные на формальной логике, для построения или отбора оптимальных методов измерения и планов эксперимента, а также для извлечения наиболее важной информации при анализе экспериментальных данных». В политологии метод главных компонент был основным инструментом проекта «Политический Атлас Современности» для линейного и нелинейного анализа рейтингов 192 стран мира по пяти специально разработанным интегральным индексам (уровня жизни, международного влияния, угроз, государственности и демократии). Для картографии результатов этого анализа разработана специальная ГИС (Геоинформационная система), объединяющая географическое пространство с пространством признаков. Также созданы карты данных политического атласа , использующие в качестве подложки двумерные главные многообразия в пятимерном пространстве стран. Отличие карты данных от географической карты заключается в том, что на географической карте рядом оказываются объекты, которые имеют сходные географические координаты, в то время как на карте данных рядом оказываются объекты (страны) с похожими признаками (индексами).
Компонента
Начальные собственные значения
Суммы квадратов нагрузок вращения
Итого
% Дисперсии
Кумулятивный %
Итого
% Дисперсии
Кумулятивный %
dimension0
5,442
90,700
90,700
3,315
55,246
55,246
,457
7,616
98,316
2,304
38,396
93,641
,082
1,372
99,688
,360
6,005
99,646
,009
,153
99,841
,011
,176
99,823
,007
,115
99,956
,006
,107
99,930
,003
,044
100,000
,004
,070
100,000
Метод выделения: Анализ главных компонент.
Матрица компонент a
Компонента
X1
,956
-,273
,084
,037
-,049
,015
X2
,986
-,138
,035
-,080
,006
,013
X3
,963
-,260
,034
,031
,060
-,010
X4
,977
,203
,052
-,009
-,023
-,040
X5
,966
,016
-,258
,008
-,008
,002
X6
,861
,504
,060
,018
,016
,023
Метод выделения: Анализ методом главных компонент.
a. Извлеченных компонент: 6
Матрица повернутых компонент a
Компонента
X1
,911
,384
,137
-,021
,055
,015
X2
,841
,498
,190
,097
,000
,007
X3
,900
,390
,183
-,016
-,058
-,002
X4
,622
,761
,174
,022
,009
,060
X5
,678
,564
,472
,007
,001
,005
X6
,348
,927
,139
,001
-,004
-,016
Метод выделения: Анализ методом главных компонент.
Метод вращения: Варимакс с нормализацией Кайзера.
a. Вращение сошлось за 4 итераций.
Модель
Нестандартизованные коэффициенты
Стандартизованные коэффициенты
t
Знч.
B
Стд. Ошибка
Бета
(Константа)
47414,184
1354,505
35,005
,001
Z1
26940,937
1366,763
,916
19,711
,001
Z2
6267,159
1366,763
,213
4,585
,001
Модель
R
R-квадрат
Скорректированный R-квадрат
Стд. ошибка оценки
dimension0
,941 a
,885
,881
10136,18468
a. Предикторы: (конст) Z1, Z2
b. Зависимая переменная: Y
Этап
Кластер объединен с
Коэффициенты
Следующий этап
Кластер 1
Кластер 2
Кластер 1
Кластер 2
,003
,004
,004
,005
,005
,005
,005
,006
,007
,007
,009
,010
,010
,010
,010
,011
,012
,012
,012
,012
,012
,013
,014
,014
,014
,014
,015
,015
,016
,017
,018
,018
,019
,019
,020
,021
,021
,022
,024
,025
,027
,030
,033
,034
,042
,052
,074
,101
,103
,126
,163
,198
,208
,583
1,072
Этап
Кластер объединен с
Коэффициенты
Этап первого появления кластера
Следующий этап
Кластер 1
Кластер 2
Кластер 1
Кластер 2
,003
,004
,004
,005
,005
,005
,005
,007
,009
,010
,010
,011
,011
,012
,012
,014
,014
,014
,017
,017
,018
,018
,019
,021
,022
,026
,026
,027
,034
,035
,035
,037
,037
,042
,044
,046
,063
,077
,082
,101
,105
,117
,126
,134
,142
,187
,265
,269
,275
,439
,504
,794
,902
1,673
2,449
Вычисление главных компонент
Формальная постановка задачи
Аппроксимация данных линейными многообразиями
Поиск ортогональных проекций с наибольшим рассеянием
Простой итерационный алгоритм сингулярного разложения
Сингулярное разложение тензоров и тензорный метод главных компонент
Матрица преобразования к главным компонентам
Остаточная дисперсия
Отбор главных компонент по правилу Кайзера
Оценка числа главных компонент по правилу сломанной трости
Нормировка
Нормировка после приведения к главным компонентам
Нормировка до вычисления главных компонент
Механическая аналогия и метод главных компонент для взвешенных данных
Специальная терминология
Пределы применимости и ограничения эффективности метода
Примеры использования
Визуализация данных
Компрессия изображений и видео
Подавление шума на изображениях
Хемометрика
Психодиагностика