Laboratorium Statystycznej Eksploracji Danych

LABORATORIUM 3

Ocena klasyfikatorów

Macierz pomyłek

Na dzisiejszych zajęciach zajmiemy się przede wszystkim metodami oceny jakości klasyfikatorów. Zgodnie z Wykładem 4 podstawowym obiektem, który ułatwi nam pracę jest macierz pomyłek, oznaczana na potrzeby tych ćwiczeń jako CM. Poszczególne pola tej macierzy

\(\mathbf{CM} = \left(\begin{array}{cc}TN & FP\\FN & TP\\\end{array}\right)\)

są określane jako TRUE NEGATIVE (dane: negatywne, predykcja: negatywna), FALSE POSITIVE (dane: negatywne, predykcja: pozytywna), FALSE NEGATIVE (dane: pozytywne, predyckja: pozytywna) oraz TRUE POSITIVE (dane: pozytywne, predykcja: pozytywna). Sumowanie rzędów \(\mathbf{CM}\) daje informacje o oryginalnych klasach, sumowanie kolumn - o przewidywanych. Wygodnym wskaźnikiem jakości klasyfikatora jest jego dokładność (ACCURACY) zdefiniowana jako

\(ACC = \frac{TP + TN}{TP + TN + FP + FN}\)

innymi są FALSE POSITIVE RATE (FPR, 1 - specificity), TRUE POSITIVE RATE (TPR, sensitivity), POSITIVE PREDICTIVE VALUE (PPV, precision), NEGATIVE PREDICTIVE VALUE (NPV)

\(FPR = \frac{FP}{FP + TN} = 1 - \frac{TN}{FP + TN}\)
\(TPR = \frac{TP}{TP + FN}\)
\(PPV = \frac{TP}{TP + FP}\)
\(NPV = \frac{TN}{TN + FN}\)

Dla ustalenia uwagi znów zaczynamy od nieśmiertelnego rozkładu Gaussa w 2D. Tym razem zdefiniujemy funkcję draw_data(), którą będziemy wykorzystwać do losowania obserwacji z tego rozkładu (dodatkowo od razu ładujemy potrzebne biblioteki)

import numpy as np
import pandas as pd

def draw_data(n1, n2, m1, m2, S1, S2):

    data1 = np.random.multivariate_normal(m1, S1, n1)
    data2 = np.random.multivariate_normal(m2, S2, n2)

    data = np.concatenate((data1, data2), axis = 0)
    data = pd.DataFrame(data, columns=["x", "y"])

    data["klasa"] = np.repeat([0, 1], [n1, n2])
    
    return data

Następnie losujemy obserwacje z rozkładów o zadanych parametrach (liczności \(n_1\), \(n_2\), wektory średnich \(\mathbf{m}_1\) i \(\mathbf{m}_2\) oraz macierze kowariancji \(\mathbf{S}_1\) i \(\mathbf{S}_2\)) i na takich danych, uważanych od tej chwili za PU (próba ucząca) budujemy klasyfikatory: LDA, QDA, a także naiwnego Bayesa

m1, m2 = [-1, -1], [2, 2]
S1 = [[4, 2], [2, 4]]
S2 = [[4, 2], [2, 4]]
n1, n2 = 30, 20

data = draw_data(n1, n2, m1, m2, S1, S2)

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
from sklearn.naive_bayes import GaussianNB

lda = LinearDiscriminantAnalysis()
qda = QuadraticDiscriminantAnalysis()
nb = GaussianNB()

for clf in [lda, qda, nb]:
    clf.fit(data[["x","y"]], data.klasa)

GaussianNB()

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

Teraz skonstruujemy funkcję CM_large(), której celem będzie wyznaczenie macierzy pomyłek, a następnie obliczenie istotnych, wyżej wymienionych cech, tzn. dokładności, FPR, TPR. Kolejnym krokiem będzie, znane już z poprzednich zajęć, dokonanie klasyfikacji za pomocą metody powtórnego podstawienia. Dodatkowo efekty działania funkcji CM_large() dla różnych klasyfikatorów umieścimy w jendej ramce danych tak, aby łatwo można było je porównać.

def CM_large(x, y):
    
    cm = pd.crosstab(x, y)
    
    acc = np.diag(cm).sum() / cm.to_numpy().sum()
    tnr, tpr = np.diag(cm) / cm.sum(axis = 1)
    npv, tpv = np.diag(cm) / cm.sum(axis = 0)
    
    return acc, 1 - tnr, tpr, npv, tpv

preds = pd.DataFrame(np.array([clf.predict(data[["x","y"]]) for clf in [lda, qda, nb]]).T)
preds.columns = ["lda", "qda", "nb"]  

