Машинное обучение с Python : кластерный анализ

 В этой статьей разберем несколько самых популярных методов кластеризации и методов оценки их эффективности. Для начала сгенерируем тестовый набор данных и рассмотрим способы его  визуализации.  Для создания тестового набора данных используем функцию make_blobs. Задаем количество элементов 100 и количество кластеров 4. Каждый элемент имеет два показателя для того, чтобы могли реализовать двумерный график. Кластеры имеют разные стандартное отклонение, то есть некоторые кластеры более рассеиваются, чем другие. Функция также возвращает метки кластеров, они понадобятся, чтобы проверить качество работы алгоритма. 

import mglearn
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
%matplotlib inline

from sklearn.datasets import make_blobs

blobs = make_blobs(n_samples=100, 

                  n_features=2, 

                  centers=4, 

                  cluster_std=[1,1.5, 2,2], 

                  shuffle=False, 

                  random_state=7)

print("{},{},{},\n{}".format(type(blobs),len(blobs),blobs[0][:5],blobs[1]))

<class 'tuple'>,2,[[ 0.46546494  3.12315514]
 [-3.83396959  8.36956175]
 [ 1.7373078   4.42546234]
 [ 1.1312175   4.68194985]
 [ 2.20656076  5.50616718]],
[0 3 0 0 0 0 2 3 0 3 3 3 3 3 3 1 1 2 2 1 0 3 2 1 0 2 2 0 1 1 1 3 1 1 2 0 3
 1 3 2 0 2 3 2 2 3 1 2 0 0 0 1 2 2 2 3 3 1 1 3 3 1 1 0 1 3 2 2 1 0 3 1 0 3
 0 0 2 2 1 1 1 3 2 0 1 2 1 1 0 0 0 2 0 2 2 3 3 2 3 0]

В качестве объекта мы получили тьюпл, состоящий из двух элементов. Первый элемент — это непосредственно набор данных, список описаний наших объектов. И второй элемент — это метки классов. Мы видим, что каждый объект представляется списком из двух элементов. Первый элемент — это x координата, второй элемент — это y координата. И метки классов, проставляются как цифры 0,1,2,3.

print("features: {}".format(blobs[0][:10]))
print("target: {}".format(blobs[1][:10]))

features: [[ 0.46546494  3.12315514]
 [-3.83396959  8.36956175]
 [ 1.7373078   4.42546234]
 [ 1.1312175   4.68194985]
 [ 2.20656076  5.50616718]
 [ 0.58894326  4.00148458]
 [-3.24935538  6.73801217]
 [ 0.65058584  8.0105625 ]
 [-0.73000011  6.25456272]
 [-1.98661945  7.35670166]]
target: [0 3 0 0 0 0 2 3 0 3]

Перенесем эти элементы в два объекта, X и y

X,y=blobs
print("{},{},{},{}".format(type(X),X.shape,type(y),y.shape))

class 'numpy.ndarray'>,(100, 2),<class 'numpy.ndarray'>,(100,)


И сделаем из них DataFrame

df_blobs = pd.DataFrame(
{
'x1': X[:,0],
'x2': X[:,1],
'y': y
}
)
df_blobs.head()

        x1            x2        y

0 -8.474725 3.843652 0
1 -7.456176 6.198874 0
2 -9.099263 5.426828 0
3 -7.968535 5.337019 0
4 -8.716583 4.145134 0

При визуализации будем также использовать библиотеку mglearn

Определим следующую функцию для графического представления кластеров в двумерной плоскости. Удобная функция, которая выдает название графика и количество кластеров. Каждая точка данных будет отмечена в соответствии с ее меткой.

def plot_2d_clusters(x, y, ax):
    y_uniques = pd.Series(y).unique()
    for y_unique_item in y_uniques:
        x[y == y_unique_item].plot(
        title=f'{len(y_uniques)} Clusters',
        kind='scatter',
        x='x1', y='x2',
        marker=f'${y_unique_item}$',
        ax=ax,
)

Рассмотрим несколько вариантов визуализации, которые будем использовать в дальнейшем

Вариант 1

from matplotlib.colors import ListedColormap
plt.figure (figsize = (8, 8)) 
colors = ListedColormap(['red', 'yellow','blue','green'])
plt.scatter(X[:, 0], X[:, 1], c = y, cmap = colors)

