Анализ главных компонент: 1я ГК имеет маленькую дисперсию

[Big_Z]

Здравствуйте, уважаемые участники форума!
В Matlab осуществляю МГК. Следующим кодом стандартизирую матрицу данных "объект-признак", делаю разложение и строю проекции на первые две ГК. Матрица имеет изначально много нулей, отрицательных нет, остальные значения примерно равномерно распределены на отрезке (0,10000). В итоге сингулярные числа убывают очень медленно. Объяснённая дисперсия на каждую ГК такая: 3%, 2,5%, 1,3% и т.д. Как это интерпретировать? Получается, если я посчитаю корреляцию ГК с признаками и сгруппирую признаки как сильно коррелированные (положительно или отрицательно) с конкретными ГК, то эти результаты нельзя считать надёжными?

%X = matr(dr,dc);
for j=1:dc
  for i=1:dr
     mu(j) = mu(j)+X(i,j);
  end;
  mu(j) = mu(j)/dr;
end;

for j=1:dc
  for i=1:dr
     si(j) = si(j)+(X(i,j)-mu(j))^2;
  end;
  si(j) = sqrt(si(j)/(dr-1));
end;

for i=1:dr
  for j=1:dc 
     X(i,j) = (X(i,j) - mu(j))/si(j);
  end;
end; 

[U,W,V] = svd(X);
lam = 100*diag(W.^2)/sum(diag(W.^2));
for i = 1:dc
    fprintf('%d %f\n',i,lam(i));
end; 

 T = U*W;
 
