In [1]:
# setup
from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"
import warnings
warnings.filterwarnings('ignore')
In [2]:
import os
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import pickle

Metody nienadzorowane - analiza skupień

Proces automatycznych klasyfikacji danych opiera się na dopasowaniu sygnatury (zestawu atrybutów) obiektu do wzorca. Proces tworzenia wzorców określa się uczeniem. Ze względu na sposób uczenia rozróżnia się dwa paradygmaty:

  • uczenie z nadzorem, gdzie proces uczenia, albo treningu związany jest z dopasowaniem klasyfikatora do zmiennej zależnej definiującej oczekiwane klasy lub wartości - tym zagadnieniem będziemy zajmować się w dalszej części;
  • uczenie bez nadzoru, gdzie w istniejącym zbiorze danych wyszukiwane są pewne wzorce bez odniesienia do zmiennej zależnej. Proces albo obejmuje cały zbiór danych i wynik jest ostateczną klasyfikacją - mówimy wtedy o analizie skupień; lub jeżeli analizowany zbiór stanowi jedynie podzbiór większej części tworzony jest klasyfikator, który może być zastosowany do całego zbioru.

W tej części kursu zajmiemy się mechanizmami nienadzorowanymi i zapoznamy się z mechanizmami tworzenia klasyfikatorów przy braku zmiennej zależnej. Zapoznamy się szczegółowo pakietem scikit-learn przeznaczonym do realizacji procesu uczenia maszynowego. Jako danych pokazowych użyjemy zbioru "Powiaty" użytego już w poprzednich częściach kursu. W tej części nie będziemy ponownie wczytywać danych bezpośrednio z bazy a wykorzystamy poznany już mechanizm pickle() to wczytania wcześniej przygotowanej dataFrame.

Do analizy skupień możemy używać jedynie danych numerycznych - wartości. W związku z tym przygotujemy zbiór o nazwie raw, zawierających wyłącznie kolumny spełniające to założenie.

In [3]:
os.chdir('/home/jarekj/Dropbox/dydaktyka/programowanie4/powiatData')
dt = pickle.load(open("dt.p","rb+"))
raw = dt.loc[:,["m","k","aprod","prod","pprod","malzenstwa","przyrost"]]

Prosta analiza skupień

W celu pokazania mechanizmu działania pakietu scikit-learn przeprowadzimy prostą (naiwną) analizę skupień przy użyciu popularnego algorytmu KMeans. Dane użyte w analizie nie zostaną w żaden sposób przekształcone (poza wcześniejszym usunięciem danych nienumerycznych). Scikit learn w najprostszej formie oferuje dwa uzupełniające się mechanizmy:

  • fit() - który służy do zbudowania modelu/klasyfikatora. Efektem działania tej funkcji jest obiekt zaiwrajacy model. W zależności od algorytmu na podstawie modelu można wywołać funkcje zwracające własności klasyfikatora
  • predict() - który służy do zastosowania klasyfikatora do nowych lub tych samych danych. Efektem działania funkcji predict() jest tablica numpy zawierające numery (indeksy) klas do których należą poszczególne obiekty
In [4]:
from sklearn.cluster import KMeans
kmeans1 = KMeans(5).fit(raw)
klasa = kmeans1.predict(raw)
kmeans1.cluster_centers_
Out[4]:
array([[ 3.10260421e+04,  3.21942069e+04,  9.72285441e+03,
         4.24188123e+04,  1.10785824e+04,  4.84061303e+00,
        -7.38314176e-01],
       [ 1.60347444e+05,  1.77370222e+05,  4.75820000e+04,
         2.23329778e+05,  6.68058889e+04,  4.40000000e+00,
        -2.22222222e-01],
       [ 6.73346095e+04,  7.15310286e+04,  2.12275524e+04,
         9.27499905e+04,  2.48880952e+04,  4.72380952e+00,
        -1.30476190e-01],
       [ 3.06817500e+05,  3.55787250e+05,  8.40875000e+04,
         4.32133250e+05,  1.46384000e+05,  4.20000000e+00,
        -1.77500000e+00],
       [ 7.91207000e+05,  9.33197000e+05,  2.40280000e+05,
         1.09752200e+06,  3.86602000e+05,  4.20000000e+00,
        -2.00000000e-01]])

Po zakończeniu procesu klasyfikacji można sprawdzić liczbę obiektów przypisanych do każdej klasy, przy pomocy mechanizmu np.unique(), który zwraca listę unikalnych wartości i ich liczbę wystąpień. Następnie warto wyniki przypisać do oryginalnej data frame. Numery klas rozpoczynają się od 0, ostateczne wyniki klasyfikacji powinny jednak rozpoczynać się od 1.

