Statistical Data Exploration LAB

Statistical Data Exploration

LAB 2

Metody LDA, QDA oraz naiwny Bayes

Głownym zdaniem zajęć będzie dokonanie dyskryminacji za pomocą metod LDA, QDA i NB, które zostały w przedstawione na Wykładzie 2 oraz Wykładzie 3. Zostaną do tego wykorzystane dedykowane funkcje z biblioteki scikit-learn, jednak aby w pełni zaznajomić się z tymi metodami warto będzie wyreślić najprostsze hiperpłaszczyzny dyskryminujące “na piechotę”. W tym celu należy poznać sposoby wyświetlania poziomic rozkładów 2D za pomocą funkcji contour(), a co za tym idzie również rozpinania siatki na zakresach współrzędnych.

Gęstość rozkładów 2D

Zaczynamy od uruchomienia dwóch modułów: numpy oraz pandas - oba już dobrze znamy. Podobnie jak podczas poprzedniego laboratorium generujemy punkty z dwóch rozkładów Gaussa o zadanych wektorach \(\mathbf{m}_1\) i \(\mathbf{m}_2\) oraz macierzach kowariancji \(\mathbf{S}_1\) i \(\mathbf{S}_1\):

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


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

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

Następnie wykreślamy punkty (potrzebujemy bibliotek seaborn oraz matplotlib.pyplot)

import seaborn as sns
import matplotlib.pyplot as plt

g = sns.relplot(data = data, x="x", y="y", hue="klasa")
g.ax.spines['top'].set_visible(True)
g.ax.spines['right'].set_visible(True)
g.ax.set_aspect('equal')

plt.axhline(0, color="black", linestyle=":")
plt.axvline(0, color="black", linestyle=":")

plt.show()

Tworzenie dwuwymiarowego rozkładu gęstości prawdopodobieństwa (PDF) jest oczywiście ideowo tożsame z przypadkiem jednowymiarowym - poniżej przypadek 1D, w którym korzystamy z pakietu scipy (funkcja norm do “obsługi” rozkładu - użyjemy norm.pdf do otrzymania gęstości prawdopodobieństwa) oraz seaborn do wykonania wykresu rozrzutu (rug):

from scipy.stats import norm

data1 = np.random.normal(0, 1, 100)
data2 = np.random.normal(-3, 0.5, 100)

g = sns.rugplot(data1)
sns.rugplot(data2)
x = np.arange(-8, 5, 0.05)
plt.plot(x, norm.pdf(x, 0, 1), color = "tab:blue")
plt.plot(x, norm.pdf(x, -3, 0.5), color = "tab:orange")

plt.show()

Różnica polega na tym, że w przypadku 2D musimy przekazać do funkcji (multivariate_normal) obie współprzędne (\(x_i\), \(y_i\)) każdej obserwacji \(\mathbf{x}_i\). Rzecz jasna, można po prostu zadeklarować macierz o odpowiednich rozmiarach wypełnioną zerami, a następnie korzystając z dwóch pętli wypełniać ją, jednak takie podejście jest mało “eleganckie”. Zamiast tego:

zdefiniujemy zakresy współprzędnych,
rozepniemy na nich siatkę (funkcja meshgrid() z numpy),
wyznaczymy wartości PDF 2D Gaussa dla każdego z tych punktów i zapiszemy je do macierzy.

Poniższy kod realizuje te procedury.

from scipy.stats import multivariate_normal as mvnorm

g = sns.relplot(data = data, x="x", y="y", hue="klasa")

x, y = np.mgrid[-5:5:.05, -5:5:.05]

contours = plt.contour(x, y, mvnorm.pdf(np.dstack([x, y]), m1, S1), 5, cmap=plt.cm.Blues, alpha = 0.5)
plt.clabel(contours, inline=True, fontsize=8)

contours1 = plt.contour(x, y, mvnorm.pdf(np.dstack([x, y]), m2, S2), 5, cmap=plt.cm.Oranges, alpha = 0.5)
plt.clabel(contours1, inline=True, fontsize=8)

g.ax.set_xlim(-5, 5);
g.ax.set_ylim(-5, 5);

plt.show()

Metoda LDA