В этом варианте используем функцию scatter. Эта функция позволяет нам вывести точки, зная их x и y координаты. На входе нужно подать два списка: первый — список x координат, второй — список y координат.  Дальше нам нужно указать, какими цветами рисовать объекты. Для этого используем объект colormap (цветовые карта), который создается функцией ListedColormap.   Создаем и говорим, что у нас есть всего 4 цвета, в данном случае красный, желтый, синий и зеленый ,при этом объекты с меткой 0 будут иметь красный цвет, с меткой 1 — желтый, 2 - синий, 3 - зеленый.

Вариант 2 с помощью определенной выше функции

fig, ax = plt.subplots(1, 1, figsize=(10, 6))
x, y = df_blobs[['x1', 'x2']], df_blobs['y']
plot_2d_clusters(x, y, ax)

Вариант 3 с помощью пакета mglearn

mglearn.discrete_scatter(X[:, 0], X[:, 1], y)
plt.legend(["cluster 0", "cluster 1", "cluster 2","cluster 3"], loc='best')

plt.xlabel("feature 0")
plt.ylabel("feature 1")

Кластеризация k-средних – один из самых простых и наиболее часто используемых алгоритмов кластеризации Применить алгоритм k-средних, воспользовавшись библиотекой scikit-learn, довольно просто. Здесь мы применяем его к синтетическим данным, которые использовали для построения предыдущих графиков. Мы создаем экземпляр класса KMeans и задаем количество выделяемых кластеров.

Для кластеризации используем метод K-means, задаем количество кластеров = 5

from sklearn.cluster import KMeans
kmeans = KMeans(n_clusters=5, random_state=7)
x, y = df_blobs[['x1', 'x2']], df_blobs['y']
y_pred = kmeans.fit_predict(x)

Cравнениваем наши предсказания c оригинальным метками :

fig, axs = plt.subplots(1, 2, figsize=(14, 6))

x, y = df_blobs[['x1', 'x2']], df_blobs['y']

plot_2d_clusters(x, y, axs[0])

plot_2d_clusters(x, y_pred, axs[1])

axs[0].set_title(f'Actuals: {axs[0].get_title()}')

axs[1].set_title(f'KMeans: {axs[1].get_title()}')

Один из первоначальных четырех кластеров был разделен на два, так как мы установили K на пять. И все же, как мы определяем значение K? У нас нет другого выбора, кроме как запустить алгоритм несколько раз с разным количеством кластеров и выбрать лучший. В следующем фрагменте кода мы перебираем три различных значения для n_clusters. У нас также есть доступ к конечным центроидам, которые вычисляются для каждого кластера после того, как алгоритм сходится. Просмотр этих центроидов проясняет, как алгоритм назначил каждую точку данных своему собственному кластеру. В последней строке нашего фрагмента кода используется треугольный маркер для построения центроидов на каждом из трех графиков:

n_clusters_options = [2, 4, 6]
fig, axs = plt.subplots(1, len(n_clusters_options), figsize=(16, 6))
for i, n_clusters in enumerate(n_clusters_options):
    x, y = df_blobs[['x1', 'x2']], df_blobs['y']
    kmeans = KMeans(n_clusters=n_clusters, random_state=7)
    y_pred = kmeans.fit_predict(x)
    plot_2d_clusters(x, y_pred, axs[i])
    axs[i].plot(
    kmeans.cluster_centers_[:,0], kmeans.cluster_centers_[:,1],
    'k^', ms=12, alpha=0.75
)

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

Для количественной оценки качества кластеризации необходимо использовать внутренние метрики, такие как внутрикластерное значение SSE (искажение), чтобы сравнивать эффективность разных результатов кластеризации K-Means. Удобно то, что в случае применения библиотеки scikit-leam явно вычислять внутрикластерное значение SSE не придется, т.к. после подгонки модели КМеаns оно уже доступно через атрибут inertia.

На основе внутрикластерного значения SSE можно воспользоваться графическим инструментом - так называемым методом локтя - для оценки оптимального количества кластеров k в имеющейся задаче. Можно сказать, что с увеличением k искажение будет понижаться. Причина в том, что образцы станут располагаться ближе к центроидам, которым они назначены. Идея метода локтя заключается в том, чтобы идентифицировать значение k, при котором искажение начинает увеличиваться быстрее всего, что прояснится, если построить график искажения для разных значений k:

distortions = []
for i in range(1, 11):
    km = KMeans(n_clusters=i, 
                init='k-means++', 
                n_init=10, 
                max_iter=300, 
                random_state=0)
    km.fit(X)
    distortions.append(km.inertia_)