In [5]:
# policzenie klas
unique, counts = np.unique(klasa, return_counts=True)
dict(zip(unique, counts))

dt['klasa'] = klasa + 1
Out[5]:
{0: 261, 1: 9, 2: 105, 3: 4, 4: 1}

Skalowanie danych

Wiele klasyfikatorów wymaga wcześniejszej transformacji danych głównie skalowania standaryzującego (z-score) lub skalowania do pewnego przedziału (najczęściej [0-1]).

Skalowanie/standaryzowanie danych jest najczęściej wykonywanych operacji na danych przed wyliczeniem niepodobieństwa. Scikit-learn oferuje kilka narzędzi skalowania danych i normalizacji. W przypadku klasycznej standaryzacji (usunięcie wartości średniej) i podzielenie przez odchylenie standardowe mamy do dyspozycji dwie funkcje: scale(), która po prostu dokonuje transformacji jednej macierzy w drugą oraz klasę StandardScaler(), która wymaga zastosowania metod fit() i transform(), dokonać transformacji odwrotnej itp. oraz może być użyta w przetwarzaniu potokowym (o czym później w ramach kursu). Podobnie w przypadku skalowania do przedziału 0-1 możemy stosować MinMaxScaler() lub funkcję minmax_scale(). Oczywiście tylko MinMaxScaler pozwala na transformację odwrtotną i może być stosowany w przekształceniach potokowych.

Wizualizacja wielowymiarowej przestrzeni

Nie jest możliwe przedstawienie wszystkich cech (wymiarów) obiektów w postaci wykresu dwuwymiarowego. Pewną możliwością jest prezentacja w postaci poznanej już macierzy wykresów relacyjnych, takie rozwiązanie jednak pozawla tylko analizować niepodobieństwo na podstawie dwóch cech (wymiarów) na raz. Jeżeli chcemy przeanalizować niepodobieństwo obiektów na podstawie wszystkich cech należ zastosować skalowanie rozwiązania oparte o rozwijanie rozmaitości (manifold). Biblioteka scikit-learn oferuje kilka metod rozwijania rozmaitości, z której najpopularniejsza to MDS - czyli skalowanie wielowymiarowe (multidimentional skaling).

Skalowanie wielowymiarowe wykorzystuje macierz niepodobieństwa wyliczaną na podstawie oryginalnych wymiarów a następnie macierz rozkładana jest do modelu w którym odległości między punktami na płaszczyźnie 2D odpowiadają "odległościom" w przestrzeni wielowymiarowej.

Do wyliczenia macierzy niepodobieństwa użyjemy funkcji pdist() z pakietu scipy. funkcja zwraca niepodobieństwo w postaci rozwiniętej (jedna połówka), którą do pracy z MDS należy przekształcić do postaci kwadratowej (squareform).

Następny krok to zbudowanie modelu transformującego a następnie zastosować model do macierzy niepodobieństwa (fit_transform()). Funkcja transform() podobnie jak poznana wcześniej funkcja predict(). O ile funkcja predict() zwraca wynik działania modelu to funkcja transform() używa klasyfikatora do przekształcenia danych z oryginalnej przestrzeni do nowej. Dodatkowo scikit-learn pozwala przyspieszyć proces transformacji/predykcji poprzez zastosowanie funkcji fit_transform/fit_predict. Funkcje te są jednak ograniczone jedynie do sytuacji gdy predykcja/transformacja dotyczy całego zbioru danych.

In [6]:
import scipy.spatial.distance as dist
distance = dist.pdist(raw)
distance = dist.squareform(distance)

from sklearn.manifold import MDS
model = MDS(n_components=2, dissimilarity='precomputed', random_state=1,eps=1e-5,n_init=10)
punkty = model.fit_transform(distance)

plt.scatter(punkty[:,0],punkty[:,1])
Out[6]:
<matplotlib.collections.PathCollection at 0x7f77979b4cf8>

Note: Niepodobieństwo a odległość

W języku angielskim w nomenklaturze informatycznej używa się pojęcia distance jako synonimu niepodobieństwa. Wynika to z popularności miary eklidesowej, która jest odległością pomiędzy obiektami w przestrzeni wielowymiarowej (informacyjnej). Ponieważ w naukach o Ziemi pojęcie odległości między obiektami jednoznacznie kojarzy się z odległością geograficzną będziemy używać pojęć niepodobieństwo (dissimilarity) i podobieństwo (similarity) zamiast odległość.