Przechodzimy teraz do liniowej analizy dyskryminacji (LDA). Będziemy korzystali z pakietu sklearn, a w szczególności na początku z funkcji LinearDiscriminantAnalysis. Sama funkcja opowiada jedynie za proces konstrukcji pustego klasyfikatora, natomiast do nauczenia go na naszych danych korzystamy z funkcji fit

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
lda = LinearDiscriminantAnalysis()
lda.fit(data[["x","y"]], data.klasa)

LinearDiscriminantAnalysis()

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.

W tym momencie możemy już odwołac się do kolejnej funkcji - predict, która na bazie zgromadzonych danych jest w stanie dokonać przewidywania klas. Na tym etapie przekażemy te same dane, na których klasyfikator LDA został nauczony:

data_comp = pd.DataFrame({'org': data.klasa, 'pred': lda.predict(data[["x","y"]])})
data_comp

##     org  pred
## 0     0     0
## 1     0     0
## 2     0     0
## 3     0     0
## 4     0     0
## 5     0     0
## 6     0     0
## 7     0     0
## 8     0     0
## 9     0     0
## 10    0     0
## 11    0     0
## 12    0     0
## 13    0     0
## 14    0     0
## 15    0     0
## 16    0     0
## 17    0     0
## 18    0     0
## 19    0     0
## 20    0     0
## 21    0     0
## 22    0     0
## 23    0     0
## 24    0     0
## 25    0     0
## 26    0     0
## 27    0     0
## 28    0     0
## 29    0     0
## 30    0     0
## 31    0     0
## 32    0     0
## 33    0     0
## 34    0     0
## 35    0     0
## 36    0     0
## 37    0     0
## 38    0     0
## 39    0     0
## 40    1     1
## 41    1     1
## 42    1     1
## 43    1     1
## 44    1     1
## 45    1     1
## 46    1     1
## 47    1     1
## 48    1     1
## 49    1     1
## 50    1     0
## 51    1     1
## 52    1     1
## 53    1     0
## 54    1     1
## 55    1     1
## 56    1     1
## 57    1     1
## 58    1     1
## 59    1     0

Oczywiście, dużo wygodniej jest porównać dane oryginalne i przewidywane za pomocą tzw. macierzy kontyngencji (współwystępowania). W ML taki obiekt nazywa się macierzą pomyłek i szerzej będziemy o tym mówić na następnych zajęciach:

pd.crosstab(data_comp.org, data_comp.pred)

## pred   0   1
## org         
## 0     40   0
## 1      3  17

Oczywiście, decyzja nt. przynalezności do klasy jets podejmowana zgodnie z obliczoną przez klasyfikator wartością prawdopodobieństwa a posteriori, które uzyskamy za pomocą funkcji predict_proba. Wartości te dołożymy do naszej ramki danych (dla wygody zaokrąglimi je do 4 miejsca po przecinku)

data_comp[["p1","p2"]] = np.round(lda.predict_proba(data[["x","y"]]), 4)
data_comp

