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