Wizualizacja macierzy niepodobieństwa

Macierz niepodobieństwa w formie numerycznej jest trudna w analizie. Najprostszą formą wizualizacji jest zastąpienie wartości kodem barwnym. Pozwala to zmniejszyć rozmiar komórki do pojedynczego piksela. W tym celu trzeba wyliczyć macierz niepodobieństwa, uwzględniając fakt że dane muszą zostać standaryzowane (skalowane) a następnie wyliczamy kwadratową macierz niepodobieństwa (squareform).

Do wyświetlenia macierzy seaborn oferuje funkcję bezpośredniego wyświetlania macierzy heatmap() - mapę ciepła.

In [7]:
import sklearn.preprocessing as preprocessing

dtm = dt.loc[dt.woj.isin(['POMORSKIE','LUBELSKIE']),["aprod","prod","pprod","malzenstwa","przyrost","m","k","woj"]]
scaled = preprocessing.scale(dtm.drop(["woj"],axis=1))

distance = dist.pdist(scaled)
distance = dist.squareform(distance)
ax = sns.heatmap(distance,cmap="Blues")

Modyfikacja etykiet Taka wersja macierzy używa domyślnych etykiet, jeżeli chcemy zmodyfikować etykiety na właściwe należy zmodyfikować parametr yticklabels i przekazać do niego listę etykiet.

In [8]:
labels = list(map(lambda x: x[:1], dtm.woj)) 
ax = sns.heatmap(distance,cmap="Blues",yticklabels=labels)

Transformacja zmiennych - nie tylko skalowanie

Zmienne w zbiorze danych mogą być różnego typu

  • ilorazowe (skala ilorazowa) występuje wtedy, gdy gdy relacje pomiędzy jej wartościami mają interpretację w świecie rzeczywistym (np temperatura w skali K).
  • interwałowe (skala ilorazowa) występuje wtedy, gdy różnica pomiędzy jej wartościami dają się obliczyć i ma interpretację w świecie rzeczywistym ale relacje (dzielenie dwóch wartości przez siebie) nie mają sensu ze względu na umowne położenie zera - wymagają przeskalowania do wartości bezwzględnych na przykład temperatura w stopniach Celsjusza.
  • wskaźnikowe np. prędkość wyrażona w m/s. Niektóre operacje dla takich danych są nieprawidłowe (np mnożenie dwóch prędkości nie ma sensu)
  • komplementarne poszczególne zmienne są składową innej większej zmiennej - na przykład liczba ludności w wielu przed-, produkcyjnym i po-produkcyjnym.

Zestaw danych w zbiorze powiaty to 7 zmiennych. Trzy z nich to typowe zmienne komplementarne sumujące się do łączne liczby ludności. Ponieważ suma to liczba ludności w powiecie należy zamienić je do postaci normalizowanej dzieląc przez łączną liczbę mieszkańców, natomiast komplementarne dane liczba kobiet i mężczyzn lepiej zmienić na wskaźnik feminizacji (liczba kobiet do mężczyzn). W ten sposób z charakterystyki powiatów wyeliminujemy wpływ liczby ludności w powiecie.

In [9]:
dtg = raw.copy()
dtg['sum'] = dtg['aprod']+dtg['prod']+dtg['pprod']
dtg['aprod'] = dtg['aprod']/dtg['sum']
dtg['prod'] = dtg['prod']/dtg['sum']
dtg['pprod'] = dtg['pprod']/dtg['sum']
dtg['fem'] = dtg['k']/dtg['m']

transf = dtg[["aprod","prod","pprod","fem","malzenstwa","przyrost"]]
pickle.dump(transf,open("transf.p","wb+"))
In [10]:
sns.boxplot(data=transf)
Out[10]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f7793672438>

Pomimo wstępnego przetworzenia dane nadal znajdują się w różnych przedziałach co będzie wpływać na wynik skalowania. Z tego powodu należy je standaryzować do ujednoliconego zakresu wartości. Do wizualizacji zastosujemy tym razem stripplot(), który ilustruje rzeczywisty rozkład każdej zmiennej.

In [11]:
scaled = preprocessing.scale(transf) # preprocessing
transf_scaled = pd.DataFrame(scaled,columns=transf.columns) # przeniesienie nazw kolumn
sns.stripplot(data=transf_scaled)
Out[11]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f7793534240>

W następnym kroku można zrealizować procedurę wielowymiarowego skalowania danych

In [12]:
distance = dist.pdist(transf_scaled)
distance = dist.squareform(distance)
punkty = MDS(n_components=2, dissimilarity='precomputed', random_state=1,eps=1e-5,n_init=10).fit_transform(distance)

