Clusterização e algoritmo K-means

Organizar dados em agrupamentos é um dos modos mais fundamentais de compreensão e aprendizado. Como por exemplo, os organismos em um sistema biologico são classificados em domínio, reino, filo, classe, etc. A análise de agrupamento é o estudo formal de métodos e algoritmos para agrupar objetos de acordo com medidas ou características semelhantes. A análise de cluster, em sua essência, não utiliza rótulos de categoria que marcam objetos com identificadores anteriores, ou seja, rótulos de classe. A ausência de informação de categoria distingue o agrupamento de dados (aprendizagem não supervisionada) da classificação ou análise discriminante (aprendizagem supervisionada). O objetivo da clusterização é encontrar estruturas em dados e, portanto, é de natureza exploratória.

A técnica de Clustering tem uma longa e rica história em uma variedade de campos científicos. Um dos algoritmos de clusterização mais populares e simples, o K-means, foi publicado pela primeira vez em 1955. Apesar do K-means ter sido proposto há mais de 50 anos e milhares de algoritmos de clustering terem sido publicados desde então, o K-means é ainda amplamente utilizado.

Fonte: Anil K. Jain, Data clustering: 50 years beyond K-means, Pattern Recognition Letters, Volume 31, Issue 8, 2010

Objetivo

  • Implementar as funções do algoritmo KMeans passo-a-passo
  • Comparar a implementação com o algoritmo do Scikit-Learn
  • Entender e codificar o Método do Cotovelo
  • Utilizar o K-means em um dataset real

Carregando os dados de teste

Carregue os dados disponibilizados, e identifique visualmente em quantos grupos os dados parecem estar distribuídos.


In [0]:
# import libraries

# linear algebra
import numpy as np 
# data processing
import pandas as pd 
# data visualization
from matplotlib import pyplot as plt 
# sys - to get maximum float value
import sys

In [117]:
# load the data with pandas
url = 'https://raw.githubusercontent.com/InsightLab/data-science-cookbook/master/2019/09-clustering/dataset.csv'
dataset = pd.read_csv(url, header=None)
dataset = np.array(dataset)

plt.scatter(dataset[:,0], dataset[:,1], s=10)
plt.show()


1. Implementar o algoritmo K-means

Nesta etapa você irá implementar as funções que compõe o algoritmo do KMeans uma a uma. É importante entender e ler a documentação de cada função, principalmente as dimensões dos dados esperados na saída.

1.1 Inicializar os centróides

A primeira etapa do algoritmo consiste em inicializar os centróides de maneira aleatória. Essa etapa é uma das mais importantes do algoritmo e uma boa inicialização pode diminuir bastante o tempo de convergência.

Para inicializar os centróides você pode considerar o conhecimento prévio sobre os dados, mesmo sem saber a quantidade de grupos ou sua distribuição.

Dica: https://docs.scipy.org/doc/numpy/reference/generated/numpy.random.uniform.html


In [0]:
def calculate_initial_centers(dataset, k):
    """
    Inicializa os centróides iniciais de maneira arbitrária 
    
    Argumentos:
    dataset -- Conjunto de dados - [m,n]
    k -- Número de centróides desejados
    
    Retornos:
    centroids -- Lista com os centróides calculados - [k,n]
    """
    
    #### CODE HERE ####
    m = dataset.shape[0]
    
    centroids = list(dataset[np.random.randint(0, m - 1, 1)])
    
    for it1 in range(k - 1):
        max_dist = -1

        for it2 in range(m):
            nrst_cent_dist = sys.float_info.max

            for it3 in range(len(centroids)):
                dist = np.linalg.norm(dataset[it2] - centroids[it3])
                # Get the distance to the nearest centroid
                if (dist < nrst_cent_dist):
                    nrst_cent_dist = dist
                    nrst_cent = dataset[it2]

            if (nrst_cent_dist > max_dist):
                max_dist = nrst_cent_dist
                new_cent = nrst_cent

        centroids.append(new_cent)

    centroids = np.array(centroids)
    ### END OF CODE ###
    
    return centroids

Teste a função criada e visualize os centróides que foram calculados.


In [119]:
k = 3
centroids = calculate_initial_centers(dataset, k)

plt.scatter(dataset[:,0], dataset[:,1], s=10)
plt.scatter(centroids[:,0], centroids[:,1], marker='^', c='red',s=100)
plt.show()


1.2 Definir os clusters

Na segunda etapa do algoritmo serão definidos o grupo de cada dado, de acordo com os centróides calculados.

1.2.1 Função de distância

Codifique a função de distância euclidiana entre dois pontos (a, b).

Definido pela equação:

$$ dist(a, b) = \sqrt{(a_1-b_1)^{2}+(a_2-b_2)^{2}+ ... + (a_n-b_n)^{2}} $$$$ dist(a, b) = \sqrt{\sum_{i=1}^{n}(a_i-b_i)^{2}} $$

In [0]:
def euclidean_distance(a, b):
    """
    Calcula a distância euclidiana entre os pontos a e b
    
    Argumentos:
    a -- Um ponto no espaço - [1,n]
    b -- Um ponto no espaço - [1,n]
    
    Retornos:
    distance -- Distância euclidiana entre os pontos
    """
    
    #### CODE HERE ####
    n = len(a)
    
    distance = 0
    for i in range(n):
        distance = distance + (a[i] - b[i])**2
    
    distance = distance**0.5
    ### END OF CODE ###
    
    return distance

Teste a função criada.


In [121]:
a = np.array([1, 5, 9])
b = np.array([3, 7, 8])

if (euclidean_distance(a,b) == 3):
    print("Distância calculada corretamente!")
else:
    print("Função de distância incorreta")


Distância calculada corretamente!

1.2.2 Calcular o centroide mais próximo

Utilizando a função de distância codificada anteriormente, complete a função abaixo para calcular o centroid mais próximo de um ponto qualquer.

Dica: https://docs.scipy.org/doc/numpy/reference/generated/numpy.argmin.html


In [0]:
def nearest_centroid(a, centroids):
    """
    Calcula o índice do centroid mais próximo ao ponto a
    
    Argumentos:
    a -- Um ponto no espaço - [1,n]
    centroids -- Lista com os centróides - [k,n]
    
    Retornos:
    nearest_index -- Índice do centróide mais próximo
    """
    
    #### CODE HERE ####
    # Check if centroids has two dimensions and, if not, convert to
    if len(centroids.shape) == 1:
        centroids = np.array([centroids])
    nrst_cent_dist = sys.float_info.max
    
    for j in range(len(centroids)):
        dist = euclidean_distance(a, centroids[j])
        if (dist < nrst_cent_dist):
            nrst_cent_dist = dist
            nearest_index = j
    ### END OF CODE ###
    
    return nearest_index

Teste a função criada


In [123]:
# Seleciona um ponto aleatório no dataset
index = np.random.randint(dataset.shape[0])
a = dataset[index,:]

# Usa a função para descobrir o centroid mais próximo
idx_nearest_centroid = nearest_centroid(a, centroids)


# Plota os dados ------------------------------------------------
plt.scatter(dataset[:,0], dataset[:,1], s=10)
# Plota o ponto aleatório escolhido em uma cor diferente
plt.scatter(a[0], a[1], c='magenta', s=30)

# Plota os centroids
plt.scatter(centroids[:,0], centroids[:,1], marker='^', c='red', s=100)
# Plota o centroid mais próximo com uma cor diferente
plt.scatter(centroids[idx_nearest_centroid,0], 
            centroids[idx_nearest_centroid,1],
            marker='^', c='springgreen', s=100)

# Cria uma linha do ponto escolhido para o centroid selecionado
plt.plot([a[0], centroids[idx_nearest_centroid,0]], 
         [a[1], centroids[idx_nearest_centroid,1]],c='orange')
plt.annotate('CENTROID', (centroids[idx_nearest_centroid,0], 
            centroids[idx_nearest_centroid,1],))
plt.show()


1.2.3 Calcular centroid mais próximo de cada dado do dataset

Utilizando a função anterior que retorna o índice do centroid mais próximo, calcule o centroid mais próximo de cada dado do dataset.


In [0]:
def all_nearest_centroids(dataset, centroids):
    """
    Calcula o índice do centroid mais próximo para cada 
    ponto do dataset
    
    Argumentos:
    dataset -- Conjunto de dados - [m,n]
    centroids -- Lista com os centróides - [k,n]
    
    Retornos:
    nearest_indexes -- Índices do centróides mais próximos - [m,1]
    """
    
    #### CODE HERE ####
    # Check if centroids has two dimensions and, if not, convert to
    if len(centroids.shape) == 1:
        centroids = np.array([centroids])

    nearest_indexes = np.zeros(len(dataset))
    
    for i in range(len(dataset)):
        nearest_indexes[i] = nearest_centroid(dataset[i], centroids)
    ### END OF CODE ###
    
    return nearest_indexes

Teste a função criada visualizando os cluster formados.


In [0]:
nearest_indexes = all_nearest_centroids(dataset, centroids)

In [126]:
plt.scatter(dataset[:,0], dataset[:,1], c=nearest_indexes)
plt.scatter(centroids[:,0], centroids[:,1], marker='^', c='red', s=100)
plt.show()


1.3 Métrica de avaliação

Após formar os clusters, como sabemos se o resultado gerado é bom? Para isso, precisamos definir uma métrica de avaliação.

