Иерархический кластерный анализ на языке программирование Python

Классификация — одна из важнейших задач, встречающихся при анализе данных. В зависимости от постановки можно различать следующие задачи: кластеризацию (классификацию в отсутствии обучающей выборки) и, собственно, классификацию (при наличии обучающей выборки),  когда данные необходимо соотнести с уже известными классами.

Первая задача возникает, как правило, когда исследуются новые объекты или явления; в этом случае кластеризация позволяет выделить однородные группы объектов, и далее, планировать последующие изыскания. Классификация при наличии обучающей выборки предполагает, основываясь на уже имеющихся данных и их соответствии известным классам, определение класса для вновь поступающих данных. Распознавание текста, речи, аутентификацияавторизация по отпечатку пальца — представляют задачи классификации этого типа.

Рассмотрим более подробно задачу кластеризации. Она предполагает разбиение массива данных на однородные группы, исходя из представлений о близости элементов между собой. При этом число кластеров может постулироваться заранее, или определяться в процессе кластерного анализа. В зависимости от того, каким образом определяются расстояния между кластеризуемыми объектами, результаты классификации будут различны. Если данные представляют собой наборы количественных признаков, вполне естественно рассматривать их как точки многомерного факторного пространства, а в качестве расстояния, например, использовать евклидову метрику. В случае трехмерного факторного пространства такой подход будет иметь ясную геометрическую интерпретацию и связь с реальным трехмерным миром (каждая точка данных - некоторая точка в трехмерном пространстве; расстояние между данными - это расстояние между такими точками).

Следует помнить, что количественные признаки могут иметь различную физическую интерпретацию: один признак может отождествляться с длиной,  другой и массой объекта, поэтому понимание "расстояния" между элементами в этом случае весьма условно. Однако, если выбранные признаки полно характеризуют исследуемые объекты, то даже такие "бессмысленные" расстояния (несмотря на различие единиц измерения) имеют важное свойство: элементы, имеющие близкие значения количественных признаков, как следствие, характеризуются малыми расстояниями между друг другом, вполне ожидаемо похожи друг на друга. Таким образом, расстояния, построенные даже на базе разнородных признаков, выполняют свою главную роль — характеризуют сходство объектов.

В задаче кластеризации признаки не всегда могут иметь количественный характер. Например, наряду с различными морфометрическими характеристиками исследуемой совокупности объектов, может отмечаться наличие или отсутствие какой-либо их особенности. Два объекта могут иметь практически идентичные морфометрические признаки, но принципиально различаться ввиду этой особенности. Есть соблазн закодировать наличие особенности числом "1", а ее отсутствие, например, "0", и далее считать, что имеются только количественных признаки. Однако, в этом случае наличие или отсутствие особенности может легко маскироваться незначительными изменениями в морфометрических признаках. Грубо говоря, если ширина двух объектов отличается на 1 см, то это может быть эквивалентным тому, что объекты имеют или не имеют важную особенность. А это может быть критичным для помещения объекта в тот или иной кластер. Данный пример показывает, что с назначением качественным признакам количественных показателей нужно быть крайне осторожным.

Общая схема кластерного анализа

Иерархический кластерный анализ состоит из двух этапов: 1) определения "расстояний" между классифицируемыми элементами; 2) объединения элементов в группы по степени близости (исходя из вычисленных расстояний).

Если на первом этапе расстояния определяются между элементами, и, таким образом, он является подготовительным; то на втором — элементы объединяются в группы. На этом этапе необходимо сравнивать группы элементов по степени сходства. Для сравнения групп элементов также используется функция "расстояния", в иерархическом кластерном анализе она называется метод кластеризации ("single", "complete" и т.д.). Например, в случае метода "single" — расстоянием между двумя кластерами является минимальное расстояния между элементами, принадлежащим им.  В случае метода "complete" в роли межгруппового (межкластерного) расстояния используется расстояние между наиболее удаленными представителями двух групп (кластеров). Таким образом, на базе расстояний между элементами, строятся расстояния между группами элементов. Если алгоритм определения расстояний между группами  определен, можно начинать объединение элементов в кластера. Это делается последовательно объединяя сначала самые близкие элементы, потом самые близкие к образовавшимся группам элементов, и так далее,  пока все элементы не будут объединены в один кластер.

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

Вычислительный пример

Рассмотрим пример кластерного анализа с использованием пакета SciPy (язык программирования Python) и визуализации полученных результатов.

Данные для кластерного анализа

Исходные данные для анализа — биологического происхождения; они представляют собой 25 описаний и количество встретившихся в этих описаниях определенных типов сообществ растений (которые зашифрованы цифрами 1-9). В каждой из колонок указываются встречаемости того или иного типа сообщества в описании; фрагмент таблицы (только 6 колонок) данных имеет следующий вид:

Задача состоит в том, чтобы провести кластерный анализ данных, определить сходные по составу сообществ описания, обосновать оптимальное количество кластеров адекватно характеризующих данные.

MS Excel документ с данными можно загрузить

