<jupyter_start><jupyter_text>Seminar 01: Naive Bayes from scratch Today we will write Naive Bayes classifier supporting different feature probabilities._Authors: [<NAME>](https://github.com/neychev), [<NAME>](https://github.com/v-goncharenko)_  Loading data<jupyter_code>import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import datasets
matplotlib.rcParams['font.size'] = 11<jupyter_output><empty_output><jupyter_text>First to load dataset we're going to use [`sklearn`](https://scikit-learn.org/stable/) package which we will extensively use during the whole course.`sklearn` implement most of classical and frequently used algorithms in Machine Learning. Also it provides [User Guide](https://scikit-learn.org/stable/user_guide.html) describing principles of every bunch of algorithms implemented.As an entry point to main `sklearn`'s concepts we recommend [getting started tutorial](https://scikit-learn.org/stable/getting_started.html) (check it out yourself). [Further tutorials](https://scikit-learn.org/stable/tutorial/index.html) can also be handy to develop your skills. First functionality we use is cosy loading of [common datasets](https://scikit-learn.org/stable/modules/classes.htmlmodule-sklearn.datasets). All we need to do is just one function call.Object generated by [`load_iris`](https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_iris.html) is described as:> Dictionary-like object, the interesting attributes are:>> ‘data’, the data to learn,>>‘target’, the classification labels,>>‘target_names’, the meaning of the labels,>>‘feature_names’, the meaning of the features,>>‘DESCR’, the full description of the dataset,>>‘filename’, the physical location of iris csv dataset (added in version 0.20)Let's see what we have<jupyter_code>dataset = datasets.load_iris()

print(dataset.DESCR)<jupyter_output>.. _iris_dataset:

Iris plants dataset
--------------------

**Data Set Characteristics:**

    :Number of Instances: 150 (50 in each of three classes)
    :Number of Attributes: 4 numeric, predictive attributes and the class
    :Attribute Information:
        - sepal length in cm
        - sepal width in cm
        - petal length in cm
        - petal width in cm
        - class:
                - Iris-Setosa
                - Iris-Versicolour
                - Iris-Virginica
                
    :Summary Statistics:

    ============== ==== ==== ======= ===== ====================
                    Min  Max   Mean    SD   Class Correlation
    ============== ==== ==== ======= ===== ====================
    sepal length:   4.3  7.9   5.84   0.83    0.7826
    sepal width:    2.0  4.4   3.05   0.43   -0.4194
    petal length:   1.0  6.9   3.76   1.76    0.9490  (high!)
    petal width:    0.1  2.5   1.20   0.76    0.9565  (high!)
    ============== ==== ==== ======= ===== ===========[...]<jupyter_text>If you aren't familiar with Iris dataset - take a minute to read description above =) (as always [more info about it in Wikipedia](https://en.wikipedia.org/wiki/Iris_flower_data_set))__TL;DR__ 150 objects equally distributed over 3 classes each described with 4 continuous featuresJust pretty table to look at:<jupyter_code># for now you don't need to understand what happens in this code - just look at the table
ext_target = dataset.target[:, None]
pd.DataFrame(
    np.concatenate((dataset.data, ext_target, dataset.target_names[ext_target]), axis=1),
    columns=dataset.feature_names + ['target label', 'target name'],
)<jupyter_output><empty_output><jupyter_text>Now give distinct names to the data we will use<jupyter_code>features = dataset.data
target = dataset.target

features.shape, target.shape<jupyter_output><empty_output><jupyter_text>__Please, remember!!!__Anywhere in our course we have an agreement to shape design matrix (named `features` in code above) as `(number_of_items, number_of_features)`  Visualize dataset Our dataset has 4 dimensions however humans are more common to 3 or even 2 dimensional data, so let's plot first 3 features colored with labels values<jupyter_code>from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(8, 8))

ax = Axes3D(fig)

ax.scatter(features[:, 0], features[:, 1], features[:, 3], c=target, marker='o')
ax.set_xlabel(dataset.feature_names[0])
ax.set_ylabel(dataset.feature_names[1])
ax.set_zlabel(dataset.feature_names[2])

plt.show()<jupyter_output><empty_output><jupyter_text>Then have a look on feature distributions<jupyter_code># remember this way to make subplots! It could be useful for you later in your work

fig, axes = plt.subplots(2, 2, figsize=(15, 10))

for i, axis in enumerate(axes.flat):
    axis.hist(features[:, i])
    axis.set_xlabel(dataset.feature_names[i])
    axis.set_ylabel('number of objects')<jupyter_output><empty_output><jupyter_text>Note that every plot above have own scale  Classifier implementation Since we aiming to implement Naive Bayes algorithm first we need some prior distribution defined.The most common distribution is (of course) Gaussian and it's params are mean and standard deviation. Let's implement class taking list of feature values, estimating distribution params and able to give probability density of any given feature value. Denote the normal distribution $\mathcal{N}(\mu, \sigma^2)$ PDF:$$f(x|\mu, \sigma^2) = \frac{1}{\sigma\sqrt{2\pi}}\exp(-\frac{(x - \mu)^2}{2\sigma^2})$$Let's implement the `GaussianDistribution` class. (Of course in practice one could always use something like `scipy.stats.norm`).Please note, that making computations with log probabilities is more stable.<jupyter_code>class GaussianDistribution:
    def __init__(self, feature):
        '''
        Args:
            feature: column of design matrix, represents all available values
                of feature to model.
                axis=0 stays for samples.
        '''
        self.mean = feature.mean(axis=0)
        self.std = feature.std(axis=0)

    def logpdf(self, value):
        '''Logarithm of probability density at value'''
        return # <YOUR CODE HERE>
    
    def pdf(self, value):
        return # <YOUR CODE HERE><jupyter_output><empty_output><jupyter_text>Let's check the result:<jupyter_code>import scipy
_test = scipy.stats.norm(loc=features[:, :2].mean(axis=0), scale=features[:, :2].std(axis=0))
assert np.allclose(
    GaussianDistribution(features[:, :2]).logpdf(features[:5, :2]),
    _test.logpdf(features[:5, :2])
)
print('Seems fine!')<jupyter_output>Seems fine!<jupyter_text>Let's focus on the classification problem now. For the case of $K$ classes label $y_i \in \{C_1, \ldots, C_k\}$. Iris classification problem has 3 classes, so $K=3$. Bayes' Theorem takes the following form:$$P(y_i = C_k|\mathbf{x}_i) = \frac{P(\mathbf{x}_i|y_i = C_k) P(y_i = C_k)}{P(\mathbf{x}_i)}$$ Please note, we prefer working with log probabilities here as well. So the equation above will take the following form:$$\log P(y_i = C_k|\mathbf{x}_i) = \log P(\mathbf{x}_i|y_i = C_k) + \log P(y_i = C_k) - \log P(\mathbf{x}_i)$$As one could mention, to find the class label with the highest probability we even do not need the last term $P(\mathbf{x}_i)$. However, we need it to get the correct estimation of the probability $P(y_i = C_k|\mathbf{x}_i)$. The $P(\mathbf{x}_i)$ term can be computed using the following property:$$P(\mathbf{x}_i) = \sum_{k=1}^K P(\mathbf{x}_i|y_i=C_k).$$It can be computed from $\log P(\mathbf{x}_i|y_i=C_k)$ values using `logsumexp` function located in `scipy.special`.Now let's implement the Naive Bayes classifier itself. The class below is inherited from `sklearn` base classes and provides all the main methods.<jupyter_code>from sklearn.base import BaseEstimator, ClassifierMixin
from scipy.special import logsumexp


class NaiveBayes(BaseEstimator, ClassifierMixin):
    '''
    Please note, using `X` and `y` for design matrix and labels in general is not a good choice,
    better stick to more informative naming conventions.
    However, to make the code consistent with sklearn implementation, we use `X` and `y` variables here.
    '''
    def fit(self, X, y, sample_weight=None, distributions=None):
        '''
        sample_weight 
            The argument is ignored. For comatibility only.
        '''
        self.unique_labels = np.unique(y)
        
        # If distributions of features are not specified, they a treated Gaussian
        if distributions is None:
            distributions = [GaussianDistribution] * X.shape[1]
        else:
            # Check whether distributions are passed for all features
            assert len(distributions) == X.shape[1]        

        # Here we find distribution parameters for every feature in every class subset
        # so P(x^i|y=C_k) will be estimated only using information from i-th feature of C_k class values
        self.conditional_feature_distributions = {} # label: [distibution for feature 1, ...]
        for label in self.unique_labels:
            feature_distribution = []
            for column_index in range(X.shape[1]):
                # `column_index` feature values for objects from `label` class
                feature_column = X[y == label, column_index]
                fitted_distr = distributions[column_index](feature_column)
                feature_distribution.append(fitted_distr)
            self.conditional_feature_distributions[label] = feature_distribution

        # Prior label distributions (unconditional probability of each class)
        self.prior_label_distibution = {
            # <YOUR CODE HERE>
        }

    def predict_log_proba(self, X):
        # Matrix of shape (n_objects : n_classes)
        class_log_probas = np.zeros((X.shape[0], len(self.unique_labels)), dtype=float)
        
        # Here we compute the class log probabilities for each class sequentially b
        for label_idx, label in enumerate(self.unique_labels):
            for idx in range(X.shape[1]):
                # All loglikelihood for every feature w.r.t. fixed label
                class_log_probas[:, label_idx] += # <YOUR CODE HERE>

            # Add log proba of label prior
            class_log_probas[:, label_idx] += # <YOUR CODE HERE>

        for idx in range(X.shape[1]):
        # If you want to get probabilities, you need to substract the log proba for every feature
            class_log_probas -= # <YOUR CODE HERE>
        return class_log_probas
    
    def predict_proba(self, X):
        return np.exp(self.predict_log_proba(X))
    
    def predict(self, X):
        log_probas = self.predict_log_proba(X)
        # we need to cast labels to their original form (they may start from number other than 0)
        return np.array([self.unique_labels[idx] for idx in log_probas.argmax(axis=1)])
nb = NaiveBayes()
nb.fit(features, target)
print('log probas:\n{}'.format(nb.predict_log_proba(features[:2])))
print('predicted labels:\n{}'.format(nb.predict(features[:2])))
print('\nIt`s alive! More tests coming.')<jupyter_output>log probas:
[[  0.         -41.14063634 -57.90531295]
 [  0.         -38.72865982 -56.71097942]]
predicted labels:
[0 0]

It`s alive! More tests coming.<jupyter_text>Now let's check our Naive Bayes classifier on the unseed data. To do so we will use `train_test_split` from `sklearn`.<jupyter_code>from sklearn.model_selection import train_test_split
features_train, features_test, target_train, target_test = train_test_split(features, target, test_size=0.25)

print(features_train.shape, features_test.shape)
nb = NaiveBayes()
nb.fit(features_train, target_train, distributions=[GaussianDistribution]*4)
nb_test_log_proba = nb.predict_log_proba(features_test)
print('Naive Bayes classifier accuracy on the train set: {}'.format(nb.score(features_train, target_train)))
print('Naive Bayes classifier accuracy on the test set: {}'.format(nb.score(features_test, target_test)))<jupyter_output>Naive Bayes classifier accuracy on the test set: 0.9473684210526315<jupyter_text>Finally, let's comapre the Naive Bayes classifier with the `sklearn` implementations.<jupyter_code>from sklearn import naive_bayes

sklearn_nb = naive_bayes.GaussianNB()
sklearn_nb.fit(features_train, target_train)
sklearn_nb_test_log_proba = sklearn_nb.predict_log_proba(features_test)
print('sklearn implementation accuracy on the train set: {}'.format(sklearn_nb.score(features_train, target_train)))
print('sklearn implementation accuracy on the test set: {}'.format(sklearn_nb.score(features_test, target_test)))<jupyter_output>sklearn implementation accuracy on the test set: 0.9473684210526315<jupyter_text>And let's even check the predictions. If you used Gaussian distribution and done everything correctly, the log probabilities should be the same.<jupyter_code>assert np.allclose(nb_test_log_proba, sklearn_nb_test_log_proba), 'log probabilities do not match'
print('Seems alright!')<jupyter_output>Seems alright!<jupyter_text>Comparing to kNN<jupyter_code>from sklearn.neighbors import KNeighborsClassifier
knn = KNeighborsClassifier(n_neighbors=3)
knn.fit(features_train, target_train)
knn.score(features_train, target_train)
knn.score(features_test, target_test)<jupyter_output><empty_output><jupyter_text>Seem like Naive Bayes classifier performance is comparable to the kNN, while Naive Bayes does not need to store all the train data (while kNN need).  Advanced distribution for NaiveBayesLet's take a look at violin plots for every feature in our dataset:<jupyter_code>plt.figure(figsize=(15, 15))
plt.violinplot(features, showmedians=True)<jupyter_output><empty_output><jupyter_text>Although we do love Gaussian distribution it is still unimodal while our features are substantially multimodal (see histograms above). So we have to implement more robust distribution estimator - Kernel Density Estimator (KDE).Idea for this method is simple: we assign some probability density to a region around actual observation. (We will return to density estimation methods to describe them carefully later in this course).Fortunately `sklearn` have KDE implemented for us already. All it needs is vector of feature values. To get probability estimations using KDE one can easily access the `sklearn.neighbors` module.<jupyter_code>from sklearn.neighbors import KernelDensity
kde = KernelDensity(bandwidth=0.28, kernel='gaussian')
feature_col = features[target==2, 2]
kde.fit(feature_col.reshape((-1, 1)))
linspace = np.linspace(feature_col.min(), feature_col.max(), 1000)
plt.plot(linspace, np.exp(kde.score_samples(linspace.reshape((-1, 1)))))
plt.grid()
plt.xlabel('feature value')
plt.ylabel('probability')<jupyter_output><empty_output><jupyter_text>To make it compatible with the Naive Bayes classifier we have implemented above, we need to create class with the same methods:<jupyter_code>class GaussianKDE:
    def __init__(self, feature):
        self.kde = KernelDensity(bandwidth=1.)
        self.kde.fit(feature.reshape((-1, 1)))

    def logpdf(self, value):
        return self.kde.score_samples(value.reshape((-1, 1)))

    def pdf(self, value):
        return np.exp(self.log_proba(value))
nb_kde = NaiveBayes()
nb_kde.fit(features, target, distributions=[GaussianKDE]*4)
print('log probas:\n{}'.format(nb_kde.predict_log_proba(features[:2])))
print('predicted labels:\n{}'.format(nb_kde.predict(features[:2])))
print('\nIt`s alive!')
print('KDE Naive Bayes classifier accuracy on the train set: {}'.format(nb_kde.score(features_train, target_train)))
print('KDE Naive Bayes classifier accuracy on the test set: {}'.format(nb_kde.score(features_test, target_test)))<jupyter_output>KDE Naive Bayes classifier accuracy on the test set: 0.9736842105263158<jupyter_text>Seems like the accuracy of the classifier has decreased. What is going on?_Hint: try varying the `bandwidth` parameter of the `KernelDensity` constructor in `GaussianKDE` class (around 0.3)._Let's take a closer look on the features distributions. Here comes the histogram:<jupyter_code>fig, axes = plt.subplots(2, 3, figsize=(15, 10))

for ax_idx, feature_idx in enumerate([2, 3]):
    for label in range(3):
        ax = axes[ax_idx, label]
        feature_col = features[target==label, feature_idx]
        ax.hist(feature_col, bins=7)
        ax.grid()
        ax.set_title('class: {}, feature: {}'.format(
            dataset.target_names[label],
            dataset.feature_names[feature_idx]
        ))<jupyter_output><empty_output><jupyter_text>We see, than the distributions within every class are unimodal. That's how KDE is approximating the PDF:<jupyter_code>fig, axes = plt.subplots(2, 3, figsize=(15, 10))
kde = KernelDensity(bandwidth=1., kernel='gaussian')

for ax_idx, feature_idx in enumerate([2, 3]):
    for label in range(3):
        ax = axes[ax_idx, label]
        feature_col = features[target==label, feature_idx]
        kde.fit(feature_col.reshape((-1, 1)))
        linspace = np.linspace(
            0.8*feature_col.min(),
            1.2*feature_col.max(),
            1000
        )
        ax.plot(linspace, np.exp(kde.score_samples(linspace.reshape((-1, 1)))))
        ax.grid()
        ax.set_title('class: {}, feature: {}'.format(
            dataset.target_names[label],
            dataset.feature_names[feature_idx]
        ))<jupyter_output><empty_output><jupyter_text>One could mention, that every feature need different `bandwidth` parameter.And that's how Gaussian distribution fits to the data:<jupyter_code>fig, axes = plt.subplots(2, 3, figsize=(15, 10))

for ax_idx, feature_idx in enumerate([2, 3]):
    for label in range(3):
        ax = axes[ax_idx, label]
        feature_col = features[target==label, feature_idx]
        gaussian_distr = GaussianDistribution(feature_col)
        linspace = np.linspace(
            feature_col.min(),
            feature_col.max(),
            1000
        )
        ax.plot(linspace, gaussian_distr.pdf(linspace.reshape((-1, 1))))
        ax.grid()
        ax.set_title('class: {}, feature: {}'.format(
            dataset.target_names[label],
            dataset.feature_names[feature_idx]
        ))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11330022075055188, 'fraction_numerical': 0.016997792494481235, 'mean_word_length': 4.113148984198646, 'pattern_counts': {'":': 0, '<': 78, '<?xml version=': 0, '>': 91, 'https://': 9, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '19', 'dir': 'jupyter-structured-clean-dedup', 'id': '80262', 'license': 'MIT', 'n_tokens_mistral': 5780, 'n_tokens_neox': 5446, 'n_words': 1661, 'path': 'week0_01_org_kNN_and_Naive_Bayes/week0_01_naive_bayes_refined.ipynb', 'repo_name': 'FabulousFatTiger/ml-mipt', 'size': '18145'}