##     org  pred      p1      p2
## 0     0     0  0.9998  0.0002
## 1     0     0  1.0000  0.0000
## 2     0     0  0.9999  0.0001
## 3     0     0  1.0000  0.0000
## 4     0     0  0.9924  0.0076
## 5     0     0  0.9970  0.0030
## 6     0     0  1.0000  0.0000
## 7     0     0  0.9991  0.0009
## 8     0     0  0.9995  0.0005
## 9     0     0  1.0000  0.0000
## 10    0     0  1.0000  0.0000
## 11    0     0  1.0000  0.0000
## 12    0     0  0.9996  0.0004
## 13    0     0  0.9385  0.0615
## 14    0     0  0.9993  0.0007
## 15    0     0  0.9968  0.0032
## 16    0     0  0.9967  0.0033
## 17    0     0  0.9994  0.0006
## 18    0     0  0.9996  0.0004
## 19    0     0  0.8551  0.1449
## 20    0     0  0.9923  0.0077
## 21    0     0  1.0000  0.0000
## 22    0     0  0.9974  0.0026
## 23    0     0  1.0000  0.0000
## 24    0     0  1.0000  0.0000
## 25    0     0  0.9965  0.0035
## 26    0     0  0.9998  0.0002
## 27    0     0  0.9134  0.0866
## 28    0     0  0.9987  0.0013
## 29    0     0  1.0000  0.0000
## 30    0     0  0.9996  0.0004
## 31    0     0  0.9985  0.0015
## 32    0     0  1.0000  0.0000
## 33    0     0  0.9992  0.0008
## 34    0     0  0.9996  0.0004
## 35    0     0  1.0000  0.0000
## 36    0     0  0.9399  0.0601
## 37    0     0  0.9995  0.0005
## 38    0     0  0.9979  0.0021
## 39    0     0  0.9999  0.0001
## 40    1     1  0.0052  0.9948
## 41    1     1  0.0245  0.9755
## 42    1     1  0.2667  0.7333
## 43    1     1  0.1276  0.8724
## 44    1     1  0.0000  1.0000
## 45    1     1  0.0000  1.0000
## 46    1     1  0.0008  0.9992
## 47    1     1  0.0599  0.9401
## 48    1     1  0.0024  0.9976
## 49    1     1  0.0001  0.9999
## 50    1     0  0.9409  0.0591
## 51    1     1  0.0007  0.9993
## 52    1     1  0.0014  0.9986
## 53    1     0  0.5770  0.4230
## 54    1     1  0.0000  1.0000
## 55    1     1  0.0008  0.9992
## 56    1     1  0.0000  1.0000
## 57    1     1  0.0000  1.0000
## 58    1     1  0.0001  0.9999
## 59    1     0  0.5779  0.4221

Spróbujmy teraz odtworzyć wartości tych prawdopodobieństw; w przypadku dwóch klas de facto wystarczy nam informacja o \(p(1|x)\), które jest dane wzorem (reguła Bayesa):

\(p(1|x)=\frac{\pi_1p(x|1)}{p(x)} = \frac{π_1p(x|1)}{\pi_1p(x|1)+\pi_2p(x|2)}\)

Prawdopodobieństwa a priori \(\pi_1\) oraz \(\pi_2\) estymujemy jako liczności klas podzielone przez całkowitą liczbę obserwacji, natomiast gęstości prawdopodobieństwa \(p(x|1)\) i \(p(x|2)\) są dwuwymiarowymi rozkładami Gaussa o średnich \(\mathbf{\hat{m}}_1\) i \(\mathbf{\hat{m}}_2\) estymowanych z danych. Istotnym założeniem metody LDA jest to, że macierze kowariancji w obu klasach powinny być takie same. W praktyce uzyskuje się to poprzez uśrednienie wyestymowanych macierzy kowariancji \(\mathbf{\hat{S}}_1\) i \(\mathbf{\hat{S}}_2\)

\(\mathbf{\hat{S}}=\frac{1}{n−2}\sum\limits_{k=1,2}(n_k−1)\mathbf{\hat{S}}_k\)

Po tych wszystkich uwagach możemy teraz wyznaczyć wartości prawdopodobieństwa a posteriori przynależności do 1 klasy i porównać je z tymi, które otrzymaliśmy bezpośrednio z funkcji predict_proba()

def own_lda(X, pi1, pi2, m1, m2, S1, S2):
  return pi1 * mvnorm.pdf(X, m1, S1) / (pi1 * mvnorm.pdf(X, m1, S1) + pi2 * mvnorm.pdf(X, m2, S2))


mh = data.groupby("klasa").mean()
Sh = data.groupby("klasa").cov()
nn = n1 + n2

Shm = ((n1 - 1) * Sh.loc[0] + (n2 - 1) * Sh.loc[1]) / (nn - 2)

data_comp["pred_own"] = np.round(own_lda(data[["x", "y"]], n1/nn, n2/nn, mh.loc[0], mh.loc[1], Shm, Shm), 4)

data_comp