по ссылке (17,0 KB)

Вычислительная среда для кластерного анализа

Импортируем необходимый набор необходимых функций SciPy, Matplotlib, Scikit-Learn.

# coding: utf8  # Общие функции кластеризации from scipy.cluster.hierarchy import * from scipy.spatial.distance import pdist  # Для оценки качества кластеров как смеси распределений from sklearn import mixture  # Анализа главных компонент from sklearn.decomposition import PCA   # Линейный дискриминантный анализ from sklearn.lda import LDA   # Для загрузки данных из xls документа import pandas as pd  # Для отрисовки графиков from pylab import * 

Далее, выполним настройку – фиксируем датчик случайных чисел (для воспроизводимости результатов при повторном запуске скрипта), настроим шрифт и его размер, используемые при построении графиков в matplotlib.

# Для вопроизодимости результатов, зависящих от генератора случайных чисел np.random.seed(1000)   # Настройка шрифтов для будущих графиков rcParams['font.family'] = 'DejaVu Sans' # Понимает русские буквы rcParams['font.size'] = 16  # numpy нужен, но уже импортирован из pytlab # import numpy as np 

Для загрузки данных из xls-файла удобно использовать пакет pandas (но для чтения xls должен быть установле xlrd).  

# Загрузим данные из xls документа data = pd.read_excel('data.xls')  # Уберем последнюю букву из наименования строк, # она в кириллице и могут быть проблемы с  # пакетами программ ниже labels = map(lambda x: x[:-1], list(data['name']))   # Выделим данные, начиная с первой колонки # Это то, что подлежит анализу data_for_clust =  data.ix[:,1:].as_matrix() 

Заметим, что в силу специфики данных — они отражают встречаемости различных типов сообществ в описаниях и, таким образом, выражены в общих единицах измерения, масштабирования данных не требуется. Переход к долям встретившихся типов сообществ также внесет ошибку, так как абсолютное число сообществ для данного описания также (вероятно; см. заключение) имеет значение. Следовательно, никакой предварительной обработки данных не требуется, и можно приступать к их кластеризации.

Выполним кластерный анализ в SciPy:

# Вычислим расстояния между каждым набором данных, # т.е. строками массива data_for_clust # Вычисляется евклидово расстояние (по умолчанию) data_dist = pdist(data_for_clust, 'euclidean')  # Главная функция иерархической кластеризии # Объедение элементов в кластера и сохранение в  # специальной переменной (используется ниже для визуализации  # и выделения количества кластеров data_linkage = linkage(data_dist, method='average') 

Процедура последовательного объединения элементов в кластера, выполняемая функцией linkage из scipy (аналогично и в MatLab), предполагает графическое представление результата в виде дендрограммы. Анализ дендрограммы, как правило, позволяет оценить оптимальное число кластеров. Однако, существуют и другие подходы к оценке их количества.

Один из подходов — метод локтя, когда строится функция, отражающая изменение суммы внутрикластерных вариаций данных от количества кластеров. Число кластеров является наиболее подходящим (оптимальным) для данной задачи в той позиции, где отмечено наибольшее изменение этой функции.

Определим данную вспомогательную функцию для выполнения кластерного анализа.

def wgss(data, groups):     '''Within groups sum of squares (wgss) #     Сумма квадратов расстояний от центроида до каждой точки данных  #     в многомерном пространстве. #     Специально на английском, чтобы объяснить название функции      '''     _data = np.array(data)     res = 0.0     for cluster in groups:         inclust = _data[np.array(groups) == cluster]         meanval = np.mean(inclust, axis=)         res += np.sum((inclust - meanval) ** 2)     return res 

Рассмотрим результаты применения метода локтя для анализируемых данных:

# -------------- Elbow method (метод локтя) ------------------------- elbow = [np.nan, wgss(data_for_clust, [1]*len(data_for_clust[:,1]))] for k in range(2, 10):     groups = fcluster(data_linkage, k, criterion='maxclust')     elbow.append(wgss(data_for_clust, groups))  fig = figure() ax = fig.add_subplot('121') # 2 графика в строке, выбираем первый график elbow = np.array(elbow) # Пусть будет numpy массив, удобней... ax.plot(elbow/np.nanmax(elbow),'o', ls='solid') ax.set_xlim([, 10]) ax.set_ylim([, 1.2]) ax.set_title(u'Сумма внутригрупповых вариаций') ax.set_xlabel(u'Число кластеров')  ax1 = fig.add_subplot('122') # выбираем второй график в строке  ax1.plot((elbow[1]-elbow)/np.nanmax(elbow),'o', ls='solid')  ax1.set_xlim([, 10]) ax1.set_ylim([, 1.2]) ax1.set_xlabel(u'Число кластеров') ax1.set_title(u'Доля объясняемой вариации') # --------------------------------------------------------------------- 

