270 KiB
Uczenie maszynowe
4. Regresja logistyczna
Uwaga: Wbrew nazwie, _regresja logistyczna jest algorytmem służącym do rozwiązywania problemów klasyfikacji (wcale nie problemów regresji!)
Do demonstracji metody regresji ligistycznej wykorzystamy klasyczny zbiór danych _Iris flower data set, składający się ze 150 przykładów wartości 4 cech dla 3 gatunków irysów (kosaćców).
_Iris flower data set
- 150 przykładów
- 4 cechy
- 3 kategorie
![]() |
![]() |
![]() |
---|---|---|
_Iris setosa | _Iris virginica | _Iris versicolor |
kosaciec szczecinkowy | kosaciec amerykański | kosaciec różnobarwny |
4 cechy:
- długość działek kielicha (_sepal length,
sl
) - szerokość działek kielicha (_sepal width,
sw
) - długość płatka (_petal length,
pl
) - szerokość płatka (_petal width,
pw
)
4.1. Dwuklasowa regresja logistyczna
Zacznijmy od najprostszego przypadku:
- ograniczmy się do 2 klas
- ograniczmy się do 1 zmiennej
→ dwuklasowa regresja logistyczna jednej zmiennej
# Przydatne importy
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import pandas
import ipywidgets as widgets
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
from IPython.display import display, Math, Latex
# Przydatne funkcje
def LatexMatrix(matrix):
"""Wyświetlanie macierzy w LaTeX-u"""
ltx = r'\left[\begin{array}'
m, n = matrix.shape
ltx += '{' + ("r" * n) + '}'
for i in range(m):
ltx += r" & ".join([('%.4f' % j.item()) for j in matrix[i]]) + r" \\\\ "
ltx += r'\end{array}\right]'
return ltx
def h(theta, X):
"""Hipoteza (wersja macierzowa)"""
return X * theta
def regdots(X, y, xlabel, ylabel):
"""Wykres danych (wersja macierzowa)"""
fig = plt.figure(figsize=(16*.6, 9*.6))
ax = fig.add_subplot(111)
fig.subplots_adjust(left=0.1, right=0.9, bottom=0.1, top=0.9)
ax.scatter([X[:, 1]], [y], c='r', s=50, label='Dane')
ax.set_xlabel(xlabel)
ax.set_ylabel(ylabel)
ax.margins(.05, .05)
plt.ylim(y.min() - 1, y.max() + 1)
plt.xlim(np.min(X[:, 1]) - 1, np.max(X[:, 1]) + 1)
return fig
def regline(fig, fun, theta, X):
"""Wykres krzywej regresji (wersja macierzowa)"""
ax = fig.axes[0]
x0 = np.min(X[:, 1]) - 1.0
x1 = np.max(X[:, 1]) + 1.0
L = [x0, x1]
LX = np.matrix([1, x0, 1, x1]).reshape(2, 2)
ax.plot(L, fun(theta, LX), linewidth='2',
label=(r'$y={theta0:.2}{op}{theta1:.2}x$'.format(
theta0=float(theta[0][0]),
theta1=(float(theta[1][0]) if theta[1][0] >= 0 else float(-theta[1][0])),
op='+' if theta[1][0] >= 0 else '-')))
def legend(fig):
"""Legenda wykresu"""
ax = fig.axes[0]
handles, labels = ax.get_legend_handles_labels()
# try-except block is a fix for a bug in Poly3DCollection
try:
fig.legend(handles, labels, fontsize='15', loc='lower right')
except AttributeError:
pass
def J(theta,X,y):
"""Wersja macierzowa funkcji kosztu"""
m = len(y)
J = 1.0 / (2.0 * m) * ((X * theta - y).T * ( X * theta - y))
return J.item()
def dJ(theta,X,y):
"""Wersja macierzowa gradientu funkcji kosztu"""
return 1.0 / len(y) * (X.T * (X * theta - y))
def GD(fJ, fdJ, theta, X, y, alpha=0.1, eps=10**-3):
"""Implementacja algorytmu gradientu prostego za pomocą numpy i macierzy"""
current_cost = fJ(theta, X, y)
while True:
theta = theta - alpha * fdJ(theta, X, y) # implementacja wzoru
current_cost, prev_cost = fJ(theta, X, y), current_cost
if current_cost > 10000:
break
if abs(prev_cost - current_cost) <= eps:
break
return theta
theta_start = np.matrix([0, 0]).reshape(2, 1)
def threshold(fig, theta):
"""Funkcja, która rysuje próg"""
x_thr = (0.5 - theta.item(0)) / theta.item(1)
ax = fig.axes[0]
ax.plot([x_thr, x_thr], [-1, 2],
color='orange', linestyle='dashed',
label=u'próg: $x={:.2F}$'.format(x_thr))
# Wczytanie pełnych (oryginalnych) danych
data_iris = pandas.read_csv("iris.csv")
print(data_iris[:6])
sl sw pl pw Gatunek 0 5.2 3.4 1.4 0.2 Iris-setosa 1 5.1 3.7 1.5 0.4 Iris-setosa 2 6.7 3.1 5.6 2.4 Iris-virginica 3 6.5 3.2 5.1 2.0 Iris-virginica 4 4.9 2.5 4.5 1.7 Iris-virginica 5 6.0 2.7 5.1 1.6 Iris-versicolor
# Ograniczenie danych do 2 klas i 1 cechy
data_iris_setosa = pandas.DataFrame()
data_iris_setosa["dł. płatka"] = data_iris["pl"] # "pl" oznacza "petal length"
data_iris_setosa["Iris setosa?"] = data_iris["Gatunek"].apply(
lambda x: 1 if x == "Iris-setosa" else 0
)
print(data_iris_setosa[:6])
dł. płatka Iris setosa? 0 1.4 1 1 1.5 1 2 5.6 0 3 5.1 0 4 4.5 0 5 5.1 0
import numpy as np
# Przygotowanie danych
m, n_plus_1 = data_iris_setosa.values.shape
n = n_plus_1 - 1
Xn = data_iris_setosa.values[:, 0:n].reshape(m, n)
X = np.matrix(np.concatenate((np.ones((m, 1)), Xn), axis=1)).reshape(m, n_plus_1)
y = np.matrix(data_iris_setosa.values[:, 1]).reshape(m, 1)
# Regresja liniowa
theta_lin = GD(J, dJ, theta_start, X, y, alpha=0.03, eps=0.000001)
fig = regdots(X, y, "x", "Iris setosa?")
legend(fig)
Próba zastosowania regresji liniowej do problemu klasyfikacji
Najpierw z ciekawości sprawdźmy, co otrzymalibyśmy, gdybyśmy zastosowali regresję liniową do problemu klasyfikacji.
fig = regdots(X, y, "x", "Iris setosa?")
regline(fig, h, theta_lin, X)
legend(fig)
/tmp/ipykernel_1480/766189158.py:53: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) theta0=float(theta[0][0]), /tmp/ipykernel_1480/766189158.py:54: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) theta1=(float(theta[1][0]) if theta[1][0] >= 0 else float(-theta[1][0])),
A gdyby tak przyjąć, że klasyfikator zwraca $1$ dla $h(x) > 0.5$ i $0$ w przeciwnym przypadku?
fig = regdots(X, y, "x", "Iris setosa?")
theta_lin = GD(J, dJ, theta_start, X, y, alpha=0.03, eps=0.000001)
regline(fig, h, theta_lin, X)
threshold(
fig, theta_lin
) # pomarańczowa linia oznacza granicę między klasą "1" a klasą "0" wyznaczoną przez próg "h(x) = 0.5"
legend(fig)
/tmp/ipykernel_1480/766189158.py:53: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) theta0=float(theta[0][0]), /tmp/ipykernel_1480/766189158.py:54: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) theta1=(float(theta[1][0]) if theta[1][0] >= 0 else float(-theta[1][0])),
- Krzywa regresji liniowej jest niezbyt dopasowana do danych klasyfikacyjnych.
- Zastosowanie progu $y = 0.5$ nie zawsze pomaga uzyskać sensowny rezultat.
- $h(x)$ może przyjmować wartości mniejsze od $0$ i większe od $1$ – jak interpretować takie wyniki?
Wniosek: w przypadku problemów klasyfikacyjnych regresja liniowa nie wydaje się najlepszym rozwiązaniem.
Wprowadźmy zatem pewne modyfikacje do naszego modelu.
Zdefiniujmy następującą funkcję, którą będziemy nazywać funkcją _logistyczną (albo sigmoidalną):
Funkcja logistyczna (sigmoidalna):
$$g(x) = \dfrac{1}{1+e^{-x}}$$
def logistic(x):
"""Funkcja logistyczna"""
return 1.0 / (1.0 + np.exp(-x))
import matplotlib.pyplot as plt
def plot_logistic():
"""Wykres funkcji logistycznej"""
x = np.linspace(-5, 5, 200)
y = logistic(x)
fig = plt.figure(figsize=(7, 5))
ax = fig.add_subplot(111)
plt.ylim(-0.1, 1.1)
ax.plot(x, y, linewidth="2")
Wykres funkcji logistycznej $g(x) = \dfrac{1}{1+e^{-x}}$:
plot_logistic()
Funkcja logistyczna przekształca zbiór liczb rzeczywistych $\mathbb{R}$ w przedział otwarty $(0, 1)$.
Funkcja regresji logistycznej dla pojedynczego przykładu o cechach wyrażonych wektorem $x$:
$$h_\theta(x) = g(\theta^T , x) = \dfrac{1}{1 + e^{-\theta^T x}}$$
Dla całej macierzy cech $X$:
$$h_\theta(X) = g(X , \theta) = \dfrac{1}{1 + e^{-X \theta}}$$
def h(theta, X):
"""Funkcja regresji logistcznej"""
return 1.0 / (1.0 + np.exp(-X * theta))
Funkcja kosztu dla regresji logistycznej:
$$J(\theta) = -\dfrac{1}{m} \left( \sum_{i=1}^{m} y^{(i)} \log h_\theta( x^{(i)} ) + \left( 1 - y^{(i)} \right) \log \left( 1 - h_\theta (x^{(i)}) \right) \right)$$
Gradient dla regresji logistycznej (wersja macierzowa):
$$\nabla J(\theta) = \frac{1}{|\vec y|} X^T \left( h_\theta(X) - \vec y \right)$$
(Jedyna różnica między gradientem dla regresji logistycznej a gradientem dla regresji liniowej to postać $h_\theta$).
def J(h, theta, X, y):
"""Funkcja kosztu dla regresji logistycznej"""
m = len(y)
h_val = h(theta, X)
s1 = np.multiply(y, np.log(h_val))
s2 = np.multiply((1 - y), np.log(1 - h_val))
return -np.sum(s1 + s2, axis=0) / m
def dJ(h, theta, X, y):
"""Gradient dla regresji logistycznej"""
return 1.0 / len(y) * (X.T * (h(theta, X) - y))
def GD(h, fJ, fdJ, theta, X, y, alpha=0.01, eps=10**-3, max_steps=10000):
"""Metoda gradientu prostego dla regresji logistycznej"""
curr_cost = fJ(h, theta, X, y)
history = [[curr_cost, theta]]
while True:
# oblicz nowe theta
theta = theta - alpha * fdJ(h, theta, X, y)
# raportuj poziom błędu
prev_cost = curr_cost
curr_cost = fJ(h, theta, X, y)
# kryteria stopu
if abs(prev_cost - curr_cost) <= eps:
break
if len(history) > max_steps:
break
history.append([curr_cost, theta])
return theta, history
# Uruchomienie metody gradientu prostego dla regresji logistycznej
theta_best, history = GD(
h, J, dJ, theta_start, X, y, alpha=0.1, eps=10**-7, max_steps=1000
)
print(f"Koszt: {history[-1][0]}")
print(f"theta = {theta_best}")
Koszt: [[0.05755617]] theta = [[ 5.02530461] [-1.99174803]]
def scalar_logistic_regression_function(theta, x):
"""Funkcja regresji logistycznej (wersja skalarna)"""
return 1.0 / (1.0 + np.exp(-(theta.item(0) + theta.item(1) * x)))
def threshold_val(fig, x_thr):
"""Rysowanie progu"""
ax = fig.axes[0]
ax.plot(
[x_thr, x_thr],
[-1, 2],
color="orange",
linestyle="dashed",
label="próg: $x={:.2F}$".format(x_thr),
)
def logistic_regline(fig, theta, X):
"""Wykres krzywej regresji logistycznej"""
ax = fig.axes[0]
x0 = np.min(X[:, 1]) - 1.0
x1 = np.max(X[:, 1]) + 1.0
Arg = np.arange(x0, x1, 0.1)
Val = scalar_logistic_regression_function(theta, Arg)
ax.plot(Arg, Val, linewidth="2")
fig = regdots(X, y, xlabel="x", ylabel="Iris setosa?")
logistic_regline(fig, theta_best, X)
threshold_val(fig, 2.5)
Traktujemy wartość $h_\theta(x)$ jako prawdopodobieństwo zdefiniowane w następujący sposób:
$$ h_\theta(x) = P(y = 1 , | , x; \theta) $$
Jeżeli $h_\theta(x) > 0.5$, to dla takiego $x$ będziemy przewidywać wartość $y = 1$. W przeciwnym wypadku uprzewidzimy $y = 0$.
Dlaczego możemy traktować wartość funkcji regresji logistycznej jako prawdopodobieństwo?
Można o tym poczytać w zewnętrznych źródłach, np. https://towardsdatascience.com/logit-of-logistic-regression-understanding-the-fundamentals-f384152a33d1
Dwuklasowa regresja logistyczna: więcej cech
Jak postąpić, jeżeli będziemy mieli więcej niż jedną cechę $x$?
Weźmy teraz wszystkie cechy występujące w zbiorze _Iris:
- długość płatków (
pl
, _petal length) - szerokość płatków (
pw
, _petal width) - długość działek kielicha (
sl
, _sepal length) - szerokość działek kielicha (
sw
, _sepal width)
data_iris_setosa_multi = pandas.DataFrame()
for feature in ["pl", "pw", "sl", "sw"]:
data_iris_setosa_multi[feature] = data_iris[feature]
data_iris_setosa_multi["Iris setosa?"] = data_iris["Gatunek"].apply(
lambda x: 1 if x == "Iris-setosa" else 0
)
print(data_iris_setosa_multi[:6])
pl pw sl sw Iris setosa? 0 1.4 0.2 5.2 3.4 1 1 1.5 0.4 5.1 3.7 1 2 5.6 2.4 6.7 3.1 0 3 5.1 2.0 6.5 3.2 0 4 4.5 1.7 4.9 2.5 0 5 5.1 1.6 6.0 2.7 0
# Przygotowanie danych
m, n_plus_1 = data_iris_setosa_multi.values.shape
n = n_plus_1 - 1
Xn = data_iris_setosa_multi.values[:, 0:n].reshape(m, n)
X = np.matrix(np.concatenate((np.ones((m, 1)), Xn), axis=1)).reshape(m, n_plus_1)
y = np.matrix(data_iris_setosa_multi.values[:, n]).reshape(m, 1)
print(X[:6])
print(y[:6])
[[1. 1.4 0.2 5.2 3.4] [1. 1.5 0.4 5.1 3.7] [1. 5.6 2.4 6.7 3.1] [1. 5.1 2. 6.5 3.2] [1. 4.5 1.7 4.9 2.5] [1. 5.1 1.6 6. 2.7]] [[1.] [1.] [0.] [0.] [0.] [0.]]
# Podział danych na zbiór trenujący i testowy
XTrain, XTest = X[:100], X[100:]
yTrain, yTest = y[:100], y[100:]
# Macierz parametrów początkowych
theta_start = np.ones(5).reshape(5, 1)
theta_best, history = GD(
h, J, dJ, theta_start, XTrain, yTrain, alpha=0.1, eps=10**-7, max_steps=1000
)
print(f"Koszt: {history[-1][0]}")
print(f"theta = {theta_best}")
Koszt: [[0.006797]] theta = [[ 1.11414027] [-2.89324615] [-0.66543637] [ 0.14887292] [ 2.13284493]]
Funkcja decyzyjna regresji logistycznej
Funkcja decyzyjna mówi o tym, kiedy nasz algorytm będzie przewidywał $y = 1$, a kiedy $y = 0$:
$$ c(x) := \left\{ \begin{array}{ll} 1, & \mbox{gdy } P(y=1 , | , x; \theta) > 0.5 \\ 0 & \mbox{w przeciwnym przypadku} \end{array}\right. $$
$$ P(y=1 ,| , x; \theta) = h_\theta(x) $$
def classifyBi(theta, X):
"""Funkcja decyzyjna regresji logistycznej"""
prob = h(theta, X).item()
return (1, prob) if prob > 0.5 else (0, prob)
print(f"theta = {theta_best}")
print(f"x0 = {XTest[0]}")
print(f"h(x0) = {h(theta_best, XTest[0]).item()}")
print(f"c(x0) = {classifyBi(theta_best, XTest[0])}")
theta = [[ 1.11414027] [-2.89324615] [-0.66543637] [ 0.14887292] [ 2.13284493]] x0 = [[1. 6.3 1.8 7.3 2.9]] h(x0) = 1.606143695982487e-05 c(x0) = (0, 1.606143695982487e-05)
Obliczmy teraz skuteczność modelu (więcej na ten temat na następnym wykładzie, poświęconym metodom ewaluacji).
correct = 0
for i, rest in enumerate(yTest):
cls, prob = classifyBi(theta_best, XTest[i])
if i < 10:
print(f"{yTest[i].item():1.0f} <=> {cls} -- prob: {prob:6.4f}")
correct += cls == yTest[i].item()
accuracy = correct / len(XTest)
print(f"\nAccuracy: {accuracy}")
0 <=> 0 -- prob: 0.0000 1 <=> 1 -- prob: 0.9816 0 <=> 0 -- prob: 0.0001 0 <=> 0 -- prob: 0.0005 0 <=> 0 -- prob: 0.0001 1 <=> 1 -- prob: 0.9936 0 <=> 0 -- prob: 0.0059 0 <=> 0 -- prob: 0.0992 0 <=> 0 -- prob: 0.0001 0 <=> 0 -- prob: 0.0001 Accuracy: 1.0
4.2. Wieloklasowa regresja logistyczna
Przykład: wszystkie cechy ze zbioru _Iris, wszystkie 3 klasy ze zbioru Iris.
import pandas
data_iris = pandas.read_csv("iris.csv")
data_iris[:6]
sl | sw | pl | pw | Gatunek | |
---|---|---|---|---|---|
0 | 5.2 | 3.4 | 1.4 | 0.2 | Iris-setosa |
1 | 5.1 | 3.7 | 1.5 | 0.4 | Iris-setosa |
2 | 6.7 | 3.1 | 5.6 | 2.4 | Iris-virginica |
3 | 6.5 | 3.2 | 5.1 | 2.0 | Iris-virginica |
4 | 4.9 | 2.5 | 4.5 | 1.7 | Iris-virginica |
5 | 6.0 | 2.7 | 5.1 | 1.6 | Iris-versicolor |
# Przygotowanie danych
import numpy as np
features = ["sl", "sw", "pl", "pw"]
m = len(data_iris)
X = np.matrix(data_iris[features])
X0 = np.ones(m).reshape(m, 1)
X = np.hstack((X0, X))
y = np.matrix(data_iris[["Gatunek"]]).reshape(m, 1)
print("X = ", X[:4])
print("y = ", y[:4])
X = [[1. 5.2 3.4 1.4 0.2] [1. 5.1 3.7 1.5 0.4] [1. 6.7 3.1 5.6 2.4] [1. 6.5 3.2 5.1 2. ]] y = [['Iris-setosa'] ['Iris-setosa'] ['Iris-virginica'] ['Iris-virginica']]
Zamieńmy etykiety tekstowe w tablicy $y$ na wektory jednostkowe (_one-hot vectors):
$$ \begin{array}{ccc} \mbox{"Iris-setosa"} & \mapsto & \left[ \begin{array}{ccc} 1 & 0 & 0 \\ \end{array} \right] \\ \mbox{"Iris-virginica"} & \mapsto & \left[ \begin{array}{ccc} 0 & 1 & 0 \\ \end{array} \right] \\ \mbox{"Iris-versicolor"} & \mapsto & \left[ \begin{array}{ccc} 0 & 0 & 1 \\ \end{array} \right] \\ \end{array} $$
Wówczas zamiast wektora $y$ otrzymamy macierz $Y$:
$$ y ; = ; \left[ \begin{array}{c} y^{(1)} \\ y^{(2)} \\ y^{(3)} \\ y^{(4)} \\ y^{(5)} \\ \vdots \\ \end{array} \right] ; = ; \left[ \begin{array}{c} \mbox{"Iris-setosa"} \\ \mbox{"Iris-setosa"} \\ \mbox{"Iris-virginica"} \\ \mbox{"Iris-versicolor"} \\ \mbox{"Iris-virginica"} \\ \vdots \\ \end{array} \right] \quad \mapsto \quad Y ; = ; \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 1 & 0 \\ \vdots & \vdots & \vdots \\ \end{array} \right] $$
def mapY(y, cls):
m = len(y)
yBi = np.matrix(np.zeros(m)).reshape(m, 1)
yBi[y == cls] = 1.0
return yBi
def indicatorMatrix(y):
classes = np.unique(y.tolist())
m = len(y)
k = len(classes)
Y = np.matrix(np.zeros((m, k)))
for i, cls in enumerate(classes):
Y[:, i] = mapY(y, cls)
return Y
# Macierz jednostkowa
Y = indicatorMatrix(y)
# Podział danych na zbiór trenujący i testowy
XTrain, XTest = X[:100], X[100:]
YTrain, YTest = Y[:100], Y[100:]
# Macierz parametrów początkowych - niech skłąda się z samych jedynek
theta_start = np.ones(5).reshape(5, 1)
Regresja logistyczna jest metodą rozwiązywania problemów klasyfikacji dwuklasowej.
Aby znaleźć rozwiązanie problemu klasyfikacji wieloklasowej metodą regresji logistycznej, trzeba przekształcić problem na zbiór problemów klasyfikacji dwuklasowej.
Alternatywnie, można użyć wielomianowej regresji logistycznej (zob. https://machinelearningmastery.com/multinomial-logistic-regression-with-python).
Od regresji logistycznej dwuklasowej do wieloklasowej
- Irysy są przydzielone do trzech klas: _Iris-setosa (0), Iris-versicolor (1), Iris-virginica (2).
- Wiemy, jak stworzyć klasyfikatory dwuklasowe typu _Iris-setosa vs. Nie-Iris-setosa (tzw. one-vs-all).
- Możemy stworzyć trzy klasyfikatory $h_{\theta_1}, h_{\theta_2}, h_{\theta_3}$ (otrzymując trzy zestawy parametrów $\theta$) i wybrać klasę o najwyższym prawdopodobieństwie.
Pomoże nam w tym funkcja _softmax, która jest uogólnieniem funkcji logistycznej na większą liczbę wymiarów.
Funkcja _softmax
Odpowiednikiem funkcji logistycznej dla wieloklasowej regresji logistycznej jest funkcja $\mathrm{softmax}$:
$$ \textrm{softmax} \colon \mathbb{R}^k \to [0,1]^k $$
$$ \textrm{softmax}(z_1,z_2,\dots,z_k) = \left( \dfrac{e^{z_1}}{\sum_{i=1}^{k}e^{z_i}}, \dfrac{e^{z_2}}{\sum_{i=1}^{k}e^{z_i}}, \ldots, \dfrac{e^{z_k}}{\sum_{i=1}^{k}e^{z_i}} \right) $$
$$ \textrm{softmax}( \left[ \begin{array}{c} \theta_1^T x \\ \theta_2^T x \\ \vdots \\ \theta_k^T x \end{array} \right] ) = \left[ \begin{array}{c} P(y=1 , | , x;\theta_1,\ldots,\theta_k) \\ P(y=2 , | , x;\theta_1,\ldots,\theta_k) \\ \vdots \\ P(y=k , | , x;\theta_1,\ldots,\theta_k) \end{array} \right] $$
def softmax(X):
"""Funkcja softmax (wersja macierzowa)"""
return np.exp(X) / np.sum(np.exp(X))
Wartości funkcji $\mathrm{softmax}$ sumują się do 1:
Z = np.matrix([[2.1, 0.5, 0.8, 0.9, 3.2]])
P = softmax(Z)
print(np.sum(P))
0.9999999999999999
def multiple_binary_classifiers(X, Y):
n = X.shape[1]
thetas = []
# Dla każdej klasy wytrenujmy osobny klasyfikator dwuklasowy.
for c in range(Y.shape[1]):
YBi = Y[:, c]
theta = np.matrix(np.random.random(n)).reshape(n, 1)
# Macierz parametrów theta obliczona dla każdej klasy osobno.
theta_best, history = GD(h, J, dJ, theta, X, YBi, alpha=0.1, eps=10**-4)
thetas.append(theta_best)
return thetas
# Macierze theta dla każdej klasy
thetas = multiple_binary_classifiers(XTrain, YTrain)
for c, theta in enumerate(thetas):
print(f"Otrzymana macierz parametrów theta dla klasy {c}:\n", theta, "\n")
Otrzymana macierz parametrów theta dla klasy 0: [[ 0.17186737] [ 0.28150146] [ 1.45560929] [-2.22394386] [-0.15036765]] Otrzymana macierz parametrów theta dla klasy 1: [[ 0.573234 ] [-0.17098277] [-0.71530812] [ 0.87669965] [-1.11447904]] Otrzymana macierz parametrów theta dla klasy 2: [[-0.32900921] [-1.66338352] [-1.78563311] [ 2.34323547] [ 2.66006068]]
Funkcja decyzyjna wieloklasowej regresji logistycznej
$$ c = \mathop{\textrm{arg},\textrm{max}}_{i \in \{1, \ldots ,k\}} P(y=i|x;\theta_1,\ldots,\theta_k) $$
def classify(thetas, X, debug=False):
regs = np.array([(X * theta).item() for theta in thetas])
if debug:
print("Po zastosowaniu regresji: ", regs)
probs = softmax(regs)
if debug:
print("Otrzymane prawdopodobieństwa: ", np.around(probs, decimals=3))
result = np.argmax(probs)
if debug:
print("Wybrana klasa: ", result)
return result
for i in range(4):
print(f"Dla x = {XTest[i]}:")
YPredicted = classify(thetas, XTest[i], debug=True)
print(f"Obliczone y = {YPredicted}")
print(f"Oczekiwane y = {np.argmax(YTest[i])}")
print()
Dla x = [[1. 7.3 2.9 6.3 1.8]]: Po zastosowaniu regresji: [-7.83341312 0.76781174 1.90044778] Otrzymane prawdopodobieństwa: [0. 0.244 0.756] Wybrana klasa: 2 Obliczone y = 2 Oczekiwane y = 2 Dla x = [[1. 4.8 3. 1.4 0.3]]: Po zastosowaniu regresji: [ 2.73127055 -1.50037187 -9.59160157] Otrzymane prawdopodobieństwa: [0.986 0.014 0. ] Wybrana klasa: 0 Obliczone y = 0 Oczekiwane y = 0 Dla x = [[1. 7.1 3. 5.9 2.1]]: Po zastosowaniu regresji: [-6.89968523 0.04545391 1.91528519] Otrzymane prawdopodobieństwa: [0. 0.134 0.866] Wybrana klasa: 2 Obliczone y = 2 Oczekiwane y = 2 Dla x = [[1. 5.9 3. 5.1 1.8]]: Po zastosowaniu regresji: [-5.4132216 -0.11638277 1.23873883] Otrzymane prawdopodobieństwa: [0.001 0.205 0.794] Wybrana klasa: 2 Obliczone y = 2 Oczekiwane y = 2