plt.scatter(punkty[:,0],punkty[:,1],s=10)
Out[12]:
<matplotlib.collections.PathCollection at 0x7f77934a1cf8>

Skalowane dane możemy następnie użyć do grupowania przy użyciu poznanego już algorytmu KMeans. Jednocześnie wraz z przyrostem doświadczenia możemy uproszczać zapis - w tym przypadku cały proces wyznaczania kas ograniczymy do jednej linii kodu: utworzeniu obiektu KMeans i zastosowaniu jego metody fit_predict(), który zwróci etykiety klas bez stosowania kroków pośrednich. Wyniki zostaną przypisane do bazowej data Frame i zaprezentowane przy pomocy poznanej już funkcji grupowania. Należy pamiętać że do wyników klasyfikacji należy dodać 1 (wyniki rozpoczynają się od 0).

In [13]:
klasa = KMeans(5).fit_predict(transf_scaled)
dt['klasa'] = klasa+1
dt.groupby("klasa").agg({"klasa":"count","malzenstwa":"mean","przyrost":"mean"})
Out[13]:
klasa malzenstwa przyrost
klasa
1 61 4.347541 -2.406557
2 83 4.751807 -1.104819
3 89 5.317978 1.711236
4 71 4.835211 -2.678873
5 76 4.523684 0.798684

Wyniki klasyfikacji można dodać jako nową kolumnę do istniejącej bazy danych przy pomocy funkcji opracoanej we wcześniejszej części kursu.

In [14]:
import sqlite3
def addColumn(frame,column,dbase,table,datatype='integer',fid='ogc_fid'):
    tmp = frame[[fid,column]]
    conn=sqlite3.connect(dbase)
    frame.to_sql('tmp', conn, if_exists='replace', index=True) #1
    conn.execute('ALTER TABLE {0} ADD COLUMN {1} {2}'.format(table,column,datatype)) #2
    query = 'UPDATE {0} SET {1} = (SELECT {1} FROM tmp WHERE tmp.{2} = centroids.{2}) WHERE {1} IS NULL'.format(table,column,fid)
    conn.execute(query) #3
    conn.commit() #3 
    conn.execute('drop table tmp') #4
In [15]:
#addColumn(dt,"klasa","database.db","centroids")

Gaussowskie modele mieszane

Gaussowskie modele mieszane (GMM) są popularną metodą klasyfikacji, która nie wymaga wyliczenia macierzy odległości. Jest to cała grupa metod, których podstawą jest założenie że skupienia w zbiorach danych są efektem zmieszania kilku zbiorów o naturalnym "normalnym" (gausowskim) rozkładzie. Tym samym rozpoznanie skupień odbywa się poprzez dopasowanie do fragmentów zbioru danych tylu modeli rozkładu normalnego ile zakładamy skupień.

Wykonanie klasyfikacji nie różni się od innych metod pakietu scikit-learn, wymaga dopasowania do zbioru uczącego (fit()) a a następnie predykcji na całym zbiorze danych lub zastosowania metody skróconej fit_predict(), jeżeli cały zbiór jest jednocześnie zbiorem uczącym.

In [16]:
from sklearn import mixture
global_gmm = mixture.GaussianMixture(n_components=5, covariance_type='full').fit(transf_scaled) # fit bo chcemy inne parametry
classes = global_gmm.predict(transf_scaled)
dt.assign(klasagauss = classes+1).groupby("klasa").agg({"klasagauss":"count"})
Out[16]:
klasagauss
klasa
1 61
2 83
3 89
4 71
5 76

Podobnie jak KMeans - GMM wymaga podania ścisłej liczby klas. W każdym systemie klasyfikacji nienadzorowanych optymalna liczba klas to taka, która minimalizuje pewien parametr będący miarą jakości (separacji) otrzymanych skupień. GMM posiada własne parametry oceny, są to Bayesowskie kryterium Informacyjne (BIC - Bayesian Information Criterion) i kryterium Aikena (AIC - Aiken Information Criterion). Optymalna liczba klas to taka, gdzie to kryterium osiąga (lokalne) minimum. Do ich wyznaczenia używamy odpowiednio metod bic() i aic() wyliczając jego wartość dla kolejno zmieniającej się liczby klas. Ostatnim krokiem jest wyrysowanie zmiany kryterium w formie wykresu liniowego.

In [17]:
gmm_aic=[]
search = range(1,10)
for n in search:
    m = mixture.GaussianMixture(n_components=n, covariance_type='full').fit(transf_scaled)
    gmm_aic.append(m.aic(transf_scaled))   