На рисунке представлены два эквивалентных по назначению графика: слева – изменение суммы внутригрупповых вариаций, справа – доля внутригрупповой вариации данных относительно общей их вариации. Видно, что при количестве кластеров k = 2 наблюдается наибольшее изменение функции, а далее, ее значения несколько стабилизируются. Также небольшое изменение наблюдается при k = 6. Таким образом,  кандидатами в оптимальные значения количества кластеров являются k = 2  и k = 6.

Также для оценки числа кластеров при большом количестве данных можно использовать информационные критерии Акаике и Байеса (на базе задачи разделения смеси распределений); поскольку имеется всего лишь 25 измерений, что, очевидно, мало, чтобы рассматривать задачу разделения смеси распределений, количество компонентов которой фактически сравнимо с числом данных (до 10 компонент смеси, т.к. исследуется ситуации до 10 кластеров), критерии Акаике и Байеса здесь не используются. При необходимости (для задачи разделения смеси) они имеются в пакете scikit-learn.

Следующий этап исследования — визуализация — весьма немаловажный этап, позволяющий подтвердить адекватность оценок количества кластеров, полученных по методу локтя, а также сделать вывод об их оптимальном количестве непосредственно.

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

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

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

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

Из проекции на главные оси видно, что число кластеров равное 2 выглядит наиболее правдоподобным, кластера хорошо разделены в факторном пространстве и проекция это отражает.

Главные компоненты остаются неизменными независимо от количества кластеров — они рассчитываются для всей совокупности данных. Зная число кластеров можно использовать несколько иной подход, а именно найти такие оси (как и в случае метода главных компонент, только в рамках  линейных преобразований исходного факторного пространства), в проекциях на которые заданные кластера будут выглядеть максимально удаленными друг от друга (т.е. разделимыми). Такой метод уже имеется в арсенале у прикладных исследователей — это линейный дискриминантный анализ. Определение главных дискриминантных осей в этом случае также выполняется на основе задачи оптимизации (здесь однако  мы его рассматривать не будем, так как это тема для отдельной статьи).

# ------------- Проекция кластеров на главные компоненты -------------- pca = PCA(n_components=2) # Будем проецировать данные на 2 главные компоненты Xt = pca.fit_transform(data_for_clust) fig = figure() for k in range(2, 8):     groups = fcluster(data_linkage, k, criterion='maxclust')     ax = fig.add_subplot(3, 2, k-1)     for j, m, c in zip(range(k), 'so><^v', 'rgbmyc'):         ax.scatter(Xt[groups==(j+1), ], Xt[groups==(j+1), 1], marker=m, s=30, label='%s'%j, facecolors=c)         ax.set_title('k=%s'%k)         ax.legend(fontsize=14, loc="lower right") fig.suptitle(u'Проекрация кластеров на главные компоненты') # ---------------------------------------------------------------------   figure() dendrogram(data_linkage, labels=labels) title(u'Дендрограмма')  show()  # ------------- Проекция на главные дискриминантны оси----------------- import seaborn as sns; sns.set(color_codes=True) rcParams['font.family'] = 'DejaVu Sans' # Импорт seaborn сбрасывает настройки, устанавливаем их снова rcParams['font.size'] = 16 lda = LDA(n_components=2) # Проецируем данные на 2 главные дискримнационные оси fig = figure() for k in range(3, 9):     groups = fcluster(data_linkage, k, criterion='maxclust')     lda.fit(data_for_clust, groups)     Xt = lda.transform(data_for_clust) # Собственно проекция данных     ax = fig.add_subplot(3, 2, k-2)     for j, m, c in zip(range(k), 'so><^v', 'rgbmyc'):         # Проекции при различном числе кластеров разные (в отличие от главных компонент!)         # Поэтому и данные выглядят на графиках различно         ax.scatter(Xt[groups==(j+1), ], Xt[groups==(j+1), 1], marker=m, s=30, label='%s'%j, facecolors=c, zorder=10)         sns.kdeplot(Xt[:, ], Xt[:, 1], shade=True, cmap="Blues")         ax.set_title('k=%s'%k)         ax.legend(fontsize=14, loc="lower right") fig.suptitle(u'Проекрация кластеров на дискриминантные оси') # --------------------------------------------------------------------- 

Поскольку число компонент, определяемое по методу линейного дискриминтного анализа, не может превосходить (этот вопрос будет подробней обсуждаться в специальной статье про линейную дискриминацию) числа классов за вычетом единицы, мы рассматриваем ситуации только с 3 и более кластерами. Голубым цветом на рисунке представлена двумерная оценка плотности проекций данных на главные оси. Случай с k=3 выглядит наиболее привлекательным и также свидетельствует об оптимальном числе классов, равном 2 (кластера под номерами "1" и "2" в этом случае мало различимы и могут быть объеденины).

Пример иерархического кластерного анализа на PythonSciPy (5,8 KB)

Таким образом, в отношении оптимального числа кластеров можно сделать следующие выводы:

• Наибольшие изменения в графике внутригрупповой вариации возникает при числе кластеров равном 2;
• Проекции данных на главные компоненты подтверждают, что при k = 2 кластера наиболее разделимы;
• Дендрограмма также констатирует явное разделение описаний на две группы;

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