##     org  pred      p1      p2  pred_own
## 0     0     0  0.9998  0.0002    0.9998
## 1     0     0  1.0000  0.0000    1.0000
## 2     0     0  0.9999  0.0001    0.9999
## 3     0     0  1.0000  0.0000    1.0000
## 4     0     0  0.9924  0.0076    0.9924
## 5     0     0  0.9970  0.0030    0.9970
## 6     0     0  1.0000  0.0000    1.0000
## 7     0     0  0.9991  0.0009    0.9991
## 8     0     0  0.9995  0.0005    0.9995
## 9     0     0  1.0000  0.0000    1.0000
## 10    0     0  1.0000  0.0000    1.0000
## 11    0     0  1.0000  0.0000    1.0000
## 12    0     0  0.9996  0.0004    0.9996
## 13    0     0  0.9385  0.0615    0.9385
## 14    0     0  0.9993  0.0007    0.9993
## 15    0     0  0.9968  0.0032    0.9968
## 16    0     0  0.9967  0.0033    0.9967
## 17    0     0  0.9994  0.0006    0.9994
## 18    0     0  0.9996  0.0004    0.9996
## 19    0     0  0.8551  0.1449    0.8551
## 20    0     0  0.9923  0.0077    0.9923
## 21    0     0  1.0000  0.0000    1.0000
## 22    0     0  0.9974  0.0026    0.9974
## 23    0     0  1.0000  0.0000    1.0000
## 24    0     0  1.0000  0.0000    1.0000
## 25    0     0  0.9965  0.0035    0.9965
## 26    0     0  0.9998  0.0002    0.9998
## 27    0     0  0.9134  0.0866    0.9134
## 28    0     0  0.9987  0.0013    0.9987
## 29    0     0  1.0000  0.0000    1.0000
## 30    0     0  0.9996  0.0004    0.9996
## 31    0     0  0.9985  0.0015    0.9985
## 32    0     0  1.0000  0.0000    1.0000
## 33    0     0  0.9992  0.0008    0.9992
## 34    0     0  0.9996  0.0004    0.9996
## 35    0     0  1.0000  0.0000    1.0000
## 36    0     0  0.9399  0.0601    0.9399
## 37    0     0  0.9995  0.0005    0.9995
## 38    0     0  0.9979  0.0021    0.9979
## 39    0     0  0.9999  0.0001    0.9999
## 40    1     1  0.0052  0.9948    0.0052
## 41    1     1  0.0245  0.9755    0.0245
## 42    1     1  0.2667  0.7333    0.2667
## 43    1     1  0.1276  0.8724    0.1276
## 44    1     1  0.0000  1.0000    0.0000
## 45    1     1  0.0000  1.0000    0.0000
## 46    1     1  0.0008  0.9992    0.0008
## 47    1     1  0.0599  0.9401    0.0599
## 48    1     1  0.0024  0.9976    0.0024
## 49    1     1  0.0001  0.9999    0.0001
## 50    1     0  0.9409  0.0591    0.9409
## 51    1     1  0.0007  0.9993    0.0007
## 52    1     1  0.0014  0.9986    0.0014
## 53    1     0  0.5770  0.4230    0.5770
## 54    1     1  0.0000  1.0000    0.0000
## 55    1     1  0.0008  0.9992    0.0008
## 56    1     1  0.0000  1.0000    0.0000
## 57    1     1  0.0000  1.0000    0.0000
## 58    1     1  0.0001  0.9999    0.0001
## 59    1     0  0.5779  0.4221    0.5779

Suche liczby nigdy nie wytrzymują porównania z wykresem. Bibliotka sklearn jest wyposażona w funkcję DescisionBoundary, która, jak sama nazwa wskazuje, umożliwia wykreślenie granicy decyzyjnej, czyli innymi słowy hiperpłaszczyzny dyskryminacyjnej. Dodatkowo na rysunke naniesiemy orygianle dane:

from sklearn.inspection import DecisionBoundaryDisplay
disp = DecisionBoundaryDisplay.from_estimator(lda, data[["x","y"]], response_method="predict", 
                                              alpha=0.4, cmap=plt.cm.coolwarm)
plt.scatter(data.x, data.y, c=data["klasa"], edgecolor="k", cmap=plt.cm.coolwarm)

plt.show()

Zastanówy się teraz, w jaki sposób powstaje taki wykres. Otóż, jeśli rozepniamy na współrzędnych \(x\) i \(y\) siatkę (z pewną dokłądnością), to w każdym punkcie tej siatki bedziemy w stanie policzyć prawdopodobieństwo przynależności obserwacji do klasy 1 (lub 2, co w problemie binarnym wychodzi na jedno). Możemy wtedy otrzymać warstwice, podobne ideowo do tych, które otrzymaliśmy dla rozkładu Gaussa:

x, y = np.mgrid[-10:10:.05, -10:10:.05]
data_grid = pd.DataFrame({'x': x.ravel(), 'y': y.ravel()})

Z = lda.predict_proba(data_grid)[:,0].reshape(x.shape)

contours = plt.contour(x, y, Z, 10)
plt.clabel(contours, inline=True, fontsize=8)

plt.show()

Nas interesuje w tym przypadku warstwica 0.5 - po jej jednej stronie prawodpodobieństwo przynalezności do klasy, powiedzmy 0, jest mniejsze od 0.5, czyli obserwacja powinna zostac zaliczona do klasy 1, pod drugiej stronie - do klasy 0. Czyli właśnie ta warstwica jest granicą decyzyjną i jeśli nałożymy ją na poprzedni rysunek, to powinniśmy otrzymać pełną zgodność.

disp = DecisionBoundaryDisplay.from_estimator(lda, data[["x","y"]], response_method="predict", 
                                              alpha=0.4, cmap=plt.cm.coolwarm, grid_resolution=200)
                                              
xmin, xmax = disp.ax_.get_xlim()
ymin, ymax = disp.ax_.get_ylim()

plt.scatter(data.x, data.y, c=data["klasa"], edgecolor="k", cmap=plt.cm.coolwarm)
plt.contour(x, y, Z, [0.5],  linewidths=2.0, colors="black")

disp.ax_.set_xlim((xmin, xmax));
disp.ax_.set_ylim((ymin, ymax));

plt.show()

Na koniec wykorzystamy zaproponowaną przez nas funkcję own_lda - efekt powinien byc ten sam:

Zlda = own_lda(data_grid, n1/nn, n2/nn, mh.loc[0], mh.loc[1], Shm, Shm).reshape(x.shape)
disp = DecisionBoundaryDisplay.from_estimator(lda, data[["x","y"]], response_method="predict", 
                                              alpha=0.4, cmap=plt.cm.coolwarm, grid_resolution=200)
                                              
xmin, xmax = disp.ax_.get_xlim()
ymin, ymax = disp.ax_.get_ylim()

plt.scatter(data.x, data.y, c=data["klasa"], edgecolor="k", cmap=plt.cm.coolwarm)
plt.contour(x, y, Zlda, [0.5],  linewidths=2.0, colors="red")

disp.ax_.set_xlim((xmin, xmax));
disp.ax_.set_ylim((ymin, ymax));

plt.show()

Warto pamiętać, że to samo powinniśmy otrzymać prosto z teorii, przekształcając równanie prostej dyskryminującej

\(\ln \frac{\pi_1}{\pi_2} - \frac{1}{2} (\mathbf{m}_1-\mathbf{m}_2)^{T} \mathbf{S}^{-1} (\mathbf{m}_1+ \mathbf{m}_2) + (\mathbf{m}_1-\mathbf{m}_2)^{T} \mathbf{S}^{-1} \mathbf{x} = 0\)

do wygodniejszej postaci

\(a_x x + a_y y + b = 0\), gdzie
\(b = \ln \frac{\pi_1}{\pi_2} - \frac{1}{2}\mathbf{a}^{T}(\mathbf{m}_1 + \mathbf{m}_2)\)
\(\mathbf{a}^{T} = (\mathbf{m}_1-\mathbf{m}_2)^{T} \mathbf{S}^{-1}\)

co skutkuje

# Podstawienie do wzorów
a = (mh.loc[0] - mh.loc[1]) @ np.linalg.inv(Shm)
b = np.log(n1 / n2) - 0.5 * a @ (mh.loc[0] + mh.loc[1])
B = -b / a[1]
A = -a[0] / a[1]


# Granica decyzyjna klas
disp = DecisionBoundaryDisplay.from_estimator(lda, data[["x","y"]], response_method="predict", 
                                              alpha=0.4, cmap=plt.cm.coolwarm, grid_resolution=200)
                                              
xmin, xmax = disp.ax_.get_xlim()
ymin, ymax = disp.ax_.get_ylim()

# Punkty
plt.scatter(data.x, data.y, c=data["klasa"], edgecolor="k", cmap=plt.cm.coolwarm)

# Linia teoretyczna
plt.axline((0,B), slope = A, ls = "--", color="red")