plt.plot(search,gmm_aic)
Out[17]:
[<matplotlib.lines.Line2D at 0x7f7791431c88>]
In [18]:
gmm_bic=[]
search = range(1,10)
for n in search:
    m = mixture.GaussianMixture(n_components=n, covariance_type='full').fit(transf_scaled)
    gmm_bic.append(m.bic(transf_scaled))   

plt.plot(search,gmm_bic)
Out[18]:
[<matplotlib.lines.Line2D at 0x7f77913efe48>]

GMM określa dla każdego obiektu prawdopodobieństwo przynależności do danej klasy. Informację na temat prawdopodobieństwa przynależności do każdego skupienia wyznaczamy metodą predict_proba(). Metoda ta występuje w wielu algorytmach klasyfikacji nadzorowanych i nienadzorowanych i jest podstawą określenia jej jakości.

In [19]:
global_gmm = mixture.GaussianMixture(n_components=5, covariance_type='full').fit(transf_scaled) # fit bo chcemy inne parametry
classes = global_gmm.predict(transf_scaled)
probas = global_gmm.predict_proba(transf_scaled)

Jeżeli chcemy określić z jakim prawdopodobieństwem dany obiekt został zaliczony do danej klasy, należy wybrać maksymalną wartość dla każdego wiersza macierzy prawdopodobieństw, przy pomocy funkcji np.max, po osi kolumn (axis=1). Taki wektor prawdopodobieństw przypisuje dla każdego obiektu jakość klasyfikacji. Ten parametr można wizualizować w postaci map podobnie jak klasy.

In [20]:
dt['gaussklasa'] = classes+1
dt['gaussprobas'] = np.max(probas,axis=1)
dt.groupby('gaussklasa').agg({'gaussklasa':'count','gaussprobas':['mean','std']})
Out[20]:
gaussklasa gaussprobas
count mean std
gaussklasa
1 86 0.789434 0.197754
2 82 0.812880 0.187041
3 55 0.739853 0.192577
4 79 0.804338 0.186758
5 78 0.907684 0.169433
In [21]:
sns.stripplot(x="gaussklasa",y="gaussprobas",data=dt.loc[:,["gaussklasa","gaussprobas"]])
Out[21]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f77913ff710>

Propagacja powinowatości (affinity propagation)

Propagacja powinowatości, to metoda która wykorzystuje tzw. przekazywanie wiadomości pomiędzy obiektami. Jej celem jest wyznaczenie tzw egzemplarów - naturalnych liderów: obiektów które lepiej sprawdzają się jako egzemplary dla innych niż jako obiekty posiadające egzemplary. Wykorzystuje macierz powinowatości - w praktyce macierz podobieństwa pomiędzy obiektami. W przypadku niepodobieństwa euklidesowego, macierz powinowatości to odwrotność macierzy niepodobieństwa: -distance. Implementacja w scikit-learn może korzystać zarówno z wbudowanej metryki euklidesowej jak i macierzy powinowatości wyliczonej zewnętrznie - na przykład w pakiecie scipy.spatial.distance. W przykładzie użyjemy gotowej macierzy powinowatości. Niestety to rozwiązanie nakłada ona pewne ograniczenia na model: nie można użyć metody predict() co oznacza że praca z dużymi zbiorami danych jest nieco skomplikowana.

Algorytm nie wymaga określenia liczby skupień, wymaga natomiast określenia preferencji. Może to być pojedyncza wartość stosowana dla całego zbioru, można też ustalać preferencje indywidualnie dla każdego obiektu. Liczba skupień jest ustalana w czasie działania algorytmu na podstawie preferencji.

In [22]:
from sklearn.cluster import AffinityPropagation
distance = dist.pdist(transf_scaled)
distance = dist.squareform(distance)
affinity = -distance
pref = -1
af = AffinityPropagation(preference=pref,affinity="precomputed").fit(affinity)
af.cluster_centers_indices_.shape[0]
Out[22]:
136

Propagacja powinowatości tworzy dużą liczbę skupień, z tego powodu trudno posługiwać się minimalizacją metryki jakości skupień do wyznaczenia ich optymalnej liczby. Do wyznaczenia liczby skupień zastosujemy metodę "łokciową", która polega na znalezieniu takiego punktu w liście wartości preferencji, do którego wraz ze zwiększaniem wartości liczba skupień rośnie nieznacznie a od której liczba skupień zaczyna gwałtownie rosnąć nawet przy niewielkich zmianach. Punkt ten wyraża się przegięciem na linii wykresu (stąd metoda "łokciowa"). DO jego wyznaczenia użyjemy zakresu przeszukiwania paremetru preferencji od -20 do -1, wyliczając dla każdej wartości liczbę skupień a następnie wynik prezentujemy w postaci wykresu liniowego.