res = pd.DataFrame([CM_large(data.klasa, preds[i]) for i in preds.columns])
res.index = preds.columns
res.columns = ["acc", "fpr", "tpr", "ppv", "npv"]
res

##       acc       fpr   tpr       ppv       npv
## lda  0.74  0.200000  0.65  0.774194  0.684211
## qda  0.78  0.166667  0.70  0.806452  0.736842
## nb   0.76  0.233333  0.75  0.821429  0.681818

Możemy też skorzystać z dedykowanej funkcji do tworzenia macierzy pomyłek

from sklearn import metrics

cm = metrics.confusion_matrix(data.klasa, preds.lda)
cm

## array([[24,  6],
##        [ 7, 13]])

lub wręcz dokonac wizualizacji:

import matplotlib.pyplot as plt

cm_display = metrics.ConfusionMatrixDisplay(cm).plot()
plt.show()

Ciekawą opcją jest też wykonanie raportu klasyfikacyjnego dla naszych modeli:

for i in preds.columns:
    print(metrics.classification_report(data.klasa, preds[i]))

##               precision    recall  f1-score   support
## 
##            0       0.77      0.80      0.79        30
##            1       0.68      0.65      0.67        20
## 
##     accuracy                           0.74        50
##    macro avg       0.73      0.73      0.73        50
## weighted avg       0.74      0.74      0.74        50
## 
##               precision    recall  f1-score   support
## 
##            0       0.81      0.83      0.82        30
##            1       0.74      0.70      0.72        20
## 
##     accuracy                           0.78        50
##    macro avg       0.77      0.77      0.77        50
## weighted avg       0.78      0.78      0.78        50
## 
##               precision    recall  f1-score   support
## 
##            0       0.82      0.77      0.79        30
##            1       0.68      0.75      0.71        20
## 
##     accuracy                           0.76        50
##    macro avg       0.75      0.76      0.75        50
## weighted avg       0.77      0.76      0.76        50

Krzywa ROC

Powyższe wartości (ACC, FPR, TPR) są punktowe, tzn. zakładają ustalony próg prawdopodobieństwa a posteriori, powodujący, że dana obserwacja zostaje przypisana do konkretnej klasy (dokładniej: próg równy 1/2). Czasem jednak istotne pytaniem staje się: jaki jest profil danego klasyfikatora, tzn. jak się zachowuje przy założeniu nierównego traktowania przynależności do klas. Taki profil często przedstawiany jest w postaci krzywej ROC (receiver-operating curve), gdzie na osi X zostają odłożone wartości FPR, natomiast na osi Y - TPR dla różnych poziomów odcięcia, czyli progów prawdopodobieństw a posteriori np. p(2|x).

W celu otrzymania krzywej ROC musimy nieco zmodyfikować funkcję odpowiedzialną za wyznaczanie wartości FPR i TPR

def CM_values(x, prob, th):
    
    cm = pd.crosstab(data.klasa, [1 if i >= th else 0 for i in y_prob])
    cm = cm.reindex(columns = np.unique(data.klasa), fill_value = 0)
    
    #print(cm)
    
    tnr, tpr = np.diag(cm) / cm.sum(axis = 1)
    
    return 1 - tnr, tpr

Pobierzemy wartości prawdopodobienstw a posteriori dla danej klasy, posortujemy je i dla każdego z nich wywołamy funkcję CM_values(). W ten sposób do powyższej funkcji przekazywana jest wartość prawdopodobieństwa, która ma być traktowana jako próg \(\tau\) - obserwacje z \(p(1|x)\) wyższą niż \(\tau\) są klasyfikowane jako pochodzące z klasy “1”.

y_prob = lda.predict_proba(data[["x","y"]])[:,1]
thr = np.unique(np.sort(y_prob))

z = np.array([CM_values(data.klasa, y_prob, t) for t in thr])

plt.plot(z[:,0],z[:,1])

Oczywiście, możemy skorzystać z gotowej funkcji RocCurveDisplay.from_estimator() i porównac ją z naszymi wynikami

metrics.RocCurveDisplay.from_estimator(lda, data[["x","y"]], data.klasa)

## <sklearn.metrics._plot.roc_curve.RocCurveDisplay object at 0x000001C7FB876BA0>

plt.scatter(z[:,0],z[:,1], color = "orange")
plt.axline((0,0), slope = 1, color = "red")

Oczywiście zawsze za pomocą tej metody możemy porównać klasyfikatory:

fig, ax = plt.subplots(figsize=(6, 6))

for clf in [lda, qda, nb]:
    metrics.RocCurveDisplay.from_estimator(clf, data[["x","y"]], data.klasa, ax = ax)

