umz21/wyk/03_Regresja_logistyczna.ipynb

Uczenie maszynowe — zastosowania

3. 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)

3.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 pl
import pandas
import ipywidgets as widgets

%matplotlib inline
%config InlineBackend.figure_format = 'svg'

from IPython.display import display, Math, Latex

# Przydatne funkcje

# Wyświetlanie macierzy w LaTeX-u
def LatexMatrix(matrix):
    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

# Hipoteza (wersja macierzowa)
def hMx(theta, X):
    return X * theta

# Wykres danych (wersja macierzowa)
def regdotsMx(X, y, xlabel, ylabel):    
    fig = pl.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)
    pl.ylim(y.min() - 1, y.max() + 1)
    pl.xlim(np.min(X[:, 1]) - 1, np.max(X[:, 1]) + 1)
    return fig

# Wykres krzywej regresji (wersja macierzowa)
def reglineMx(fig, fun, theta, X):
    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 '-')))

# Legenda wykresu
def legend(fig):
    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

# Wersja macierzowa funkcji kosztu
def JMx(theta,X,y):
    m = len(y)
    J = 1.0 / (2.0 * m) * ((X * theta - y).T * ( X * theta - y))
    return J.item()

# Wersja macierzowa gradientu funkcji kosztu
def dJMx(theta,X,y):
    return 1.0 / len(y) * (X.T * (X * theta - y)) 

# Implementacja algorytmu gradientu prostego za pomocą numpy i macierzy
def GDMx(fJ, fdJ, theta, X, y, alpha=0.1, eps=10**-3):
    current_cost = fJ(theta, X, y)
    logs = [[current_cost, theta]]
    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
        logs.append([current_cost, theta]) 
    return theta, logs

thetaStartMx = np.matrix([0, 0]).reshape(2, 1)

# Funkcja, która rysuje próg
def threshold(fig, theta):
    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)

XMx3 = np.matrix(np.concatenate((np.ones((m, 1)), Xn), axis=1)).reshape(m, n_plus_1)
yMx3 = np.matrix(data_iris_setosa.values[:, 1]).reshape(m, 1)

fig = regdotsMx(XMx3, yMx3, '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.

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)

XMx3 = np.matrix(np.concatenate((np.ones((m, 1)), Xn), axis=1)).reshape(m, n_plus_1)
yMx3 = np.matrix(data_iris_setosa.values[:, 1]).reshape(m, 1)

# Regresja liniowa
theta_e3, logs3 = GDMx(JMx, dJMx, thetaStartMx, XMx3, yMx3, alpha=0.03, eps=0.000001)

fig = regdotsMx(XMx3, yMx3, 'x', 'Iris setosa?')
reglineMx(fig, hMx, theta_e3, XMx3)
legend(fig)

A gdyby tak przyjąć, że klasyfikator zwraca $1$ dla $h(x) > 0.5$ i $0$ w przeciwnym przypadku?

fig = regdotsMx(XMx3, yMx3, 'x', 'Iris setosa?')
theta_e3, logs3 = GDMx(JMx, dJMx, thetaStartMx, XMx3, yMx3, alpha=0.03, eps=0.000001)
reglineMx(fig, hMx, theta_e3, XMx3)
threshold(fig, theta_e3)  # pomarańczowa linia oznacza granicę między klasą "1" a klasą "0" wyznaczoną przez próg "h(x) = 0.5"
legend(fig)

• 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}}$$

# Funkjca logistycza

def logistic(x):
    return 1.0 / (1.0 + np.exp(-x))

import matplotlib.pyplot as plt
def plot_logistic():
    x = np.linspace(-5,5,200)
    y = logistic(x)

    fig = plt.figure(figsize=(7,5))
    ax = fig.add_subplot(111)
    plt.ylim(-.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}}$$

# Funkcja regresji logistcznej
def h(theta, X):
    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$).

# Funkcja kosztu dla regresji logistycznej
def J(h, theta, X, y):
    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

# Gradient dla regresji logistycznej
def dJ(h, theta, X, y):
    return 1.0 / len(y) * (X.T * (h(theta, X) - y))

# Metoda gradientu prostego dla regresji logistycznej
def GD(h, fJ, fdJ, theta, X, y, alpha=0.01, eps=10**-3, maxSteps=10000):
    errorCurr = fJ(h, theta, X, y)
    errors = [[errorCurr, theta]]
    while True:
        # oblicz nowe theta
        theta = theta - alpha * fdJ(h, theta, X, y)
        # raportuj poziom błędu
        errorCurr, errorPrev = fJ(h, theta, X, y), errorCurr
        # kryteria stopu
        if abs(errorPrev - errorCurr) <= eps:
            break
        if len(errors) > maxSteps:
            break
        errors.append([errorCurr, theta]) 
    return theta, errors

# Uruchomienie metody gradientu prostego dla regresji logistycznej
thetaBest, errors = GD(h, J, dJ, thetaStartMx, XMx3, yMx3, 
                       alpha=0.1, eps=10**-7, maxSteps=1000)
print("error =", errors[-1][0])
print("theta =", thetaBest)

error = [[0.05755617]]
theta = [[ 5.02530461]
 [-1.99174803]]

# Funkcja regresji logistycznej (wersja skalarna)
def scalar_logistic_regression_function(theta, x):
    return 1.0/(1.0 + np.exp(-(theta.item(0) + theta.item(1) * x)))

# Rysowanie progu
def threshold_val(fig, x_thr):
    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))

# Wykres krzywej regresji logistycznej
def logistic_regline(fig, theta, X):
    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 = regdotsMx(XMx3, yMx3, xlabel='x', ylabel='Iris setosa?')