In [23]:
prefs = list(range(-20,-1))
nclusts=[]
for i in prefs:
    af = AffinityPropagation(preference=i,affinity="precomputed").fit(affinity)
    nclusts.append(af.cluster_centers_indices_.shape[0])
    
sns.lineplot(x=prefs,y=nclusts,markers=True)
Out[23]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f7791407320>

Wyliczenie statystyk podsumowujących dla egzemplarów

In [24]:
af = AffinityPropagation(preference=-7,affinity="precomputed").fit(affinity)
af.cluster_centers_indices_
Out[24]:
array([  7,  15,  51,  66, 107, 127, 129, 162, 163, 165, 233, 300, 308,
       318, 336, 379])
In [25]:
exemplars = transf_scaled.iloc[af.cluster_centers_indices_,]
exemplars.head()
Out[25]:
aprod prod pprod fem malzenstwa przyrost
7 -1.453340 -0.154035 1.113961 0.682458 -0.639678 -1.630103
15 -0.981295 0.625054 0.349274 0.151909 -0.639678 -0.821936
51 -1.445907 -1.226670 1.701428 2.068342 -1.081637 -0.481654
66 -0.740713 -1.452280 1.326881 0.072878 -0.639678 -1.332357
107 0.208531 1.832015 -1.159979 -0.571363 -0.197719 0.454119

Dendrogram

Po wyznaczeniu egzemplarów pozostajemy z dużą liczbą potencjalnych skupień. Ich liczbę możemy zredukować stosując dendrogram na wyznaczonych egzemlarach. Do wyznaczenia dendrogramu żyjemy metod z pakietu scipy.cluster.hierarchy: dendrogram do wizualizacji, linkage do zbudowania i fcluster do wyznaczenia liczby skupień. Dendrogramy mozemy stosować do dowolnych danych, nie tylko wyników propagacji powinowatości.

Do wyznaczenia skupień potrzebujemy macierzy odległości. Pakiet scipy korzysta z macierzy odległości w formie rozwiniętej (połówki macierzy) tak więc musimy wyznaczyć pod-macierz na podstawie indeksów egzemplarów a następnie przekształcić do postaci rozwiniętej. wyznaczenie podmacierzy wymaga zastosowania funkcji np.ix_ do zbudowania prawidłowego systemu indeksowania (macierz[indexy,indexy] nie zadziała). Następnie przy pomocy funkcji linkage budujemy strukturę drzewiastą i wizualizujemy go funkcją dendrogram()

In [26]:
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
In [27]:
subdistance = distance[np.ix_(af.cluster_centers_indices_,af.cluster_centers_indices_)]
subdistance_vect = dist.squareform(subdistance)
coml_linkage = linkage(subdistance_vect,'complete')
dn = dendrogram(coml_linkage)

Dendrogram, którego etykiety to jedynie numery obiektów jest trudny w interpretacji. Dlatego dendrogram warto wyposażyć w etykiety, które pobierzemy z oryginalnej data Frame. Dendrogram ma wiele parametrów modyfikujących jego wygląd, tak więc oprócz dodania etykiet, dokonamy ich rotacji oraz zmienimy próg koloru w celu przetestowania potencjalnie innej liczby skupień.

In [28]:
labels = list(dt.powiat[af.cluster_centers_indices_])
dn = dendrogram(coml_linkage,labels=labels,leaf_font_size=8, leaf_rotation=90, 
                color_threshold = 3,get_leaves=True,truncate_mode='level')
In [29]:
clusters = fcluster(coml_linkage,3.5,'distance')
In [30]:
clusters
Out[30]:
array([3, 3, 2, 3, 1, 3, 1, 1, 1, 2, 1, 3, 1, 1, 1, 1], dtype=int32)

Do wyznaczenia ostatecznej liczby skupień musimy przeklasyfikować etykiety będące wynikiem działania propagacji powinowatości (długa lista, wiele klas) na etykiety będące efektem działania dendrogramu (krótka lista, kilka klas). Zadanie to wydaje się skomplikowane, jednak w praktyce wystarczy przeindeksować listę etykiet z dendrogramu (określoną tu jako clusters) listą etykiet z propagacji powinowatości (tu: _af.lables__). Wykonujemy to prostym ale nieoczywistym działaniem, polegającym na użyciu orginalnych numerów jako indeksów. Uwaga! skupienia muszą się zaczynać od 0.

