27 августа 2012 г.

Алгоритм 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.
исходники тут

2 комментария:

Анонимный комментирует...

При экспорте результатов кластеризации. В то время, когда заполняется список

res[l.tolist().index(l.min())].append((labels[i], data[i]))

Как количество названий(столбцов) может равняться количеству значений массива данных. Эти значения сойдутся только при переборе по второму уровню массива данных.

esemi комментирует...

Слушай, вот чесно, восемь лет прошло))
Я ен помню что где куда сходится))