O algoritmo K-means tem como objetivo escolher centróides que minimizem a soma quadrática das distância entre os dados de um cluster e seu centróide. Essa métrica é conhecida como inertia.

$$\sum_{i=0}^{n}\min_{c_j \in C}(||x_i - c_j||^2)$$

A inertia, ou o critério de soma dos quadrados dentro do cluster, pode ser reconhecido como uma medida de o quão internamente coerentes são os clusters, porém ela sofre de alguns inconvenientes:

  • A inertia pressupõe que os clusters são convexos e isotrópicos, o que nem sempre é o caso. Desta forma, pode não representar bem em aglomerados alongados ou variedades com formas irregulares.
  • A inertia não é uma métrica normalizada: sabemos apenas que valores mais baixos são melhores e zero é o valor ótimo. Mas em espaços de dimensões muito altas, as distâncias euclidianas tendem a se tornar infladas (este é um exemplo da chamada “maldição da dimensionalidade”). A execução de um algoritmo de redução de dimensionalidade, como o PCA, pode aliviar esse problema e acelerar os cálculos.

Fonte: https://scikit-learn.org/stable/modules/clustering.html

Para podermos avaliar os nosso clusters, codifique a métrica da inertia abaixo, para isso você pode utilizar a função de distância euclidiana construída anteriormente.

$$inertia = \sum_{i=0}^{n}\min_{c_j \in C} (dist(x_i, c_j))^2$$

In [0]:
def inertia(dataset, centroids, nearest_indexes):
    """
    Soma das distâncias quadradas das amostras para o 
    centro do cluster mais próximo.
    
    Argumentos:
    dataset -- Conjunto de dados - [m,n]
    centroids -- Lista com os centróides - [k,n]
    nearest_indexes -- Índices do centróides mais próximos - [m,1]
    
    Retornos:
    inertia -- Soma total do quadrado da distância entre 
    os dados de um cluster e seu centróide
    """
    
    #### CODE HERE ####
    # Check if centroids has two dimensions and, if not, convert to
    if len(centroids.shape) == 1:
        centroids = np.array([centroids])

    inertia = 0
    
    for i in range(len(dataset)):
        inertia = inertia + euclidean_distance(dataset[i], centroids[int(nearest_indexes[i])])**2
    ### END OF CODE ###
    
    return inertia

Teste a função codificada executando o código abaixo.


In [128]:
tmp_data = np.array([[1,2,3],[3,6,5],[4,5,6]])
tmp_centroide = np.array([[2,3,4]])

tmp_nearest_indexes = all_nearest_centroids(tmp_data, tmp_centroide)
if inertia(tmp_data, tmp_centroide, tmp_nearest_indexes) == 26:
    print("Inertia calculada corretamente!")
else:
    print("Função de inertia incorreta!")


Inertia calculada corretamente!

In [129]:
# Use a função para verificar a inertia dos seus clusters
inertia(dataset, centroids, nearest_indexes)


Out[129]:
2402.4442703793466

1.4 Atualizar os clusters

Nessa etapa, os centróides são recomputados. O novo valor de cada centróide será a media de todos os dados atribuídos ao cluster.


In [0]:
def update_centroids(dataset, centroids, nearest_indexes):
    """
    Atualiza os centroids
    
    Argumentos:
    dataset -- Conjunto de dados - [m,n]
    centroids -- Lista com os centróides - [k,n]
    nearest_indexes -- Índices do centróides mais próximos - [m,1]
    
    Retornos:
    centroids -- Lista com centróides atualizados - [k,n]
    """
    
    #### CODE HERE ####
    # Check if centroids has two dimensions and, if not, convert to
    if len(centroids.shape) == 1:
        centroids = np.array([centroids])
    
    sum_data_inCentroids = np.zeros((len(centroids), len(centroids[0])))
    num_data_inCentroids = np.zeros(len(centroids))
    
    for i in range(len(dataset)):
        cent_idx = int(nearest_indexes[i])
        sum_data_inCentroids[cent_idx] += dataset[i]
        num_data_inCentroids[cent_idx] += 1
    
    for i in range(len(centroids)):
        centroids[i] = sum_data_inCentroids[i]/num_data_inCentroids[i]
    ### END OF CODE ###
    
    return centroids

Visualize os clusters formados


In [131]:
nearest_indexes = all_nearest_centroids(dataset, centroids)

# Plota os os cluster ------------------------------------------------
plt.scatter(dataset[:,0], dataset[:,1], c=nearest_indexes)

# Plota os centroids
plt.scatter(centroids[:,0], centroids[:,1], marker='^', c='red', s=100)
for index, centroid in enumerate(centroids):
    dataframe = dataset[nearest_indexes == index,:]
    for data in dataframe:
        plt.plot([centroid[0], data[0]], [centroid[1], data[1]], 
                 c='lightgray', alpha=0.3)