In [31]:
labels = clusters[af.labels_]
dt['afklasa'] = clusters[af.labels_]

Ocena jakości klasyfikacji

Do oceny jakości klasyfikacji zarówno nadzorowanych jak i nienadzorowanych służą metryki. W przypadku klasyfikacji nienadzorowanych, gdzie nie ma zmiennej zawierającej klasy referencyjne scikit-learn dostarcza tylko jedną metrykę - sylwetkę (silhouette). Można wyliczyć wartość sumaryczną metryki (silhouette.score) lub wyliczyć wartość metryki dla każdego z obiektu z osobna.

Wartość metryki dla każdego obiektu z osobna pozwala zbudować diagram sylwetkowy. Wartość metryki oznacza relacje podobieństwa obiektu do własnej klasy do pierwszej klasy pod względem podobieństwa do której dany obiekt nie należy. Wartość ujemna oznacza że obiekt nie należy do tej klasy do której powinien.

In [32]:
from sklearn import metrics
score = metrics.silhouette_score(distance,labels = labels, metric = "precomputed")
sampscore = metrics.silhouette_samples(distance,labels = labels, metric = "precomputed")
score
Out[32]:
0.25848074222411915

Do wykonania diagramu:

  1. budujemy data frame z etykiet klas i wartości metryki
  2. sortujemy wartości wg klas i podobieństwa
  3. Dodajemy zmienną porządkową (po posortowaniu)
  4. Rysujemy w formie diagramu liniowego (barplot rysuje się zbyt długo)
In [33]:
sh = pd.DataFrame({'class':labels+1,'scores':sampscore})
sh_sorted = sh.sort_values(['class','scores'],ascending=[True,False])
sh_sorted['id']=np.arange(1,len(sh_sorted)+1)
sns.lineplot(x='id', y='scores', hue='class', data=sh_sorted)
plt.axhline(0, color='k',linewidth=2)
Out[33]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f77936bbcf8>
Out[33]:
<matplotlib.lines.Line2D at 0x7f77936bb518>

Analiza skupień w przestrzeni składowych głównych

Porównajmy wykresy skalowania wielowymiarowego dla danych surowych skalowanych a następnie i przekształconych do przestrzeni składowych głównych. Użyjemy do tego poznanej metody subplot.

In [34]:
from sklearn.decomposition import PCA 
from sklearn.preprocessing import scale
distance = dist.pdist(raw)
distance = dist.squareform(distance)
punkty = MDS(n_components=2, dissimilarity='precomputed', random_state=1,eps=1e-5,n_init=10).fit_transform(distance)

pca_raw = PCA(6,whiten=True).fit(raw)
pca_raw.explained_variance_ratio_
pca_components = pca_raw.transform(raw)
Out[34]:
array([9.98356701e-01, 1.36338704e-03, 2.11751664e-04, 6.81603012e-05,
       1.97566827e-10, 1.14546477e-11])
In [35]:
scaled = scale(raw)
distance = dist.pdist(scaled)
distance = dist.squareform(distance)
punkty = MDS(n_components=2, dissimilarity='precomputed', random_state=1,eps=1e-5,n_init=10).fit_transform(distance)

pca_scaled = PCA(6,whiten=True).fit(scaled)
pca_scaled.explained_variance_ratio_
pca_components_scales = pca_scaled.transform(scaled)
Out[35]:
array([7.16477118e-01, 1.91805435e-01, 8.87303094e-02, 2.52861192e-03,
       4.21719126e-04, 3.68059525e-05])
In [36]:
fig, ax = plt.subplots(2,2)
ax[0,0].scatter(pca_components[:,0],pca_components[:,1])
ax[0,1].scatter(punkty[:,0],punkty[:,1])
ax[1,0].scatter(pca_components_scales[:,0],pca_components_scales[:,1])
ax[1,1].scatter(punkty[:,0],punkty[:,1])
Out[36]:
<matplotlib.collections.PathCollection at 0x7f77912934a8>
Out[36]:
<matplotlib.collections.PathCollection at 0x7f7791293978>
Out[36]:
<matplotlib.collections.PathCollection at 0x7f77912933c8>
Out[36]:
<matplotlib.collections.PathCollection at 0x7f7791293d30>

A następnie dane po transformacji:

In [37]:
transf_scaled = scale(transf)

distance = dist.pdist(transf_scaled)
distance = dist.squareform(distance)
punkty = MDS(n_components=2, dissimilarity='precomputed', random_state=1,eps=1e-5,n_init=10).fit_transform(distance)

