Алгоритм k-means [Кластерный анализ и Python]

Наконец дошли руки дописать цикл заметок о кластерном анализе. 

Тема данной заметки алгоритм k-means и пример работы с ним на языке python.

В первой заметке уже был описан набор данных, который мы исследуем, а также приведён список литературы на тему, поэтому перейдём сразу к делу.

Теория

Данный алгоритм разбивает исходное множество объектов на k кластеров так, что средние в кластере (для всех переменных) максимально возможно отличаются друг от друга.
Выбор числа k может базироваться на результатах предшествующих исследований, теоретических соображениях или интуиции. Если нет предположений относительно этого числа, рекомендуют создать 2 кластера, затем 3, 4, 5 и т.д., сравнивая полученные результаты.

Работа алгоритма делится на несколько этапов:

• случайно выбрать k точек, являющихся начальными «центрами масс» кластеров (число k выбирается исследователем);
• отнести каждый объект к кластеру с ближайшим «центром масс»;
• пересчитать «центры масс» кластеров согласно их текущему составу;
• если критерий остановки алгоритма не удовлетворен, вернуться к п. 2.

В качестве критерия остановки работы алгоритма обычно выбирают одно из двух условий:

• «центром масс» кластеров стабилизировались, т.е. все наблюдения принадлежат кластеру, которому принадлежали до текущей итерации;
• число итераций равно максимальному числу итераций.

Данный алгоритм также имеет недостатки:

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

Практика

Сам алгоритм реализован в старой доброй библиотечке scipy, в модуле cluster.vq. Функция kmeans(obs, k_or_guess, iter=20, thresh=1e-05)

принимает на вход набор данных и количество искомых кластеров. 

from scipy.cluster import *
names, data = get_data()
centroids = vq.kmeans(numpy.array(data), 5, iter=200)[0]

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

import numpy
from scipy.spatial.distance import cdist
def kmeans_export(centroids, data, labels):
    """Export kmeans result"""
    res = [[] for i in xrange(len(centroids))]
    d = cdist(numpy.array(data), centroids, 'euclidean')
    for i, l in enumerate(d):
        res[l.tolist().index(l.min())].append((labels[i], data[i]))
    return res
K_res = kmeans_export(centroids, data, names)

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

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

from collections import deque
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
def kmeans_draw(clusters):
    """Drawing kmeans clustering result"""
    colors = deque(['r', 'g', 'b', 'c', 'm', 'y', 'k'])
    fig = plt.figure()
    # Prior to version 1.0.0, the method of creating a 3D axes was different. For those using older versions of matplotlib,
    # change ax = fig.add_subplot(111, projection='3d') to ax = Axes3D(fig).
    ax = Axes3D(fig)
    for cluster in clusters:
        color = colors.popleft()
        for name, coord in cluster:
            x, y, z = coord
            ax.plot3D([x], [y], [z], marker='o', c=color)
    ax.set_xlabel(u'Белки')
    ax.set_ylabel(u'Жиры')
    ax.set_zlabel(u'Углеводы')
    plt.show()
kmeans_draw(K_res)

Данная функция перебирает все полученные кластеры и рисует каждый входящий в него объект в виде цветной точки на трёхмерном графике.
Вместо сохранения картинки мне показалось логичным предоставить возможность покрутить график в пространстве (plt.show()) и вот, что получилось.

Послесловие 

В следующий раз пободаемся с алгоритмом MST.
исходники тут