plt.show()


Execute a função de atualização e visualize novamente os cluster formados


In [0]:
centroids = update_centroids(dataset, centroids, nearest_indexes)

2. K-means

2.1 Algoritmo completo

Utilizando as funções codificadas anteriormente, complete a classe do algoritmo K-means!


In [0]:
class KMeans():
    
    def __init__(self, n_clusters=8, max_iter=300):
        self.n_clusters = n_clusters
        self.max_iter = max_iter
    
    def fit(self,X):
        
        # Inicializa os centróides
        self.cluster_centers_ = calculate_initial_centers(X, self.n_clusters)
        
        # Computa o cluster de cada amostra
        self.labels_ = all_nearest_centroids(X, self.cluster_centers_)
        
        # Calcula a inércia inicial
        old_inertia   = inertia(X, self.cluster_centers_, self.labels_)
        self.inertia_ = old_inertia
        
        for index in range(self.max_iter):
            
            #### CODE HERE ####
            self.cluster_centers_ = update_centroids(X, self.cluster_centers_, self.labels_)
            self.labels_          = all_nearest_centroids(X, self.cluster_centers_)
            self.inertia_         = inertia(X, self.cluster_centers_, self.labels_)
            
            if (self.inertia_ == old_inertia):
                break
            else:
                old_inertia = self.inertia_
            ### END OF CODE ###
                    
        return self
    
    def predict(self, X):
        
        return all_nearest_centroids(X, self.cluster_centers_)

Verifique o resultado do algoritmo abaixo!


In [134]:
kmeans = KMeans(n_clusters=3)
kmeans.fit(dataset)

print("Inércia = ", kmeans.inertia_)

plt.scatter(dataset[:,0], dataset[:,1], c=kmeans.labels_)
plt.scatter(kmeans.cluster_centers_[:,0], 
            kmeans.cluster_centers_[:,1], marker='^', c='red', s=100)
plt.show()


Inércia =  608.6035508327781

2.2 Comparar com algoritmo do Scikit-Learn

Use a implementação do algoritmo do scikit-learn do K-means para o mesmo conjunto de dados. Mostre o valor da inércia e os conjuntos gerados pelo modelo. Você pode usar a mesma estrutura da célula de código anterior.

Dica: https://scikit-learn.org/stable/modules/generated/sklearn.cluster.KMeans


In [135]:
#### CODE HERE ####
from sklearn.cluster import KMeans as sk_KMeans

skkmeans = sk_KMeans(n_clusters=3).fit(dataset)

print("Scikit-Learn KMeans' inertia: ", skkmeans.inertia_)
print("My KMeans inertia:            ", kmeans.inertia_)


Scikit-Learn KMeans' inertia:  608.6035508327782
My KMeans inertia:             608.6035508327781

3. Método do cotovelo

Implemete o método do cotovelo e mostre o melhor K para o conjunto de dados.


In [169]:
#### CODE HERE ####

# Initialize array of Ks
ks = np.array(range(1, 11))
# Create array to receive the inertias for each K
inertias = np.zeros(len(ks))

for i in range(len(ks)):
    # Compute inertia for K
    kmeans = KMeans(ks[i]).fit(dataset)
    inertias[i] = kmeans.inertia_
    
    # Best K is the last one to improve the inertia in 30%
    if (i > 0 and (inertias[i - 1] - inertias[i])/inertias[i] > 0.3):
        best_k_idx = i

print("Best K: {}\n".format(ks[best_k_idx]))
plt.plot(ks, inertias, marker='o')
plt.plot(ks[best_k_idx], inertias[best_k_idx], 'ro')


Best K: 3

Out[169]:
[<matplotlib.lines.Line2D at 0x7f9fe0e2be48>]

4. Dataset Real

Exercícios

1 - Aplique o algoritmo do K-means desenvolvido por você no datatse iris [1]. Mostre os resultados obtidos utilizando pelo menos duas métricas de avaliação de clusteres [2].

Dica: você pode utilizar as métricas completeness e homogeneity.

2 - Tente melhorar o resultado obtido na questão anterior utilizando uma técnica de mineração de dados. Explique a diferença obtida.

Dica: você pode tentar normalizar os dados [3].

3 - Qual o número de clusteres (K) você escolheu na questão anterior? Desenvolva o Método do Cotovelo sem usar biblioteca e descubra o valor de K mais adequado. Após descobrir, utilize o valor obtido no algoritmo do K-means.

4 - Utilizando os resultados da questão anterior, refaça o cálculo das métricas e comente os resultados obtidos. Houve uma melhoria? Explique.


In [0]:
#### CODE HERE ####