plt.plot(range(1, 11), distortions, marker='o')
plt.xlabel('Number of clusters')
plt.ylabel('Distortion')
plt.tight_layout()
plt.show()

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

Этот коэффициент не предполагает знания истинных меток объектов, и позволяет оценить качество кластеризации, используя только саму (неразмеченную) выборку и результат кластеризации. Сначала силуэт определяется отдельно для каждого объекта. Обозначим через a — среднее расстояние от данного объекта до объектов из того же кластера, через b — среднее расстояние от данного объекта до объектов из ближайшего кластера (отличного от того, в котором лежит сам объект). Тогда силуэтом данного объекта называется величина: s=(b-a)/max(a,b) Силуэтом выборки называется средняя величина силуэта объектов данной выборки. Таким образом, силуэт показывает, насколько среднее расстояние до объектов своего кластера отличается от среднего расстояния до объектов других кластеров. Коэффициент силуэта ограничен диапазоном от -1 до 1. Значения, близкие к -1, соответствуют плохим (разрозненным) кластеризациям, значения, близкие к нулю, говорят о том, что кластеры пересекаются и накладываются друг на друга, значения, близкие к 1, соответствуют "плотным" четко выделенным кластерам. Таким образом, чем больше силуэт, тем более четко выделены кластеры, и они представляют собой компактные, плотно сгруппированные облака точек. С помощью силуэта можно выбирать оптимальное число кластеров k (если оно заранее неизвестно) — выбирается число кластеров, максимизирующее значение силуэта. Коэффициент силуэта доступен в виде функции silhouette_samples из модуля metrics библиотеки scikit-leam Теперь, вместо того чтобы выполнять визуальное исследование кластеров, мы будем перебирать несколько значений для n_clusters и сохранять оценку силуэта после каждой итерации. Как вы можете видеть, silhouette_score принимает два параметра – точки данных (x) и предсказанные метки кластера (y_pred):

from sklearn.metrics import silhouette_score
n_clusters_options = [2, 3, 4, 5, 6, 7, 8]
silhouette_scores = []
for i, n_clusters in enumerate(n_clusters_options):
    x, y = df_blobs[['x1', 'x2']], df_blobs['y']
    kmeans = KMeans(n_clusters=n_clusters, random_state=7)
    y_pred = kmeans.fit_predict(x)
    silhouette_scores.append(silhouette_score(x, y_pred))

Мы можем просто выбрать значение n_clusters, которое дает лучший результат. Здесь мы помещаем рассчитанные баллы в фрейм данных и используем гистограмму для их сравнения:

fig, ax = plt.subplots(1, 1, figsize=(12, 6), sharey=False)
pd.DataFrame(
{
'n_clusters': n_clusters_options,
'silhouette_score': silhouette_scores,
}
).set_index('n_clusters').plot(
title='KMeans: Silhouette Score vs # Clusters chosen',
kind='bar',
ax=ax
)

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

Можно более детально рассмотреть значения силуэта для каждого элемента каждого кластера. 

import numpy as np
from matplotlib import cm
from sklearn.metrics import silhouette_samples

km = KMeans(n_clusters=4, 
            init='k-means++', 
            n_init=10, 
            max_iter=300,
            tol=1e-04,
            random_state=0)
y_km = km.fit_predict(X)
cluster_labels = np.unique(y_km)
n_clusters = cluster_labels.shape[0]
silhouette_vals = silhouette_samples(X, y_km, metric='euclidean')
y_ax_lower, y_ax_upper = 0, 0
yticks = []
for i, c in enumerate(cluster_labels):
    c_silhouette_vals = silhouette_vals[y_km == c]
    c_silhouette_vals.sort()
    y_ax_upper += len(c_silhouette_vals)
    color = cm.jet(float(i) / n_clusters)
    plt.barh(range(y_ax_lower, y_ax_upper), c_silhouette_vals, height=1.0, 
             edgecolor='none', color=color)
    yticks.append((y_ax_lower + y_ax_upper) / 2.)
    y_ax_lower += len(c_silhouette_vals)
    
silhouette_avg = np.mean(silhouette_vals)
plt.axvline(silhouette_avg, color="red", linestyle="--") 
plt.yticks(yticks, cluster_labels + 1)
plt.ylabel('Cluster')
plt.xlabel('Silhouette coefficient')
plt.tight_layout()
plt.show()

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

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

Алгомеративная кластеризация