## <sklearn.metrics._plot.roc_curve.RocCurveDisplay object at 0x000001C7FBB274D0>
## <sklearn.metrics._plot.roc_curve.RocCurveDisplay object at 0x000001C7FBB88410>
## <sklearn.metrics._plot.roc_curve.RocCurveDisplay object at 0x000001C7FB9A9CD0>

    
plt.show()

Warto przy tym wspomnieć, że na rysunkach jest podana wartość AUC - jest pole pod krzywą (area under curve) - kolejna miara jakości klasyfikatora. Możemy ją także otrzymać za pomocą funkcji roc_auc_score()

for clf in [lda, qda, nb]:
    print(clf, metrics.roc_auc_score(data.klasa, clf.predict_proba(data[["x","y"]])[:,1]))

## LinearDiscriminantAnalysis() 0.8300000000000001
## QuadraticDiscriminantAnalysis() 0.8433333333333333
## GaussianNB() 0.8366666666666666

Podział danych

W poprzednich punktach milcząco zakładaliśmy, że testujemy klasyfikatory na tych samych danych, na których zostały wyuczone, co jest oczywistym mankamentem - skutecznośc tych metod jest w efekcie przeszacowana. Poprawnie powinniśmy podzielić zbiór danych na część terningową i testową:

from sklearn.model_selection import train_test_split

x_train, x_test, y_train, y_test = train_test_split(data[["x","y"]], data.klasa, test_size=0.4)

Możemy teraz porównać jakoś klasyfikacji na tych zbiorach:

lda.fit(x_train, y_train)

LinearDiscriminantAnalysis()

print(metrics.classification_report(y_train, lda.predict(x_train)))

##               precision    recall  f1-score   support
## 
##            0       0.78      0.90      0.84        20
##            1       0.71      0.50      0.59        10
## 
##     accuracy                           0.77        30
##    macro avg       0.75      0.70      0.71        30
## weighted avg       0.76      0.77      0.75        30

print(metrics.classification_report(y_test, lda.predict(x_test)))

##               precision    recall  f1-score   support
## 
##            0       0.70      0.70      0.70        10
##            1       0.70      0.70      0.70        10
## 
##     accuracy                           0.70        20
##    macro avg       0.70      0.70      0.70        20
## weighted avg       0.70      0.70      0.70        20

Warto zwrócić uwagę, że rozkład klas w tych próbkach niekoniecznie jest równy:

print(y_train.groupby(y_train).count() / y_train.count())

## klasa
## 0    0.666667
## 1    0.333333
## Name: klasa, dtype: float64

print(y_test.groupby(y_test).count() / y_test.count())

## klasa
## 0    0.5
## 1    0.5
## Name: klasa, dtype: float64

aby sobie z poradzić, używamy stratyfikacji, czyli żądamy, aby rozkłąd był (w miarę) taki sam zarówno w próbce uczącej jak i testowej:

x_train, x_test, y_train, y_test = train_test_split(data[["x","y"]], data.klasa, test_size=0.4, stratify = data.klasa)


print(y_train.groupby(y_train).count() / y_train.count())

## klasa
## 0    0.6
## 1    0.4
## Name: klasa, dtype: float64

print(y_test.groupby(y_test).count() / y_test.count())

## klasa
## 0    0.6
## 1    0.4
## Name: klasa, dtype: float64

Kroswalidacja

W przypadku braku możliwości skorzystania z próby walidacyjnej (tzn. braku możliwości wydzielenia części danych) korzysta się z metody kroswalidacji (walidacji krzyżowej). W tym celu dzieli się dane na kilka (np. 5 lub 10) części, trenuje się klasyfikator na próbie powstałej przez usunięcie jednej części a testuje się na tej usuniętej. Procedurę powtarza się dla każdej częsci, zliczając pomyłki, a następnie sumując i dzieląc przez liczbę elementów, w ten sposób otrzymując oszacowanie błędu klasyfikatora.

Odpowiednią funkcją z pakietu sklearn jest np. KFold()

from sklearn.model_selection import KFold
kf = KFold(n_splits = 5)

for train, test in kf.split(data[["x","y"]]):
    print(test)

## [0 1 2 3 4 5 6 7 8 9]
## [10 11 12 13 14 15 16 17 18 19]
## [20 21 22 23 24 25 26 27 28 29]
## [30 31 32 33 34 35 36 37 38 39]
## [40 41 42 43 44 45 46 47 48 49]

Tu wypisaliśmy jedynie indeksy zbioru testowego, możemy też pokusić się o wizualizację:

fig, ax = plt.subplots(figsize=(10,5))