disp.ax_.set_xlim((xmin, xmax));
disp.ax_.set_ylim((ymin, ymax));

plt.show()

potwierdzając, że wszystkie trzy (a w zasadzie cztery) podejścia są tożsame.

QDA

W przypadku kwadratowej analizy dyskryminacji (QDA) korzystamy z funkcji QuadraticDiscriminantAnalysis(), natomiast dalsze instrukcje są analogiczne do poprzednich.

from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis

qda = QuadraticDiscriminantAnalysis()
qda.fit(data[["x","y"]], data.klasa)

QuadraticDiscriminantAnalysis()

disp = DecisionBoundaryDisplay.from_estimator(qda, data[["x","y"]], response_method="predict", 
                                              alpha=0.4, cmap=plt.cm.coolwarm, grid_resolution=200)
                                        
plt.show()

Aby znów wykorzystać podejście bezpośredniego wyznaczania prawdopodobieństwa a posteriori, należy pamiętać, że w metodzie QDA zakłada się różne macierze kowariancji \(\mathbf{S}_1\) i \(\mathbf{S}_2\). Nasza funkcja own_lda() jest już na to przygotowana - wystarczy po prostu tym razem podstawić wyestymowane macierze:

Zqda = own_lda(data_grid, n1/nn, n2/nn, mh.loc[0], mh.loc[1], Sh.loc[0], Sh.loc[1]).reshape(x.shape)
disp = DecisionBoundaryDisplay.from_estimator(qda, data[["x","y"]], response_method="predict", 
                                              alpha=0.4, cmap=plt.cm.coolwarm, grid_resolution=200)
                                              
xmin, xmax = disp.ax_.get_xlim()
ymin, ymax = disp.ax_.get_ylim()

plt.scatter(data.x, data.y, c=data["klasa"], edgecolor="k", cmap=plt.cm.coolwarm)
plt.contour(x, y, Zqda, [0.5],  linewidths=2.0, colors="red")

disp.ax_.set_xlim((xmin, xmax));
disp.ax_.set_ylim((ymin, ymax));

plt.show()

Naiwny Bayes

Warto wspomnieć jaka jest (bardzo istotna!) różnica pomiędzy klasyfikatorem Bayesa, który stanowi podstawę metod LDA i QDA oraz naiwnym klasyfikatorem Bayesowskim. Otóż w tym ostatnim przypadku zakładamy niezależność poszczególnych składowych wektora \(\mathbf{x} = \left(x^1, x^2, ..., x^p \right)^{T}\). W efekcie prowadzi to do następującego wzoru na gęstości prawdopodobieństw \(p(\mathbf{x}|k)\):

\(p(\mathbf{x}|k) = \prod_{i=1}^{p}p(x^i|k)\)

Jeśli chodzo o kwestie programistyczne, to cała prodcedura przebiega podobnie, z tą różnicą, że musimy zaimportować funkcję GaussianNB() z pakietu sklearn.naive_bayes. Reszta wygląda tak samo:

from sklearn.naive_bayes import GaussianNB

nb = GaussianNB()
nb.fit(data[["x","y"]], data.klasa)

GaussianNB()

disp = DecisionBoundaryDisplay.from_estimator(nb, data[["x","y"]], response_method="predict", 
                                              alpha=0.4, cmap=plt.cm.coolwarm, grid_resolution=200)
                                              
plt.scatter(data.x, data.y, c=data["klasa"], edgecolor="k", cmap=plt.cm.coolwarm)
                                    
plt.show()

Wiele klas

Nic nie stoi na przeszkodzie, aby powyższe rozważania zastosować do przykładu wielu klas. Poniżej przykładowy kod dla 4 klas. Najpierw modyfikujemy trochę funkcję do losowania, aby była bardziej elastyczna:

def draw_datan(n, m, S):

    g = len(n)
    d = len(m[0])
    data = np.random.multivariate_normal(m[0], S[0], n[0])
    
    for i in range(1, g):
        data = np.concatenate([data, np.random.multivariate_normal(m[i], S[i], n[i])])
                     
    data = pd.DataFrame(data, columns = ['x'+str(i) for i in range(1, d + 1)])

    data["klasa"] = np.repeat(range(g), n)
    
    return data