Алгомеративная кластеризация (agglomerative clustering) относится к семейству алгоритмов кластеризации, в основе которых лежат одинаковые принципы: алгоритм начинает свою работу с того, что каждую точку данных заносит в свой собственный кластер и по мере выполнения объединяет два наиболее схожих между собой кластера до тех пор, пока не будет удовлетворен определенный критерий остановки. 

Критерий остановки, реализованный в scikit-learn – это количество кластеров, поэтому схожие между собой кластеры объединяются до тех пор, пока не останется заданное число кластеров. Есть несколько критериев связи (linkage), которые задают точный способ измерения «наиболее схожего кластера». В основе этих критериев лежит расстояние между двумя существующими кластерами. 

В scikit-learn реализованы следующие три критерия: 

    ward - метод по умолчанию ward (метод Варда) выбирает и объединяет два кластера так, чтобы прирост дисперсии внутри кластеров был минимальным. Часто этот критерий приводит к получению кластеров относительно одинакового размера. 

    average - метод average (метод средней связи) объединяет два кластера, которые имеют наименьшее среднее значение всех расстояний, измеренных между точками двух кластеров. complete метод complete (метод полной связи или метод максимальной связи) объединяет два кластера, которые имеют наименьшее расстояние между двумя их самыми удаленными точками. 

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

Применить алгоритм алгомеративной кластеризации с помощью библиотеки scikit-learn, также довольно просто. 

from sklearn.cluster import AgglomerativeClustering
agg = AgglomerativeClustering(n_clusters=4)
assignment = agg.fit_predict(X)
mglearn.discrete_scatter(X[:, 0], X[:, 1], assignment)
plt.xlabel("feature 0")
plt.ylabel("feature 1")

Есть еще один инструмент для визуализации результатов иерархической кластеризации, называемый дендрограммой (dendrogram) и позволяющий обрабатывать многомерные массивы данных. К сожалению, в scikit-learn нет инструментов, позволяющих рисовать дендрограммы. Однако их легко можно создать с помощью SciPy. По сравнению с алгоритмами кластеризации scikit-learn алгоритмы кластеризации SciPy имеют немного другой интерфейс. В SciPy используется функция, которая принимает массив данных X в качестве аргумента и вычисляет массив связей (linkage array) с записанными сходствами между кластерами. Затем мы можем передать этот массив функции SciPy dendrogram, чтобы построить дендрограмму

импортируем функцию dendrogram и функцию кластеризации ward из SciPy

from scipy.cluster.hierarchy import dendrogram, ward

linkage_array = ward(X)

dendrogram(linkage_array)

делаем отметки на дереве, соответствующие двум, трем или четырем кластерам

ax = plt.gca()

bounds = ax.get_xbound()

ax.plot(bounds, [80, 80], '--', c='k')

ax.plot(bounds, [60, 60], '--', c='k')

ax.plot(bounds, [35, 35], '--', c='k')

ax.text(bounds[1], 80, ' two clusters', va='center', fontdict={'size': 15})

ax.text(bounds[1], 60, ' three clusters', va='center', fontdict={'size': 15})

ax.text(bounds[1], 35, ' four clusters', va='center', fontdict={'size': 15})

plt.xlabel("observation index")

plt.ylabel("cluster distance")

plt.savefig('Ris10.png', dpi=300)

DBSCAN

Еще один очень полезный алгоритм кластеризации – DBSCAN (density-based spatial clustering of applications with noise – плотностный алгоритм кластеризации пространственных данных с присутствием шума). Основные преимущества алгоритма DBSCAN заключаются в том, что пользователю не нужно заранее задавать количество кластеров, алгоритм может выделить кластеры сложной формы и способен определить точки, которые не принадлежат какому-либо кластеру. 

DBSCAN работает немного медленнее, чем алгоритм агломеративной кластеризации и алгоритм k-средних, но также может масштабироваться на относительно большие наборы данных. 

DBSCAN определяет точки, расположенные в «густонаселенных» областях пространства характеристик, когда многие точки данных расположены близко друг к другу. Эти области называются плотными (dense) областями пространства характеристик. 

Идея алгоритма DBSCAN заключается в том, что кластеры образуют плотные области данных, которые отделены друг от друга относительно пустыми областями. Точки, находящиеся в плотной области, называются ядровыми примерами (core samples) или ядровыми точками (core points). 

Алгоритм DBSCAN имеет два параметра: min_samples и eps. Если по крайней мере min_samples точек находятся в радиусе окрестности eps рассматриваемой точки, то эта точка классифицируется как ядровая. Ядровые точки, расстояния между которыми не превышают радиус окрестности eps, помещаются алгоритмом DBSCAN в один и тот же кластер. 

