Teoria związana z tym zdaniem jest dostępna m.in. na stronie wykładu Statistical Data Exploration (Statystyczna Eksploracja Danych): Wykład 1 (dyskryminacja Fishera dla dwóch klas) oraz Wykład 2 (dyskryminacja Fishera dla wielu klas).
Pierwszym zadaniem, związanym z dyskryminacją Fishera jest znalezienie kierunku (wektora) a, który optymalnie rozdziela zrzutowane na niego obserwacje.
Dwie klasy
Załóżmy, że chcemy dokonać dyskryminacji zbiorów przedstawionych na poniższym rysunku:
Zaczniemy od uruchomienia bibliotek pandas oraz
numpy, które będa nam potrzebne, odpowiednio, do stworzenia
zbiorów danych oraz wyznaczenia statystyk.
Następnie tworzymy dwie ramki danych \(X1\) oraz \(X2\), w których umieszczone zostają współprzędne obserwacji należących do poszczególnych klas (odpowiednio “niebieskiej” i “pomarańczowej”). Kolejnym krokiem jest wyznaczenia wartości średnich \(\mathbf{m}_1\), \(\mathbf{m}_2\)
import pandas as pd
import numpy as np
X1 = pd.DataFrame({'x': [2,2,2,1,3], 'y': [2,1,3,2,2]})
X2 = pd.DataFrame({'x': [6,6,6,5,7], 'y': [0,1,-1,0,0]})
X = pd.concat([X1, X2], ignore_index=True)
X["class"] = np.repeat(["A","B"], 5)
m = X.groupby("class").mean()
m## x y
## class
## A 2.0 2.0
## B 6.0 0.0oraz macierzy kowariancji \(\mathbf{S}_1\), \(\mathbf{S}_2\).
S = X.groupby("class").cov()
S## x y
## class
## A x 0.5 0.0
## y 0.0 0.5
## B x 0.5 0.0
## y 0.0 0.5Zgodnie z Wykładem 1 wektor a otrzymujemy dokonując maksymalizacji wyrażenia
\(J = \frac{\left(\mathbf{a}^T\mathbf{m}_2 - \mathbf{a}^T\mathbf{m}_1\right)^2}{\mathbf{a}^T\mathbf{W}\mathbf{a}}\)
gdzie \(\mathbf{W}\) (macierz zmienności wewnątrzgrupowej) jest zdefiniowana jako
\(\mathbf{W} = \frac{1}{n - 2} \sum\limits_{k=1}^{2}(n_k - 1)\mathbf{S_k} = \frac{1}{n-2}\sum\limits_{k=1}^{2} \sum\limits_{i=1}^{n_k}(\mathbf{x}_{ki} - \mathbf{m}_k)(\mathbf{x}_{ki} - \mathbf{m}_k)^{T}\)
a \(n_1\) i \(n_2\) to liczby obserwacji w odpowiednich klasach, natomiast \(n\) jest ich całkowitą sumą (\(n=n_1+n_2\)). Maksymalizacja \(J\) prowadzi do następującego wyrażenia na wektor \(\mathbf{a}\)
\(\mathbf{a} \sim \mathbf{W}^{-1} (\mathbf{m}_2-\mathbf{m}_1)\)
Korzystając z tych informacji łatwo znajdujemy odpowiednie wartości,
używając funkcji invert() z modułu
numpy.linalg do odwrócenia macierzy oraz operatora
@ do wykonania mnożenia wektorowego jak również funkcji
count() z pakietu pandas, zliczającej rzędy w
ramce danych:
n = X.groupby("class")["x"].count()
nn = n.sum()
W = ((n["A"] - 1) * S.loc["A"] + (n["B"] - 1) * S.loc["B"])/(nn - 2)
a = np.linalg.inv(W) @ (m.loc["B"] - m.loc["A"])
a## array([ 8., -4.])Teraz możemy już wykreślić nasze punkty korzystając z funkcji
relplot() z pakietu seaborn, dla wygody
umieścić dodatkowe osie X i Y na górze i po prawej oraz wykonać
“krzyżak” związany z początkiem układu współrzędnych:
import matplotlib.pyplot as plt
import seaborn as sns
g = sns.relplot(data = X, x="x", y="y", s = 200, hue = "class")
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()
Rzecz jasna, sam wektor \(\mathbf{a}=(a_x,a_y)^T\) wskazuje jedynie
optymalny kierunek, natomiast nam zależy na wykreśleniu prostej. W tym
celu skorzystamy z funkcji axline(), której argumentami są
współczynniki np. punkt przez który pzrechodzi prosta oraz współczynnik
kierunkowy, w naszym przypadku wynosi on \(a_y/a_x\), natomiast punkt to \((0,0)\):
g = sns.relplot(data = X, x="x", y="y", s = 200, hue = "class")
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.axline((0,0), slope = a[1]/a[0])
g.ax.set_xlim(-1, 8);
g.ax.set_ylim(-3, 5);
plt.show()
Pozostaje kwestia wyznaczenia płaszczyzny dyskryminującej spełniającej zależność
\((\mathbf{m}_2-\mathbf{m}_1)^{T} \mathbf{W}^{-1} \left[\mathbf{x} - \frac{1}{2}(\mathbf{m}_1+\mathbf{m}_2) \right] = 0\)
którą w rozpatrywanym przypadku można przedstawić w trochę “łatwiejszy” sposób jako
\(a_x x + a_y y + b = 0\), gdzie \(b = -\frac{1}{2}\mathbf{a}^{T}(\mathbf{m}_1 + \mathbf{m}_2)\)
Finalnie
g = sns.relplot(data = X, x="x", y="y", s = 200, hue = "class")
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=":")
b = -0.5 * a @ (m.loc["A"] + m.loc["B"])
plt.axline((0,0), slope = a[1]/a[0], ls = "--", color="red")
plt.axline((0,-b/a[1]), slope = -a[0]/a[1], ls = "-", color="red")
g.ax.set_xlim(-1, 8);
g.ax.set_ylim(-3, 5);
plt.show()
Dwie klasy, rozkład Gaussa
Korzystnie jest sprawdzić działanie dyskryminacji Fisherowskiej dla
większe liczby danych. W tym celu wykorzystamy funkcje
numpy.random.multivariate() do losowania wartości z
wielowymiarowego rozkładu Gaussa
\(p(\mathbf{x}) = \frac{1}{(2 \pi)^\frac{p}{2}|\mathbf{S}|^{\frac{1}{2}}} \exp \left[-\frac{1}{2} (\mathbf{x} - \mathbf{m})^{T} \mathbf{S}^{-1} (\mathbf{x} - \mathbf{m}) \right]\)
mu1 = np.array([2, 2])
mu2 = np.array([6, 0])
sd1 = np.array([[1, 0], [0, 1]])
n1, n2 = 30, 20
data1 = np.random.multivariate_normal(mu1, sd1, n1)
data2 = np.random.multivariate_normal(mu2, sd1, n2)
data = np.concatenate((data1, data2), axis = 0)
data = pd.DataFrame(data, columns=["x", "y"])
data["class"] = np.repeat(["A", "B"], [n1, n2])
data## x y class
## 0 2.799958 1.866527 A
## 1 2.984211 1.522112 A
## 2 1.739409 3.390852 A
## 3 1.552810 3.209152 A
## 4 1.450630 2.322322 A
## 5 2.596080 3.468396 A
## 6 3.530678 1.822342 A
## 7 3.556398 2.423483 A
## 8 1.639491 1.103416 A
## 9 2.776634 2.217972 A
## 10 1.691018 1.926216 A
## 11 0.954718 2.594532 A
## 12 1.714227 2.378221 A
## 13 2.340956 2.997295 A
## 14 3.197227 1.657765 A
## 15 2.852307 2.482576 A
## 16 2.497191 1.318704 A
## 17 2.783735 1.859940 A
## 18 1.373503 3.399012 A
## 19 0.763892 1.744744 A
## 20 2.549815 0.920486 A
## 21 1.202263 1.344912 A
## 22 0.598853 2.410519 A
## 23 2.530399 3.340088 A
## 24 1.256073 3.512527 A
## 25 1.780183 3.508131 A
## 26 1.076478 1.725421 A
## 27 2.290766 2.198084 A
## 28 0.107845 1.331910 A
## 29 1.216992 1.448033 A
## 30 6.167218 1.726474 B
## 31 6.668141 -1.290404 B
## 32 7.094656 0.184525 B
## 33 5.947336 -1.291892 B
## 34 7.006920 -0.768779 B
## 35 7.773712 0.994549 B
## 36 5.502288 -0.107843 B
## 37 5.676746 -0.047088 B
## 38 5.927647 -1.163406 B
## 39 6.240993 -0.013153 B
## 40 6.095486 1.090713 B
## 41 5.601793 0.325502 B
## 42 6.866705 1.478801 B
## 43 7.154981 0.565399 B
## 44 5.446326 -1.346438 B
## 45 5.224730 1.618120 B
## 46 5.208416 0.201986 B
## 47 4.674512 -0.028047 B
## 48 5.830912 -0.184659 B
## 49 7.007596 -0.499889 BEfekt powinien wyglądać podobnie jak na poniższym rysunku
Wiele klas
Metodę Fishera można uogólnić na wiele klas. Posłuży nam ona do analizy poniższego przykładu
Dla wielu klas (Wykład 2) wektor \(\mathbf{a}\) maksymalizuje wyrażenie
\(J = \frac{\mathbf{a}^{T} \mathbf{B} \mathbf{a}}{\mathbf{a}^{T} \mathbf{W} \mathbf{a}}\)
gdzie \(\mathbf{B}\) jest macierzą zmienności międzygrupowej, zdefiniowaną jako
\(\mathbf{B} =\frac{1}{g-1}\sum\limits_{k=1}^g n_k (\mathbf{m}_k-\mathbf{m})(\mathbf{m}_k-\mathbf{m})^{T}\)
przy czym \(\mathbf{m}\) jest średnią po wszystkich punktach (bez rozróżniania klas), \(g\) jest liczbą klas, natomiast \(\mathbf{W}\) zostaje uogólnione jako
\(\mathbf{W} =\frac{1}{n-g}\sum\limits_{k=1}^g(n_k - 1) \mathbf{S}_k\)
Wektor \(\mathbf{a}\) maksymalizujący \(J\) jest wektorem własnym macierzy \(\mathbf{W}^{-1}\mathbf{B}\) odpowiadającym największej wartości własnej.
Podobnie jak w przypadku dwóch klas, zaczynamy od przypisania danych, naniesienia punktów na wykres oraz wyznaczenia średnich, kowariancji oraz liczności klas.
from string import ascii_uppercase
neb = np.array([[-1, 0, 0, 0, 1], [0, 0, 1, -1, 0]])
m = np.array([[2, 2], [4, 4], [5, 6], [9, 8], [10, 10]])
gg = m.shape[0]
npp = neb.shape[1]
data = np.tile(neb, gg) + np.repeat(m, npp, axis = 0).T
data = pd.DataFrame(data.T, columns = ["x", "y"])
data["class"] = np.repeat([i for i in ascii_uppercase[0:gg]], npp)
mm = data.groupby("class").mean()
S = data.groupby("class").cov()
mmall = data.drop(columns = ["class"]).mean()Następnie należy wyznaczyć ogólną średnią \(\mathbf{m}\) oraz macierze \(\mathbf{B}\) i \(\mathbf{W}\)
mmall = data.drop(columns = ["class"]).mean()
n = data.groupby("class")["x"].count()
ntot = n.sum()
W = np.array([(n.loc[i]-1) * S.loc[i] for i in ascii_uppercase[0:gg]]).sum(axis = 0) / (ntot - gg)
mtab = np.array((mm-mmall))
B = np.array([n.iloc[i] * (mtab[[i],].T @ mtab[[i],]) for i in range(gg)]).sum(axis = 0)
B## array([[230., 210.],
## [210., 200.]])Wreszcie pozostaje wyznaczyć wartości i wektory własne maceirzy \(\mathbf{U} = \mathbf{W}^{-1}\mathbf{B}\), wybrać wektor własny odpowiadający największej wartości własnej oraz wykreślić kierunek.
U = np.linalg.inv(W) @ B
lam, u = np.linalg.eig(U)
# chcemy wybrać nawiększy
idx = lam.argsort()[::-1]
lam = lam[idx]
u = u[:,idx]
g = sns.relplot(data = data, x="x", y = "y", s = 200, hue="class")
g.ax.set_aspect('equal')
g.ax.spines['top'].set_visible(True)
g.ax.spines['right'].set_visible(True)
plt.axhline(0, color="black", linestyle=":")
plt.axvline(0, color="black", linestyle=":")
plt.axline((0,0), slope = u[1,0]/u[0,0], ls = "--", color="red")
plt.show()
Pozostaje wreszcie kwestia rzutowania punktów na wybrany kierunek. Można
tego dokonać, korzystając z następującej transformacji
\(\hat{x} = \frac{A y + x}{A^2 + 1}\) \(\hat{y} = A \hat{x}\)
gdzie \(A = a_y / a_x\) to współczynnik kierunkowy powyższej prostej
A = u[1,0] / u[0,0]
data["x_r"] = (A*data.y+data.x)/(A**2 + 1)
data["y_r"] = A * data.x_r
g = sns.relplot(data = data, x="x", y = "y", s = 100, hue="class")
g.ax.set_aspect('equal')
g.ax.spines['top'].set_visible(True)
g.ax.spines['right'].set_visible(True)
plt.axhline(0, color="black", linestyle=":")
plt.axvline(0, color="black", linestyle=":")
plt.axline((0,0), slope = u[1,0]/u[0,0], ls = "--", color="red")
sns.scatterplot(data = data, x="x_r", y = "y_r", s = 100, hue="class", legend = False)
plt.show()