m = [[-3, 1], [2, 2], [0, 0], [3, -2]]
S = [[[2, 0], [0, 2]], [[4, 0], [0, 4]], [[1, 0], [0, 1]], [[1, -0.5], [-0.5, 0.5]]]
n = 40, 20, 10, 30

data = draw_datan(n, m, S)

a następnie obok siebie rysujemy granice dla metod LDA, QDA i NB:

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

X = data.drop(columns=["klasa"])
lda.fit(X, data.klasa)

LinearDiscriminantAnalysis()

qda.fit(X, data.klasa)

QuadraticDiscriminantAnalysis()

nb.fit(X, data.klasa)

GaussianNB()

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize = (12,3))
DecisionBoundaryDisplay.from_estimator(lda, X, response_method="predict", 
                                              alpha=0.4, cmap=plt.cm.coolwarm, ax = ax1);
ax1.scatter(X.iloc[:,0], X.iloc[:,1], c=data["klasa"], edgecolor="k", cmap=plt.cm.coolwarm)

DecisionBoundaryDisplay.from_estimator(qda, X, response_method="predict", 
                                              alpha=0.4, cmap=plt.cm.coolwarm, ax = ax2);
ax2.scatter(X.iloc[:,0], X.iloc[:,1], c=data["klasa"], edgecolor="k", cmap=plt.cm.coolwarm)

DecisionBoundaryDisplay.from_estimator(nb, X, response_method="predict", 
                                              alpha=0.4, cmap=plt.cm.coolwarm, ax = ax3);
ax3.scatter(X.iloc[:,0], X.iloc[:,1], c=data["klasa"], edgecolor="k", cmap=plt.cm.coolwarm)

plt.show()

	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

	solver solver: {'svd', '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. - 'eigen': Eigenvalue decomposition. Can be combined with shrinkage or custom covariance estimator.	'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. Enabling shrinkage is expected to improve the model when some classes have a relatively small number of training data points compared to the number of features by mitigating overfitting during the covariance estimation step. This should be left to `None` if `covariance_estimator` is used. Note that shrinkage works only with 'eigen' solver.	None
	priors priors: array-like of shape (n_classes,), default=None Class priors. By default, the class proportions are inferred from the training data.	None
	reg_param reg_param: float, default=0.0 Regularizes the per-class covariance estimates by transforming S2 as ``S2 = (1 - reg_param) * S2 + reg_param * np.eye(n_features)``, where S2 corresponds to the `scaling_` attribute of a given class.	0.0
	store_covariance store_covariance: bool, default=False If True, the class covariance matrices are explicitly computed and stored in the `self.covariance_` attribute. .. versionadded:: 0.17	False
	tol tol: float, default=1.0e-4 Absolute threshold for the covariance matrix to be considered rank deficient after applying some regularization (see `reg_param`) to each `Sk` where `Sk` represents covariance matrix for k-th class. This parameter does not affect the predictions. It controls when a warning is raised if the covariance matrix is not full rank. .. 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 the 'eigen' solver.	None

	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

	solver solver: {'svd', '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. - 'eigen': Eigenvalue decomposition. Can be combined with shrinkage or custom covariance estimator.	'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. Enabling shrinkage is expected to improve the model when some classes have a relatively small number of training data points compared to the number of features by mitigating overfitting during the covariance estimation step. This should be left to `None` if `covariance_estimator` is used. Note that shrinkage works only with 'eigen' solver.	None
	priors priors: array-like of shape (n_classes,), default=None Class priors. By default, the class proportions are inferred from the training data.	None
	reg_param reg_param: float, default=0.0 Regularizes the per-class covariance estimates by transforming S2 as ``S2 = (1 - reg_param) * S2 + reg_param * np.eye(n_features)``, where S2 corresponds to the `scaling_` attribute of a given class.	0.0
	store_covariance store_covariance: bool, default=False If True, the class covariance matrices are explicitly computed and stored in the `self.covariance_` attribute. .. versionadded:: 0.17	False
	tol tol: float, default=1.0e-4 Absolute threshold for the covariance matrix to be considered rank deficient after applying some regularization (see `reg_param`) to each `Sk` where `Sk` represents covariance matrix for k-th class. This parameter does not affect the predictions. It controls when a warning is raised if the covariance matrix is not full rank. .. 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 the 'eigen' solver.	None