На старте алгоритм выбирает произвольную точку. Затем он находит все точки, удаленные от стартовой точки на расстоянии, не превышающем радиус окрестности eps. Если множество точек, находящихся в пределах радиуса окрестности eps, меньше значения min_samples, стартовая точка помечается как шум (noise), это означает, что она не принадлежит какому-либо кластеру. 

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

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

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

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

Применим алгоритм DBSCAN к нашему набору данных :

from sklearn.cluster import DBSCAN

eps_options = [0.1, 1.0, 2.0, 5.0]

fig, axs = plt.subplots(1, len(eps_options) + 1, figsize=(14, 6))

x, y = df_blobs[['x1', 'x2']], df_blobs['y']

plot_2d_clusters(x, y, axs[0])

axs[0].set_title(f'{axs[0].get_title()}\nActuals')

for i, eps in enumerate(eps_options, 1):

    y_pred = DBSCAN(eps=eps, min_samples=3,

    metric='euclidean').fit_predict(x)

    plot_2d_clusters(x, y_pred, axs[i])

    axs[i].set_title(f'{axs[i].get_title()}\nDBSCAN\neps = {eps}')

Полученные кластеры для набора данных blobs помогают нам идентифицировать эффект гиперпараметра eps: Очень маленький ЭПС не позволяет образоваться никаким образцам керна. Когда eps был установлен на 0.1, почти все точки рассматривались как шум. Основные точки начали формироваться по мере того, как мы увеличивали значение eps. Однако в какой-то момент, когда eps был установлен на 0,5, два кластера были ошибочно объединены. Аналогично, значение min_samples может сделать или сломать наш алгоритм кластеризации.

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

Индекс Калински-Харабаса представляет собой отношение суммы дисперсии между кластерами и дисперсии между кластерами для всех кластеров (где дисперсия определяется как сумма квадратов расстояний): Оценка выше, когда кластеры плотные и хорошо разделены, таким образом более высокое значение Калински-Харабаса относится к модели с более определенными кластерами.

Индекс Дэвиcа-Болдуина  вычисляет компактность как расстояние от объектов кластера до их центроидов, а отделимость - как расстояние между центроидами. Индекс Дэвиcа-Болдуина должен минимизироваться для роста кластеризации, более низкий индекс Дэвиса-Болдина относится к модели с лучшим разделением между кластерами. Ноль - это наименьший возможный результат. Значения, близкие к нулю, указывают на лучшее разделение. 

from sklearn.cluster import DBSCAN

from sklearn.metrics import calinski_harabasz_score

from sklearn.metrics import davies_bouldin_score

ch = []

db = []

no = []

ei = []

for e in np.arange(1.0, 5.0, 0.5):

    dbscan = DBSCAN(eps=e, min_samples=3, leaf_size=50)

    Y = dbscan.fit_predict(X)

    ch.append(calinski_harabasz_score(X, Y))

    db.append(davies_bouldin_score(X, Y))

    no.append(np.sum(Y == -1))

    ei.append(e)

fig, axs = plt.subplots(1, 3, figsize=(15, 4))

axs[0].plot(ei, ch, '-')

axs[0].set_title('calinski-harabasz score')

axs[1].plot(ei, db, '--')

axs[1].set_title('davies-bouldin score')

axs[2].plot(ei, no, '-.')

axs[2].set_title('number of noisy points')

По всем трем показателям оптимальное значение eps=2,5.

Применим алгоритм DBSCAN кластеризации с помощью библиотеки scikit-learn к нашим данным и выведем рассмотренные выше метрики и визуализируем результат

from sklearn.cluster import DBSCAN
dbscan = DBSCAN(eps=2.5, min_samples=3, leaf_size=50)
Y = dbscan.fit_predict(X)
print("No. clusters: {}".format(np.unique(dbscan.labels_).shape))
print("No. noisy points: {}".format(np.sum(Y == -1)))
print("CH = {}".format(calinski_harabasz_score(X, Y)))
print("DB = {}".format(davies_bouldin_score(X, Y)))

No. clusters: (4,)
No. noisy points: 0
CH = 416.61835189638856
DB = 0.4019974549644243

mglearn.discrete_scatter(X[:, 0], X[:, 1], Y)
plt.legend(["cluster 0", "cluster 1", "cluster 2","cluster 3"], loc='best')
plt.xlabel("feature 0")
plt.ylabel("feature 1")