logistic_regline(fig, thetaBest, XMx3)
threshold_val(fig, 2.5)

Traktujemy wartość $h_\theta(x)$ jako prawdopodobieństwo, że cecha przyjmie wartość pozytywną:

$$ 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.

data_iris_setosa_multi = pandas.DataFrame()
data_iris_setosa_multi['dł. płatków'] = data_iris['pl']  # "pl" oznacza "petal length" (długość płatków)
data_iris_setosa_multi['szer. płatków'] = data_iris['pw']  # "pw" oznacza "petal width" (szerokość płatków)
data_iris_setosa_multi['dł. dz. k.'] = data_iris['sl']  # "sl" oznacza "sepal length" (długość działek kielicha)
data_iris_setosa_multi['szer. dz. k.'] = data_iris['sw']  # "sw" oznacza "sepal width" (szerokość działek kielicha)
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])

   dł. płatków  szer. płatków  dł. dz. k.  szer. dz. k.  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

%matplotlib inline

import matplotlib.pyplot as plt
import seaborn

# 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)

XMx4 = np.matrix(np.concatenate((np.ones((m, 1)), Xn), axis=1)).reshape(m, n_plus_1)
yMx4 = np.matrix(data_iris_setosa_multi.values[:, n]).reshape(m, 1)

print(XMx4[:6])
print(yMx4[: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 = XMx4[:100], XMx4[100:]
yTrain, yTest = yMx4[:100], yMx4[100:]

# Macierz parametrów początkowych
thetaTemp = np.ones(5).reshape(5,1)

thetaBest, errors = GD(h, J, dJ, thetaTemp, XTrain, yTrain, 
                       alpha=0.1, eps=10**-7, maxSteps=1000)
print("error =", errors[-1][0])
print("theta =", thetaBest)

error = [[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 = \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):
    prob = h(theta, X).item()
    return (1, prob) if prob > 0.5 else (0, prob)

print("theta =", thetaBest)
print("x0 =", XTest[0])
print("h(x0) =", h(thetaBest, XTest[0]).item())
print("c(x0) =", classifyBi(thetaBest, XTest[0]), "\n")

theta = [[ 1.11414027]
 [-2.89324615]
 [-0.66543637]
 [ 0.14887292]
 [ 2.13284493]]
x0 = [[1.  6.3 1.8 7.3 2.9]]
h(x0) = 1.6061436959824844e-05
c(x0) = (0, 1.6061436959824844e-05) 

Obliczmy teraz skuteczność modelu (więcej na ten temat na następnym wykładzie, poświęconym metodom ewaluacji).

acc = 0.0
for i, rest in enumerate(yTest):
    cls, prob = classifyBi(thetaBest, XTest[i])
    if i < 10:
        print(int(yTest[i].item()), "<=>", cls, "-- prob:", round(prob, 4))
    acc += cls == yTest[i].item()

print("\nAccuracy:", acc / len(XTest))

0 <=> 0 -- prob: 0.0
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

3.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.
    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

# one-hot matrix
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
thetaTemp = np.ones(5).reshape(5,1)

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] $$

# Zapis macierzowy funkcji softmax
def softmax(X):
    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

# Dla każdej klasy wytrenujmy osobny klasyfikator dwuklasowy.

def trainMaxEnt(X, Y):
    n = X.shape[1]
    thetas = []
    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.
        thetaBest, errors = GD(h, J, dJ, theta, 
                               X, YBi, alpha=0.1, eps=10**-4)
        thetas.append(thetaBest)
    return thetas

# Macierze theta dla każdej klasy
thetas = trainMaxEnt(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.49207837]
 [ 0.32297899]
 [ 1.27518552]
 [-2.10699503]
 [-0.55887263]] 

Otrzymana macierz parametrów theta dla klasy 1:
 [[ 0.49531628]
 [ 0.36093408]
 [-1.27157499]
 [ 0.28948428]
 [-0.38775578]] 

Otrzymana macierz parametrów theta dla klasy 2:
 [[-0.14368885]
 [-1.59697421]
 [-2.03887846]
 [ 2.40910042]
 [ 2.54991032]] 

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(5):
    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.73217642  0.56835819  2.05282303]
Otrzymane prawdopodobieństwa:  [0.    0.185 0.815]
Wybrana klasa:  2
Obliczone y = 2
Oczekiwane y = 2

Dla x = [[1.  4.8 3.  1.4 0.3]]:
Po zastosowaniu regresji:  [ 2.75047926 -1.29797382 -9.78808679]
Otrzymane prawdopodobieństwa:  [0.983 0.017 0.   ]
Wybrana klasa:  0
Obliczone y = 0
Oczekiwane y = 0

Dla x = [[1.  7.1 3.  5.9 2.1]]:
Po zastosowaniu regresji:  [-6.99411744  0.13689343  1.96966296]
Otrzymane prawdopodobieństwa:  [0.    0.138 0.862]
Wybrana klasa:  2
Obliczone y = 2
Oczekiwane y = 2

Dla x = [[1.  5.9 3.  5.1 1.8]]:
Po zastosowaniu regresji:  [-5.52843441 -0.41148816  1.19377859]
Otrzymane prawdopodobieństwa:  [0.001 0.167 0.832]
Wybrana klasa:  2
Obliczone y = 2
Oczekiwane y = 2

Dla x = [[1.  6.1 2.6 5.6 1.4]]:
Po zastosowaniu regresji:  [-6.80386129  0.4691731   1.87452121]
Otrzymane prawdopodobieństwa:  [0.    0.197 0.803]
Wybrana klasa:  2
Obliczone y = 2
Oczekiwane y = 2