for i, (train, test) in enumerate(kf.split(data[["x","y"]])):
    ax.scatter(test, np.tile(i, len(test)), color = "blue", marker = "_", lw = 5)
    ax.scatter(train, np.tile(i, len(train)), color = "orange", marker = "_", lw = 5)

ax.set(yticks=np.arange(5),yticklabels=np.arange(1,6));

plt.show()

Inna opcją jest wykonanie kroswalidacji z przetasowaniem danych:

kf = KFold(n_splits = 5, shuffle = True)

fig, ax = plt.subplots(figsize=(10,5))

for i, (train, test) in enumerate(kf.split(data[["x","y"]])):
    ax.scatter(test, np.tile(i, len(test)), color = "blue", marker = "_", lw = 5)
    ax.scatter(train, np.tile(i, len(train)), color = "orange", marker = "_", lw = 5)

ax.set(yticks=np.arange(5),yticklabels=np.arange(1,6));

plt.show()

Wypada przy tym ocenić wyniki, wyznaczając srednią wartość błędów:

from sklearn.model_selection import cross_val_score

cvs = cross_val_score(lda, data[["x","y"]], data.klasa, cv = kf)
print(cvs)

## [0.9 0.7 0.5 0.6 0.7]

print(cvs.mean())

## 0.68

Jeszcze inną opcją jest wykorzystanie wielokrotenj kroswalidacji za pomocą RepeatedKFold() - tu z definicji mamy przetasowanie i możemy się przekonać, czy liczba podziałów wpływa na wyniki estymacji

from sklearn.model_selection import RepeatedKFold

def test(clf, data, m, n):
    kf = RepeatedKFold(n_repeats = m, n_splits = n)
    score = cross_val_score(clf, data[["x","y"]], data.klasa, cv = kf)
    return score.mean(), score.std()
    
z = np.array([test(lda, data, 5, i) for i in range(2, 40)])

plt.errorbar(range(2,40), z[:,0], yerr = z[:,1], fmt="o", capsize=4)

plt.show()

	priors priors: array-like of shape (n_classes,), default=None Prior probabilities of the classes. If specified, the priors are not adjusted according to the data.	None
	var_smoothing var_smoothing: float, default=1e-9 Portion of the largest variance of all features that is added to variances for calculation stability. .. versionadded:: 0.20	1e-09

	solver solver: {'svd', 'lsqr', 'eigen'}, default='svd' Solver to use, possible values: - 'svd': Singular value decomposition (default). Does not compute the covariance matrix, therefore this solver is recommended for data with a large number of features. - 'lsqr': Least squares solution. Can be combined with shrinkage or custom covariance estimator. - 'eigen': Eigenvalue decomposition. Can be combined with shrinkage or custom covariance estimator. .. versionchanged:: 1.2 `solver="svd"` now has experimental Array API support. See the :ref:`Array API User Guide ` for more details.	'svd'
	shrinkage shrinkage: 'auto' or float, default=None Shrinkage parameter, possible values: - None: no shrinkage (default). - 'auto': automatic shrinkage using the Ledoit-Wolf lemma. - float between 0 and 1: fixed shrinkage parameter. This should be left to None if `covariance_estimator` is used. Note that shrinkage works only with 'lsqr' and 'eigen' solvers. For a usage example, see :ref:`sphx_glr_auto_examples_classification_plot_lda.py`.	None
	priors priors: array-like of shape (n_classes,), default=None The class prior probabilities. By default, the class proportions are inferred from the training data.	None
	n_components n_components: int, default=None Number of components (<= min(n_classes - 1, n_features)) for dimensionality reduction. If None, will be set to min(n_classes - 1, n_features). This parameter only affects the `transform` method. For a usage example, see :ref:`sphx_glr_auto_examples_decomposition_plot_pca_vs_lda.py`.	None
	store_covariance store_covariance: bool, default=False If True, explicitly compute the weighted within-class covariance matrix when solver is 'svd'. The matrix is always computed and stored for the other solvers. .. versionadded:: 0.17	False
	tol tol: float, default=1.0e-4 Absolute threshold for a singular value of X to be considered significant, used to estimate the rank of X. Dimensions whose singular values are non-significant are discarded. Only used if solver is 'svd'. .. versionadded:: 0.17	0.0001
	covariance_estimator covariance_estimator: covariance estimator, default=None If not None, `covariance_estimator` is used to estimate the covariance matrices instead of relying on the empirical covariance estimator (with potential shrinkage). The object should have a fit method and a ``covariance_`` attribute like the estimators in :mod:`sklearn.covariance`. if None the shrinkage parameter drives the estimate. This should be left to None if `shrinkage` is used. Note that `covariance_estimator` works only with 'lsqr' and 'eigen' solvers. .. versionadded:: 0.24	None