pca_transf = PCA(6,whiten=True).fit(transf_scaled)
pca_transf.explained_variance_ratio_
pca_components = pca_transf.transform(transf_scaled)
Out[37]:
array([5.88167048e-01, 1.68889911e-01, 1.27865726e-01, 9.86808491e-02,
       1.63964658e-02, 1.87235629e-31])
In [38]:
fig, ax = plt.subplots(1,2)
ax[0].scatter(pca_components[:,0],pca_components[:,1])
ax[1].scatter(punkty[:,0],punkty[:,1])
Out[38]:
<matplotlib.collections.PathCollection at 0x7f77912089b0>
Out[38]:
<matplotlib.collections.PathCollection at 0x7f7791208e80>
In [39]:
pc_dt = pd.DataFrame(pca_components,columns=["PC"+str(i) for i in range(1,7)])
fig, ax = plt.subplots(1,2)
sns.boxplot(data=pc_dt,ax=ax[0])
sns.boxplot(data=transf_scaled,ax=ax[1])
Out[39]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f77911b6898>
Out[39]:
<matplotlib.axes._subplots.AxesSubplot at 0x7f7791159160>

Po przekształceniu do przestrzeni składowych głównych możemy usunąć zbędne wymiary, wyjaśniające niewielką część zmienności a następnie możemy przeprowadzić dowolną procedurę klasyfikacyjną używając składowych zamiast oryginalnych zmiennych.

In [40]:
gmm_aic=[]
for n in range(1,15):
    m = mixture.GaussianMixture(n_components=n, covariance_type='full').fit(pca_components[:,0:3])
    gmm_aic.append(m.aic(pca_components[:,0:3]))   

plt.plot(range(1,15),gmm_aic)
Out[40]:
[<matplotlib.lines.Line2D at 0x7f77910c9a58>]

Przetwarzanie potokowe (pipeline)

Poprzednie przykłady pokazały, że procedury analityczne składają się z reguły z kilku kroków, gdzie wyniki poprzedniego kroku są wartościami wejściowymi w następnym. Taki zestaw procesów nosi nazwę potoku (pipeline) i może w sposób znaczący ułatwić pracę, zwłaszcza w sytuacjach gdzie procedury są powtarzalne. W niniejszym przykładzie omówioną wyżej procedurę klasyfikacyjną składającą się z trzech osobnych pod-procedur zastąpimy jedną procedurą składającym się z trzech następujących kroków: skalowania, składowych głównych i klasyfikatora GMM. Następnie te trzy procedury zostaną połączone w jeden potok. Każdy krok potoku składa się z dwuelementowej krotki: nazwy kroku i wywoływanej w danym kroku funkcji. Parametry poszczególnych funkcji można określić na etapie definiowania funkcji. dict(nazwa_funkcji__parametr)

In [41]:
from sklearn.pipeline import Pipeline
from sklearn import decomposition
from sklearn import mixture
from sklearn import preprocessing

scale = preprocessing.StandardScaler()
gmm = mixture.GaussianMixture(n_components=6, covariance_type='full')
pca = decomposition.PCA(n_components=2)

pipe = Pipeline(steps=[('scale',scale),('pca',pca),('gmm',gmm)])

Gdy potok jest już przygotowany, możemy go używać jako pojedynczego klasyfikatora z wykorzystaniem funkcji fit(), predict(), transform(), predict_proba() itp.

Wewnętrznie wszystkie kroki potoku z wyjątkiem ostatniego zwracają wynik metody transform(). Jeżeli z jakiegoś powodu potrzebujemy stanu pośredniego (na przykład w celu wyliczenia wartości miar jakości klasyfikacji) możemy wywołać pośredni krok potoku przy pomocy funkcji named_steps['krok'].metoda()

In [42]:
model = pipe.fit(raw)
classes = model.predict(raw)+1

# dostęp do pośrednikch kroków
trans = model.named_steps['pca'].transform(raw)
metrics.silhouette_score(trans, classes)
Out[42]:
0.08148975179288753

Jeżeli parametry jakiegoś kroku chcemy zmienić po zdefiniowaniu potoku można przekazać je do potoku w postaci słownika:

dict(nazwa_funkcji__parametr=wartosc)

Uwaga: słownik definiujemy jawnie przez funkcję dict! Na przykład jeżeli chcemy zmienić ilość komponentów w modelu GMM do 4, możemy to zrobić w następujący sposób:

In [43]:
model = pipe.fit(raw,dict(gmm__n_components=4))