Autorzy¶
Studia 2 stopnia - Informatyka (zaoczne) 1 rok
+-
+
- Krzysztof Bojakowski +
- Adam Stelmaszyk +
- Patryk Osiński +
- Marcin Jakubik +
Zadanie¶
Klasyfikacja za pomocą naiwnej metody bayesowskiej (rozkłady ciągłe). Implementacja powinna założyć, że cechy są ciągłe (do wyboru rozkład normalny i jądrowe wygładzenie). Na wejściu oczekiwany jest zbiór, który zawiera p-cech ciągłych, wektor etykiet oraz wektor prawdopodobieństw a priori dla klas. Na wyjściu otrzymujemy prognozowane etykiety oraz prawdopodobieństwa a posteriori. Dodatkową wartością może być wizualizacja obszarów decyzyjnych w przypadku dwóch cech.
+Zbiór danych¶
Do prezentacji wyników użyliśmy danych pochodzących z UCI ML Breast Cancer Wisconsin (Diagnostic) dataset. Zbiór zawiera 30 cech, na podstawie których klasyfikowany jest nowotwór - jako złośliwy lub niezsłośliwy.
+import matplotlib.pyplot as plt
+from sklearn.datasets import load_breast_cancer
+import numpy as np
+import pandas as pd
+from sklearn.model_selection import train_test_split
+
+breast_cancer = load_breast_cancer()
+
+# Sprawdzenie poprawności rozmiarów
+assert len(breast_cancer.data[0]) == len(breast_cancer.feature_names)
+assert len(breast_cancer.target) == len(breast_cancer.data)
+
+values = np.c_[breast_cancer.data, breast_cancer.target]
+columns = np.array(list(breast_cancer.feature_names) + ['target'])
+
+df = pd.DataFrame(values, columns=columns)
+
+X = df.iloc[:, 0:-1]
+y = df.iloc[:,-1]
+
+df.head()
+
+ | mean radius | +mean texture | +mean perimeter | +mean area | +mean smoothness | +mean compactness | +mean concavity | +mean concave points | +mean symmetry | +mean fractal dimension | +... | +worst texture | +worst perimeter | +worst area | +worst smoothness | +worst compactness | +worst concavity | +worst concave points | +worst symmetry | +worst fractal dimension | +target | +
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | +17.99 | +10.38 | +122.80 | +1001.0 | +0.11840 | +0.27760 | +0.3001 | +0.14710 | +0.2419 | +0.07871 | +... | +17.33 | +184.60 | +2019.0 | +0.1622 | +0.6656 | +0.7119 | +0.2654 | +0.4601 | +0.11890 | +0.0 | +
1 | +20.57 | +17.77 | +132.90 | +1326.0 | +0.08474 | +0.07864 | +0.0869 | +0.07017 | +0.1812 | +0.05667 | +... | +23.41 | +158.80 | +1956.0 | +0.1238 | +0.1866 | +0.2416 | +0.1860 | +0.2750 | +0.08902 | +0.0 | +
2 | +19.69 | +21.25 | +130.00 | +1203.0 | +0.10960 | +0.15990 | +0.1974 | +0.12790 | +0.2069 | +0.05999 | +... | +25.53 | +152.50 | +1709.0 | +0.1444 | +0.4245 | +0.4504 | +0.2430 | +0.3613 | +0.08758 | +0.0 | +
3 | +11.42 | +20.38 | +77.58 | +386.1 | +0.14250 | +0.28390 | +0.2414 | +0.10520 | +0.2597 | +0.09744 | +... | +26.50 | +98.87 | +567.7 | +0.2098 | +0.8663 | +0.6869 | +0.2575 | +0.6638 | +0.17300 | +0.0 | +
4 | +20.29 | +14.34 | +135.10 | +1297.0 | +0.10030 | +0.13280 | +0.1980 | +0.10430 | +0.1809 | +0.05883 | +... | +16.67 | +152.20 | +1575.0 | +0.1374 | +0.2050 | +0.4000 | +0.1625 | +0.2364 | +0.07678 | +0.0 | +
5 rows × 31 columns
+Implementacja naiwnej metody Bayesowskiej¶
Implementacja jest w dwóch wariantach:
+-
+
- Rozkład normalny - Gaussian (Normal) Distribution +
- Rozkład log-normalny - Log-Gaussian (Log-Normal) Distribution +
Matematyczny opis algorytmu¶
Twierdzenie Bayesa¶
${X}$ będzie pewnym zdarzeniem o ${n}$ cechach $$ X = (x_1, x_2, \ldots, x_n)$$ +prawdopodobieństwo a posteriori ${P(C_k | X)}$ dla klasy ${C_k}$ jest dane przez: $$P(C_k | X) = \frac{P(X | C_k) P(C_k)}{P(X)}$$
+Gdzie: +$$ +\begin{aligned} +& \text{(1)}\hspace{1cm}P(C_k | X)\text{ to prawdopodobieństwo (probability) a posteriori klasy } C_k \text{ pod warunkiem zajścia zdarzenia } X \\ +& \text{(2)}\hspace{1cm}P(X | C_k) \text{ to prawdopodobieństwo (likelihood) obserwacji zdarzenia } X \text{ pod warunkiem klasy } C_k \\ +& \text{(3)}\hspace{1cm}P(C_k) \text{ to prawdopodobieństwo apriori klasy } C_k \\ +& \text{(4)}\hspace{1cm}P(X) \text{ to całkowite prawdopodobieństwo obserwacji zdarzenia X we wszystkich klasach} +\end{aligned} +$$
+Wzory na średnią i odchylenie standardowe¶
Średnia (Mean)¶
Średnia ${\mu}$ zestawu danych ${X}$ składającego się z ${n}$ obserwacji jest obliczana jako:
+$$ +\mu = \frac{1}{n} \sum_{i=1}^{n} x_i +$$
+Gdzie:
+-
+
- ${\mu}$ to średnia (mean) +
- ${n}$ to liczba obserwacji w zestawie danych +
- ${x_i}$ to pojedyncza obserwacja w zestawie danych +
Odchylenie standardowe (Standard Deviation)¶
Odchylenie standardowe ${ \sigma }$ zestawu danych ${ X }$ składającego się z ${ n }$ obserwacji jest obliczane jako pierwiastek kwadratowy z wariancji:
+$$ +\sigma = \sqrt{\frac{1}{n} \sum_{i=1}^{n} (x_i - \mu)^2} +$$
+Gdzie:
+-
+
- ${\sigma}$ to odchylenie standardowe (standard deviation) +
- ${\mu}$ to średnia (mean) zestawu danych +
- ${n}$ to liczba obserwacji w zestawie danych +
- ${x_i}$ to pojedyncza obserwacja w zestawie danych +
Prawdpodobiestwo (Likelihood) obserwacji zdarzenia - rozkład normalny¶
$$ +\text{Likelihood}(x | \mu, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x - \mu)^2}{2\sigma^2}\right) +$$
+$$ +\text{Posterior} = \text{Likelihood} \times \text{Prior} +$$
+Prawdopodobieństwo (Log-Likelihood) obserwacji zdarzenia - rozkład log-normalny¶
$$ +\text{Log-Likelihood}(x | \mu, \sigma) = -\frac{1}{2} \log(2\pi\sigma^2) - \frac{(x - \mu)^2}{2\sigma^2} +$$
+$$ +\text{Log-Posterior} = \text{Log-Likelihood} + \text{Log-Prior} +$$
+import numpy as np
+from scipy.stats import norm
+
+class GaussianNaiveBayes:
+ def __init__(self, priors=None, var_smoothing=1e-8, log=False):
+ self.priors = priors
+ self.classes = None
+ self.means = None
+ self.stds = None
+ self.log = log
+ self.var_smoothing = var_smoothing
+
+ def fit(self, X, y):
+ self.epsilon = self.var_smoothing * np.var(X, axis=0).max()
+
+ self.classes = np.unique(y)
+ self.means = np.zeros((len(self.classes), X.shape[1]))
+ self.stds = np.zeros((len(self.classes), X.shape[1]))
+
+ for idx, cls in enumerate(self.classes):
+ X_c = X[y == cls]
+ self.means[idx, :] = X_c.mean(axis=0)
+ self.stds[idx, :] = X_c.std(axis=0)
+
+ if self.priors is None:
+ self.priors = np.bincount(y) / len(y)
+
+ self.stds += self.epsilon
+
+ def predict(self, X):
+ posteriors = self.predict_proba(X)
+ return np.argmax(posteriors, axis=1)
+
+ def score(self, X, y_true):
+ posteriors = self.predict_proba(X)
+ y_pred = np.argmax(posteriors, axis=1)
+
+ correct_predictions = np.sum(y_true == y_pred)
+ accuracy = correct_predictions / len(y_true)
+ return accuracy
+
+ def score_visualize(self, X, y_true):
+ """
+ Funkcja rysuje dwa wykresy słupkowe przedstawiające rzeczywiste wartości i przewidywane wartości.
+ """
+ posteriors = self.predict_proba(X)
+ y_pred = np.argmax(posteriors, axis=1)
+
+ # Liczymy ilość wystąpień każdej klasy dla rzeczywistych i przewidywanych wartości
+ true_counts = np.bincount(y_true)
+ pred_counts = np.bincount(y_pred)
+
+ # Tworzymy subplots - dwa wykresy obok siebie
+ fig, axes = plt.subplots(1, 2, figsize=(12, 6))
+
+ # Pierwszy wykres - rzeczywiste wartości
+ axes[0].bar([0, 1], true_counts, color=['blue', 'orange'])
+ axes[0].set_title('Rzeczywiste wartości')
+ axes[0].set_xticks([0, 1])
+ axes[0].set_xticklabels(['Nowotwór złośliwy', 'Nowotwór łagodny'])
+ axes[0].set_ylabel('Liczba wystąpień')
+
+ # Drugi wykres - przewidywane wartości
+ axes[1].bar([0, 1], pred_counts, color=['blue', 'orange'])
+ axes[1].set_title('Przewidywane wartości')
+ axes[1].set_xticks([0, 1])
+ axes[1].set_xticklabels(['Nowotwór złośliwy', 'Nowotwór łagodny'])
+ axes[1].set_ylabel('Liczba wystąpień')
+
+ # Wyświetlamy wykresy
+ plt.show()
+
+ def predict_proba(self, X):
+ if self.log:
+ return self._calculate_posterior_log(X)
+ else:
+ return self._calculate_posterior_normal_vectorized(X)
+
+ def _normal_gaussian_pdf(self, x, mean, std):
+ """
+ Implementacja Normal Gaussian (likelihood) z użyciem SciKit - używana do testów
+ """
+ return np.prod(norm.pdf(x, mean, std))
+
+ def normal_gaussian_pdf(self, x, mean, std):
+ """
+ implementacja Normal Gaussian (likelihood)
+
+ kształt x jest zmieniany na (n_samples, 1, n_features) aby numpy mógł broadcastować
+ operacje z mediąną, której kształt to (n_classes, n_features)
+ """
+ coefficient = 1 / (np.sqrt(2 * np.pi * std**2))
+ exponent = np.exp(-(np.square(x[:, np.newaxis] - mean)) / (2 * np.square(std)))
+ return coefficient * exponent
+
+ def log_gaussian_pdf(self, x, mean, std):
+ """
+ implementacja Log-Normal Gaussian (likelihood)
+ """
+ return -0.5 * np.log(2 * np.pi * std**2) - ((x - mean)**2 / (2 * std**2))
+
+ def _calculate_posterior_log(self, x):
+ """
+ Log-normal
+ Użyliśmy log-sum-exp trick dla numerycznej stabilności.
+ """
+ log_likelihoods = np.zeros((x.shape[0], len(self.classes)))
+ for idx, _ in enumerate(self.classes):
+ prior = np.log(self.priors[idx])
+ log_likelihood = np.sum(self.log_gaussian_pdf(x, self.means[idx, :], self.stds[idx, :]), axis=1)
+ log_likelihoods[:, idx] = log_likelihood + prior
+
+ # The Log-Sum-Exp Trick
+ max_log_likelihoods = np.max(log_likelihoods, axis=1, keepdims=True)
+ log_posteriors = log_likelihoods - max_log_likelihoods
+ posteriors = np.exp(log_posteriors)
+ posteriors = posteriors / np.sum(posteriors, axis=1, keepdims=True)
+ return posteriors
+
+ def _calculate_posterior_normal_vectorized(self, X):
+ """
+ Implementacja z użyciem wektorów - zwraca takie same wyniki jak _calculate_posterior_normal_iterative
+ Działa szybciej i pomija konieczność iteracji w pętli.
+ """
+ likelihoods = self.normal_gaussian_pdf(X, self.means, self.stds)
+ posteriors = likelihoods.prod(axis=2) * self.priors
+ posteriors = posteriors / np.sum(posteriors, axis=1, keepdims=True)
+ return posteriors
+
+ def _calculate_posterior_normal_iterative(self, X):
+ """
+ Implementacja iteracyjna - zwraca takie same wyniki _calculate_posterior_normal_vectorized
+ Używana do testów, aby sprawdzić czy wersja z użyciem wektorów działa poprawnie
+ """
+ posteriors = []
+ for x in X:
+ likelihoods = np.zeros((len(self.classes)))
+ for idx, _ in enumerate(self.classes):
+ likelihoods[idx] = self._normal_gaussian_pdf(x, self.means[idx, :], self.stds[idx, ]) * self.priors[idx]
+ posteriors.append(likelihoods)
+ posteriors = np.array(posteriors)
+ posteriors = posteriors / np.sum(posteriors, axis=1, keepdims=True)
+ return posteriors
+
+def pp_predictions(preds, probs):
+ i = 1
+ for pred, prob in zip(preds, probs):
+ print(f"\n=== {i} ===")
+
+ if pred == 0:
+ print(f"Predykcja: {pred} - (Nowotwór złośliwy)")
+ else:
+ print(f"Predykcja: {pred} - (Nowotwór łagodny)")
+
+ print(f"Prawdpodobieństwo (Złośliwy, Łagodny): {np.round(prob, 4)}")
+ i += 1
+
Prezentacja działania algorytmu¶
+# Podział zbioru danych na treningowy i testowy
+X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.1, random_state=48)
+
+gnb = GaussianNaiveBayes(log=False)
+gnb.fit(X_train.values, y_train)
+
+# Prawdopodobieństwa prior
+priors = np.array(gnb.priors)
+counts = [np.sum(y_train == 0.0), np.sum(y_train == 1.0)]
+
+# Etykiety klas
+labels = ['Nowotwór złośliwy', 'Nowotwór łagodny']
+
+# Tworzenie diagramu kołowego
+fig, axes = plt.subplots(1, 2, figsize=(11, 6), facecolor='white')
+
+wedges, _, _ = axes[0].pie(priors, autopct='%1.1f%%', startangle=90, colors=['skyblue', 'orange'])
+axes[0].axis('equal') # Równe proporcje, aby koło było kołem
+axes[0].set_title('Prawdopodobieństwa prior')
+
+# Funkcja zwracająca liczbę dla autopct
+def func(pct, allvals):
+ absolute = int(np.round(pct/100.*np.sum(allvals)))
+ return f"{absolute}"
+
+axes[1].pie(counts, autopct=lambda pct: func(pct, counts), startangle=90, colors=['skyblue', 'orange'])
+axes[1].axis('equal') # Równe proporcje, aby koło było kołem
+axes[1].set_title('Liczba wystąpień każdej klasy')
+
+# Dodanie jednej legendy dla obu wykresów
+fig.legend(wedges, labels, title="Klasy", loc="center right", bbox_to_anchor=(0.6, 0.75))
+
+# Ustawienie białego tła dla wykresów
+fig.patch.set_facecolor('white')
+
+# Wyświetlenie wykresu
+plt.show()
+
print(f"Prawdopodobieństwa prior: {gnb.priors}")
+score = gnb.score(X_test.values, y_test)
+print("Celność na zbiorze testowym:", score)
+
+# Predykcje
+preds = gnb.predict(X_test.values[15:20])
+probs = gnb.predict_proba(X_test.values[15:20])
+pp_predictions(preds, probs)
+
Prawdopodobieństwa prior: [0.38085938 0.61914062] +Celność na zbiorze testowym: 0.9824561403508771 + +=== 1 === +Predykcja: 1 - (Nowotwór łagodny) +Prawdpodobieństwo (Złośliwy, Łagodny): [0. 1.] + +=== 2 === +Predykcja: 1 - (Nowotwór łagodny) +Prawdpodobieństwo (Złośliwy, Łagodny): [0.1798 0.8202] + +=== 3 === +Predykcja: 1 - (Nowotwór łagodny) +Prawdpodobieństwo (Złośliwy, Łagodny): [0.0067 0.9933] + +=== 4 === +Predykcja: 1 - (Nowotwór łagodny) +Prawdpodobieństwo (Złośliwy, Łagodny): [0.0472 0.9528] + +=== 5 === +Predykcja: 0 - (Nowotwór złośliwy) +Prawdpodobieństwo (Złośliwy, Łagodny): [1. 0.] ++
plt.figure(figsize=(14, 14))
+score = gnb.score_visualize(X_test.values, y_test)
+
<Figure size 1008x1008 with 0 Axes>+
Porównanie wyników z gotową implementacją z użyciem biblioteki sklearn¶
+from sklearn.naive_bayes import GaussianNB
+
+clf = GaussianNB()
+clf.fit(X_train.values, y_train)
+print("Celność na zbiorze testowym:", clf.score(X_test.values, y_test))
+preds = clf.predict(X_test.values[15:20])
+probs = clf.predict_proba(X_test.values[15:20])
+pp_predictions(preds, probs)
+
Celność na zbiorze testowym: 0.9824561403508771 + +=== 1 === +Predykcja: 1.0 - (Nowotwór łagodny) +Prawdpodobieństwo (Złośliwy, Łagodny): [0. 1.] + +=== 2 === +Predykcja: 1.0 - (Nowotwór łagodny) +Prawdpodobieństwo (Złośliwy, Łagodny): [0.1123 0.8877] + +=== 3 === +Predykcja: 1.0 - (Nowotwór łagodny) +Prawdpodobieństwo (Złośliwy, Łagodny): [0. 1.] + +=== 4 === +Predykcja: 1.0 - (Nowotwór łagodny) +Prawdpodobieństwo (Złośliwy, Łagodny): [0.047 0.953] + +=== 5 === +Predykcja: 0.0 - (Nowotwór złośliwy) +Prawdpodobieństwo (Złośliwy, Łagodny): [1. 0.] ++
Wizualizacja obszarów decyzyjnych¶
+from matplotlib.colors import ListedColormap
+from sklearn.decomposition import PCA
+
+# Wizualizacja obszarów decyzyjnych
+def plot_decision_boundaries(X, y, model, title, h=0.9):
+ x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
+ y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
+ xx, yy = np.meshgrid(np.arange(x_min, x_max, h),np.arange(y_min, y_max, h))
+
+ Z = model.predict(np.c_[xx.ravel(), yy.ravel()])
+ Z = Z.reshape(xx.shape)
+ plt.figure(figsize=(12, 9))
+ colors_background = ['#FFFFFF', '#00AACC']
+ colors_points = ['#0000FF', '#00FF00']
+ plt.contourf(xx, yy, Z, alpha=0.8, cmap=ListedColormap(colors_background))
+ plt.scatter(X[:, 0], X[:, 1], c=y, s=50, edgecolor='b', cmap=ListedColormap(colors_points))
+ plt.title(title)
+ plt.show()
+
+# 30 cech zostało zredukowane do dwóch
+# PCA dokonuje liniowej redukcji wymiarów za pomocą algorytmu Singular Value Decomposition (SVD).
+pca = PCA(n_components=2)
+X_pca = pca.fit_transform(X_train.values)
+
+gnb = GaussianNaiveBayes(log=True)
+gnb.fit(X_pca, y_train)
+plot_decision_boundaries(X_pca, y_train, gnb, "Wizualizacja obszarów decyzyjnych")
+