figure
scatter(T(:,1), T(:,2), 'g');
labels = num2str((1:dr)','%d');    %'
text(T(:,1), T(:,2), labels, 'horizontal','left', 'vertical','bottom');
grid on 

[Andrey_Kireew]

На мой взгляд вы просто некорректно применяете метод главных компонент, отсюда все недоразумения.
Поиск главных компонент действительно основывается на сингулярном разложении, но не исходной матрицы данных, а ковариационной матрицы признаков. Т.е. Вам нужно сделать вот так:

[U,W,V] = svd(X'*X);

Тогда матрица W будет диагональной матрицей сингулярных чисел (для чего Вы их возводите в квадрат - не совсем понятно).
Чтобы найти сами главные компоненты, нужно сделать так:

T = X*U;

Обратите внимание, что при нахождении Г.К. сингулярные числа не используются. Эти числа всего лишь показывают какая доля от общей дисперсии данных приходится на каждую главную компоненту.
По поводу медленного убывания сингулярных чисел - ответ вполне определённый: так ведёт себя шум. В данных, содержащих полезную информацию присутствует несколько больших сингулярных чисел, значения которых резко убывают. Связанные с ними главные компоненты интерпретируются как латентные факторы, в которых сосредоточены вся полезная информация. Остальные сингулярные числа убывают медленно, связанные с ними главные компоненты интерпретируются как шум.
В вашем случае чтобы говорить о чём то конкретном, нужно сначала правильно найти собственные числа. Скорее всего их распределение будет выглядеть так как я описал.

Кстати,нормирование данных, как это сделано у Вас (деление на si(j) в 18 строке) далеко не всегда оправдано. Если Вы предполагаете, что дисперсия шума во всех признаках примерно одинакова, то этого лучше не делать.

Хотелось бы предостеречь Вас так же от попыток осуществить отбор главных компонент на основании их корреляций с результативной переменной. Наиболее адекватным критерием значимости главных компонент является их мощность, которая определяется величиной соответствующих собственных чисел. Дело в том что на слабых главных компонентах могут, а скорее всего и будут, наблюдаться сильные корреляции. Это поначалу внушает оптимизм, но в действительности это ложные корреляции, которые характеризуют особенности обучающей выборки и не воспроизводятся на неизвестных данных.

В некоторых литературных источниках рекомендуется отбирать главные компоненты по t - критерию. В других напротив - этот подход критикуется. Лично у меня наиболее адекватные результаты всегда получались при отборе по мощности. Возможно это связано с тем, что применение t - критерия корректно только в случае нормального распределения, а распределения мощных главных компонент как правило сильно отличаеются от него.

[Big_Z]

Спасибо за Ваш ответ, но, насколько я понимаю, Вы не правы. Во-первых, матрицу Х можно сразу раскладывать в svd, а результаты интерпретировать так, как я писал ранее (см., например, на стр. 90-92 здесь http://www.machinelearning.ru/wiki/images/6/6d/voron-ml-1.pdf, проверял и по другим источникам); во-вторых, доля дисперсии считается как отношение по квадратам сингулярных чисел, квадраты как раз будут собственными числами Transpose(Х)*Х. В-третьих, нормирование на с.к.о. нужно не потому, что разная дисперсия шума, а потому что сами проекции на разные признаковые переменные (то есть исходные координаты) дают разное рассеяние (дисперсию) - здесь в теории нет речи о шуме.

Какая-то беда, имхо, что матрица разреженная сильно :(

Может, есть у кого-нить ещё мысли. Прикладываю файл с данными и код Matlab для него.
Код: https://drive.google.com/file/d/0ByyhgdYYEKaZMEltMS02bkd5eFk/view?usp=sharing
Матричка: https://drive.google.com/file/d/0ByyhgdYYEKaZRV9KWkptZHFWRGs/view?usp=sharing

[Andrey_Kireew]

Действительно, можно делать и так как у Вас. В принципе результат будет точно таким же (я правда так никогда не делал).
Посмотрел для интереса Вашу матрицу. Уж не знаю что это за данные, но полезная информация в них безусловно есть - PCA находит около 50 независимых признаков. В принципе на них можно по построить регрессионное уравнение или классификатор. Конечно для этого вам всё же неплохо было бы исправить свою главную ошибку. Повторяться я не буду, тем более Вы всё равно считаете меня не правым. Могу только посоветовать посмотреть и более серьёзную литературу, а не ограничиваться обрывками сомнительных лекций из сети. Возможно разберётесь сами.

[Big_Z]

Ещё раз спасибо, буду рад, если подскажете книжку полезную.
Правильно я понимаю, что 50 признаков оказались потому, что есть 50 хорошо отделимых от нуля собств. чисел?
И что значит T = Х*U - это будет матрица координат ГК в старом признаковом базисе или это уже матрица проекций, то есть координат объектов в новом базисе - базисе, образованном ГК?

[Andrey_Kireew]

Вижу всё таки Вы нашли свою ошибку. Т - это координаты объектов в новом базисе, т. е. это и есть главные компоненты (у Воронцова, которого я успел прочитал, главными компонентами почему то называются собственные векторы ковариационной матрицы данных, но на самом деле это лишь базис).
Вам нужно брать только первые 50

T=X*U(:,1:50);

только с ними и работайте. Всё остальное это - шум. Если Вам больше нравится - называйте это остатками. Шум не обязательно физический. Это может быть шум ошибок округления, который неизбежно возникает при SVD, может быть шум выборки. И много ещё чего. Но всё это нужно отбрасывать, так как к устойчивым закономерностям оно не имеют никакого отношения.

Кстати эти 50 компонент я выделял без нормирования, так что смотрите как вам лучше. Нормирование может исказить результаты.

Что касается литературы, то всё это относится к факторному анализу. У меня на компе быстро удалось найти только 2 старые книжки: Харман. Современный факторный анализ. 1972г., К. Иберла Факторный анализ, 1980 - есть они и в сети. Сделайте в поисковике запрос "pdf djvu факторный анализ" и найдёте ещё много интересного. Обычно в таких книгах есть раздел про компонентный анализ.
Ещё можно посмотреть в википедии https://ru.wikipedia.org/wiki/%D0%9C%D0%B5%D1%82%D0%BE%D0%B4_%D0%B3%D0%BB%D0%B0%D0%B2%D0%BD%D1%8B%D1%85_%D0%BA%D0%BE%D0%BC%D0%BF%D0%BE%D0%BD%D0%B5%D0%BD%D1%82
написано вроде интересно, не не факт, что без ошибок. Такие источники нужно перепроверять.

[Big_Z]

И вновь спасибо!

[Евгений Машеров]

Интерпретировать это просто. У Вас столбцы некоррелированы. Можно предположить, что это независимые величины, а наблюдаемые отличия корреляций от нуля лишь игра выборки. Если бы имела место точная некоррелированность - все сингулярные значения были бы равны (при условии равенства дисперсий), и соответственно равны дисперсии главных компонент, а поскольку это выборка из независимых случайных величин - строгой ортогональности не будет, будут случайные отклонения коэффициентов корреляции от нуля, и в силу их разброс сингулярных значений и дисперсий. Конечно, возможна и ситуация нелинейной зависимости при некоррелированности, но независимость вероятнее.

[Big_Z]

Спасибо, Евгений Леонидович! Если интересно, эти данные - численность планктона. Под svd идёт матрица "объект-признак" Х в виде "станция сбора-вид планктона". В таблицу занесён свод результатов нерегулярных по пространству и времени рейсов за 10 лет в одном море. Матрица разреженная. Численность варьируется от 0 до 10000. Задача стояла выяснить, есть ли связи между распределением пар, троек и т.д. видов планктона. К сожалению, невозможно заполучить хорошие регулярные данные (их просто нет нигде). В итоге работаем со "средней температурой по больнице". Вот хотел рассматривать станции как объекты, а виды как признаки, а потом, выделив главные компоненты (вслед за Воронцовым, ГК я понимаю как новый базис) и рассчитав проекции на них, посмотреть положительные и отрицательные сильные корреляции столбцов исходной матрицы со столбцами матрицы проекций объектов на ГК (то есть координат в новом базисе).

В таких ситуациях, корреляции между векторами численности по станциям для пары видов не считают, а берут расстояние Брея-Кёртиса (в числителе сумма модулей разностей по всем компонентам двух векторов, а в знаменателе - сумма пар значений по всем компонентам). Для кластеризации тоже это расстояние используют, судя по многочисленным публикациям. Думал над кластеризацией (прям этих гигантских столбцов матрицы Х), но что-то сомневаюсь - слишком топорно.

Буду рад, если подскажете какие-то методы, которые можно применить для установления интересующих меня зависимостей. Или же здесь лучше признать (судя по маленьким медленно убывающим сингулярным числам) что зависимостей нет и смотреть уже более локализованные по времени и пространству подматрицы?

[Andrey_Kireew]

Посмотрите многомерное шкалирование, оно позволяет оптимально представить данные в пространстве меньшей размерности. Если нет задачи строить модель, а просто проанализировать имеющиеся данные - это лучший выбор. Кластер - анализ в Вашем случае лучше сразу не применять, при таком большом количестве признаков объекты скорее всего будут неразличимы, выбор метрики скорее всего не решит проблемы.

Посмотрите М. Дэйвисон Многомерное шкалирование - М, 1988. (должна быть в сети)

[Big_Z]

Спасибо, Andrey_Kireew, буду смотреть! Думал, что и МГК сойдёт как средство уменьшения размерности (желательно до 2 ГК :)), но, похоже, надо другими средствами.

[Евгений Машеров]

Первое, что приходит в голову - дело в разреженности. В силу которой действительные корреляции занижены.

[Andrey_Kireew]

У МГК есть ограничения. В идеале он предназначен для разделения компонент, имеющих нормальное многомерное совместное распределение. Если это не так - оптимальный результат уже не получить. Более слабое требование - это однородность анализируемых признаков. Именно в таких условиях МГК чаще всего и используется. При этом результаты получаются заведомо не оптимальные, но всё же приемлемые. Для анализа неоднородных данных, как у Вас, МГК вообще говоря не предназначен, на получение с его помощью положительного результата можно только наедятся, гарантий никаких нет.

Если Вы в matlab, то посмотрите ещё самоорганизующиеся карты Кохонена (функция learnsom - описание есть в help), они приводят к результату - похожему на многомерное шкалирование, нормальных реализаций которого в matlab нет (cmdscale - это функция линейного шкалирования, которое мало чем отличается от МГК).

[Big_Z]

Спасибо, я несколько лет назад плотно разбирался (в том числе кодил) с SOM и GHSOM - growing hierarchical SOM. Это будет та же кластеризация, то есть слишком нажористо по ресурсам и, имхо, трудно интерпретируемо: итераций много, векторы длинные, центры кластеров (веса нейронов) непонятно как задавать (конкретно в этой задаче побаиваюсь рандома) и вообще, каких размеров карту брать. Много вопросов с Кохоненом тут. Но если в Matlab-e уже есть реализация, то попробую и напишу, что получилось

[Andrey_Kireew]

Посмотрите - это визуализация Вашей матрицы:

Распределение здесь близко к биноминальному. Думаю есть смысл выполнить бинаризацию данных.

for i=1:427;
  for j=1:849;
      if M(i,j)~=0
         M(i,j)=1;
      end
 end
end

Для таких данных потом использовать соответствующие методы. Например ограниченные машины Больцмана. Что касается метрики - то можно будет использовать расстояние Хемминга.