08 Naive Bayes — First Principles¶

Goal: Build a Gaussian Naive Bayes classifier from scratch using only NumPy, verify it against sklearn and the src/ library implementation, and explore when the naive independence assumption helps and hurts.

Prerequisites: Probability & Statistics, theory.md

Outline:

1. Problem Setup — WHY
2. Mathematical Core — WHAT
3. Solution Method — HOW
4. Implementation — BUILD
5. Library Comparison
6. Experiments and Failures — VERIFY
7. Connections

In [1]:

import random

import matplotlib.pyplot as plt
import numpy as np

%load_ext autoreload
%autoreload 2

SEED = 42
random.seed(SEED)
rng = np.random.default_rng(SEED)

import random import matplotlib.pyplot as plt import numpy as np %load_ext autoreload %autoreload 2 SEED = 42 random.seed(SEED) rng = np.random.default_rng(SEED)

1. Problem Setup — WHY¶

The Spam Classification Intuition¶

Imagine classifying emails as spam or not-spam based on word presence. With thousands of word features and limited labeled data, flexible models overfit. Naive Bayes makes a radical simplification — treat every word as independent given the class — and trades model fidelity for dramatic variance reduction.

Let's visualize the core idea with a 2D continuous example first.

In [2]:

# Generate two Gaussian classes with known parameters
n_per_class = 100
X_class0 = rng.normal(loc=[1.0, 2.0], scale=[0.8, 0.6], size=(n_per_class, 2))
X_class1 = rng.normal(loc=[3.0, 4.0], scale=[1.0, 0.8], size=(n_per_class, 2))

X_demo = np.vstack([X_class0, X_class1])
y_demo = np.array([0] * n_per_class + [1] * n_per_class)

fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# Scatter plot
for c, color, label in [(0, 'tab:blue', 'Class 0'), (1, 'tab:red', 'Class 1')]:
    mask = y_demo == c
    axes[0].scatter(X_demo[mask, 0], X_demo[mask, 1], c=color, alpha=0.5, label=label)
axes[0].set(title='2D Classification Problem', xlabel='Feature 1', ylabel='Feature 2')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Feature 1 distributions per class
for c, color, label in [(0, 'tab:blue', 'Class 0'), (1, 'tab:red', 'Class 1')]:
    mask = y_demo == c
    axes[1].hist(X_demo[mask, 0], bins=20, alpha=0.5, color=color, label=label, density=True)
axes[1].set(title='Feature 1 | Class', xlabel='Feature 1', ylabel='Density')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# Feature 2 distributions per class
for c, color, label in [(0, 'tab:blue', 'Class 0'), (1, 'tab:red', 'Class 1')]:
    mask = y_demo == c
    axes[2].hist(X_demo[mask, 1], bins=20, alpha=0.5, color=color, label=label, density=True)
axes[2].set(title='Feature 2 | Class', xlabel='Feature 2', ylabel='Density')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print('Key insight: Naive Bayes models each feature distribution SEPARATELY per class,')
print('then multiplies them (or adds log-densities) to get the class posterior.')

# Generate two Gaussian classes with known parameters n_per_class = 100 X_class0 = rng.normal(loc=[1.0, 2.0], scale=[0.8, 0.6], size=(n_per_class, 2)) X_class1 = rng.normal(loc=[3.0, 4.0], scale=[1.0, 0.8], size=(n_per_class, 2)) X_demo = np.vstack([X_class0, X_class1]) y_demo = np.array([0] * n_per_class + [1] * n_per_class) fig, axes = plt.subplots(1, 3, figsize=(15, 4)) # Scatter plot for c, color, label in [(0, 'tab:blue', 'Class 0'), (1, 'tab:red', 'Class 1')]: mask = y_demo == c axes[0].scatter(X_demo[mask, 0], X_demo[mask, 1], c=color, alpha=0.5, label=label) axes[0].set(title='2D Classification Problem', xlabel='Feature 1', ylabel='Feature 2') axes[0].legend() axes[0].grid(True, alpha=0.3) # Feature 1 distributions per class for c, color, label in [(0, 'tab:blue', 'Class 0'), (1, 'tab:red', 'Class 1')]: mask = y_demo == c axes[1].hist(X_demo[mask, 0], bins=20, alpha=0.5, color=color, label=label, density=True) axes[1].set(title='Feature 1 | Class', xlabel='Feature 1', ylabel='Density') axes[1].legend() axes[1].grid(True, alpha=0.3) # Feature 2 distributions per class for c, color, label in [(0, 'tab:blue', 'Class 0'), (1, 'tab:red', 'Class 1')]: mask = y_demo == c axes[2].hist(X_demo[mask, 1], bins=20, alpha=0.5, color=color, label=label, density=True) axes[2].set(title='Feature 2 | Class', xlabel='Feature 2', ylabel='Density') axes[2].legend() axes[2].grid(True, alpha=0.3) plt.tight_layout() plt.show() print('Key insight: Naive Bayes models each feature distribution SEPARATELY per class,') print('then multiplies them (or adds log-densities) to get the class posterior.')

Key insight: Naive Bayes models each feature distribution SEPARATELY per class,
then multiplies them (or adds log-densities) to get the class posterior.

2. Mathematical Core — WHAT¶

Bayes' Theorem¶

$$P(y = k \mid x) = \frac{P(x \mid y = k) \, P(y = k)}{P(x)}$$

Naive Independence Assumption¶

$$P(x \mid y = k) = \prod_{j=1}^{d} P(x_j \mid y = k)$$

Gaussian Feature Likelihood¶

$$P(x_j \mid y = k) = \frac{1}{\sqrt{2\pi \sigma_{kj}^2}} \exp\!\left(-\frac{(x_j - \mu_{kj})^2}{2\sigma_{kj}^2}\right)$$

Decision Rule (Log-Space)¶

$$\hat{y} = \arg\max_k \left[ \log P(y = k) + \sum_{j=1}^{d} \log P(x_j \mid y = k) \right]$$

Let's verify Bayes' theorem numerically on a tiny example.

In [3]:

# Tiny example: 1 feature, 2 classes, Gaussian likelihoods
# Class 0: mu=1, sigma=0.5, prior=0.6
# Class 1: mu=3, sigma=1.0, prior=0.4
from scipy.stats import norm

prior = np.array([0.6, 0.4])
mu = np.array([1.0, 3.0])
sigma = np.array([0.5, 1.0])

x_test = 2.0

# P(x | y=k) for each class
likelihood = np.array([norm.pdf(x_test, mu[k], sigma[k]) for k in range(2)])

# P(x) = sum_k P(x|y=k) P(y=k)
evidence = np.sum(likelihood * prior)

# P(y=k | x) = P(x|y=k) P(y=k) / P(x)
posterior = likelihood * prior / evidence

print('Bayes\' Theorem — step by step')
print(f'  x_test = {x_test}')
print(f'  P(x | y=0) = {likelihood[0]:.4f},  P(x | y=1) = {likelihood[1]:.4f}')
print(f'  P(y=0) = {prior[0]},  P(y=1) = {prior[1]}')
print(f'  P(x) = {evidence:.4f}')
print(f'  P(y=0 | x) = {posterior[0]:.4f},  P(y=1 | x) = {posterior[1]:.4f}')
print(f'  Sum of posteriors = {posterior.sum():.4f}')

assert np.isclose(posterior.sum(), 1.0, atol=1e-12)
print('\nPosteriors sum to 1 ✓')

# Tiny example: 1 feature, 2 classes, Gaussian likelihoods # Class 0: mu=1, sigma=0.5, prior=0.6 # Class 1: mu=3, sigma=1.0, prior=0.4 from scipy.stats import norm prior = np.array([0.6, 0.4]) mu = np.array([1.0, 3.0]) sigma = np.array([0.5, 1.0]) x_test = 2.0 # P(x | y=k) for each class likelihood = np.array([norm.pdf(x_test, mu[k], sigma[k]) for k in range(2)]) # P(x) = sum_k P(x|y=k) P(y=k) evidence = np.sum(likelihood * prior) # P(y=k | x) = P(x|y=k) P(y=k) / P(x) posterior = likelihood * prior / evidence print('Bayes\' Theorem — step by step') print(f' x_test = {x_test}') print(f' P(x | y=0) = {likelihood[0]:.4f}, P(x | y=1) = {likelihood[1]:.4f}') print(f' P(y=0) = {prior[0]}, P(y=1) = {prior[1]}') print(f' P(x) = {evidence:.4f}') print(f' P(y=0 | x) = {posterior[0]:.4f}, P(y=1 | x) = {posterior[1]:.4f}') print(f' Sum of posteriors = {posterior.sum():.4f}') assert np.isclose(posterior.sum(), 1.0, atol=1e-12) print('\nPosteriors sum to 1 ✓')

Bayes' Theorem — step by step
  x_test = 2.0
  P(x | y=0) = 0.1080,  P(x | y=1) = 0.2420
  P(y=0) = 0.6,  P(y=1) = 0.4
  P(x) = 0.1616
  P(y=0 | x) = 0.4010,  P(y=1 | x) = 0.5990
  Sum of posteriors = 1.0000

Posteriors sum to 1 ✓

3. Solution Method — HOW¶

Training (Parameter Estimation)¶

1. Priors: $\hat{\pi}_k = n_k / n$
2. Means: $\hat{\mu}_{kj} = \frac{1}{n_k} \sum_{i: y_i=k} x_{ij}$
3. Variances: $\hat{\sigma}^2_{kj} = \frac{1}{n_k} \sum_{i: y_i=k} (x_{ij} - \hat{\mu}_{kj})^2 + \epsilon$

Prediction¶

For each test point $x$, compute the log-posterior for every class:

$$a_k = \log \hat{\pi}_k + \sum_{j=1}^{d} \left[ -\frac{1}{2} \log(2\pi \hat{\sigma}^2_{kj}) - \frac{(x_j - \hat{\mu}_{kj})^2}{2\hat{\sigma}^2_{kj}} \right]$$

Predict: $\hat{y} = \arg\max_k \, a_k$

No iteration. No gradient descent. One pass through data.

4. Implementation — BUILD¶

Build a Gaussian Naive Bayes classifier from scratch.

In [4]:

class GaussianNBScratch:
    """Gaussian Naive Bayes from scratch — log-space implementation."""

    def __init__(self, var_smoothing: float = 1e-9):
        self.var_smoothing = var_smoothing
        self.classes_ = None
        self.mean_ = None       # shape (n_classes, n_features)
        self.var_ = None        # shape (n_classes, n_features)
        self.class_log_prior_ = None  # shape (n_classes,)

    def fit(self, X, y):
        """Estimate class priors and per-class Gaussian parameters."""
        X = np.asarray(X, dtype=float)
        y = np.asarray(y)

        self.classes_ = np.unique(y)
        n_classes = len(self.classes_)
        n_features = X.shape[1]

        self.mean_ = np.zeros((n_classes, n_features))
        self.var_ = np.zeros((n_classes, n_features))
        counts = np.zeros(n_classes)

        # Variance smoothing: small epsilon relative to largest variance
        epsilon = self.var_smoothing * max(float(np.var(X, axis=0).max()), 1.0)

        for idx, c in enumerate(self.classes_):
            X_c = X[y == c]
            counts[idx] = X_c.shape[0]
            self.mean_[idx] = X_c.mean(axis=0)
            self.var_[idx] = X_c.var(axis=0) + epsilon  # population variance + smoothing

        self.class_log_prior_ = np.log(counts / X.shape[0])
        return self

    def _joint_log_likelihood(self, X):
        """Compute log P(y=k) + sum_j log P(x_j | y=k) for all classes.

        Returns array of shape (n_samples, n_classes).
        """
        # X: (n, d) -> (n, 1, d)
        # mean_: (K, d) -> (1, K, d)
        diff = X[:, None, :] - self.mean_[None, :, :]

        # Log Gaussian density: -0.5 * [log(2*pi*var) + (x-mu)^2 / var]
        log_density = -0.5 * (
            np.log(2.0 * np.pi * self.var_)[None, :, :]
            + diff ** 2 / self.var_[None, :, :]
        )

        # Sum over features, add log prior
        # Result shape: (n_samples, n_classes)
        return self.class_log_prior_[None, :] + log_density.sum(axis=2)

    def predict(self, X):
        """Predict class with highest joint log-likelihood."""
        X = np.asarray(X, dtype=float)
        jll = self._joint_log_likelihood(X)
        return self.classes_[np.argmax(jll, axis=1)]

    def score(self, X, y):
        """Return classification accuracy."""
        return float(np.mean(self.predict(X) == np.asarray(y)))

class GaussianNBScratch: """Gaussian Naive Bayes from scratch — log-space implementation.""" def __init__(self, var_smoothing: float = 1e-9): self.var_smoothing = var_smoothing self.classes_ = None self.mean_ = None # shape (n_classes, n_features) self.var_ = None # shape (n_classes, n_features) self.class_log_prior_ = None # shape (n_classes,) def fit(self, X, y): """Estimate class priors and per-class Gaussian parameters.""" X = np.asarray(X, dtype=float) y = np.asarray(y) self.classes_ = np.unique(y) n_classes = len(self.classes_) n_features = X.shape[1] self.mean_ = np.zeros((n_classes, n_features)) self.var_ = np.zeros((n_classes, n_features)) counts = np.zeros(n_classes) # Variance smoothing: small epsilon relative to largest variance epsilon = self.var_smoothing * max(float(np.var(X, axis=0).max()), 1.0) for idx, c in enumerate(self.classes_): X_c = X[y == c] counts[idx] = X_c.shape[0] self.mean_[idx] = X_c.mean(axis=0) self.var_[idx] = X_c.var(axis=0) + epsilon # population variance + smoothing self.class_log_prior_ = np.log(counts / X.shape[0]) return self def _joint_log_likelihood(self, X): """Compute log P(y=k) + sum_j log P(x_j | y=k) for all classes. Returns array of shape (n_samples, n_classes). """ # X: (n, d) -> (n, 1, d) # mean_: (K, d) -> (1, K, d) diff = X[:, None, :] - self.mean_[None, :, :] # Log Gaussian density: -0.5 * [log(2*pi*var) + (x-mu)^2 / var] log_density = -0.5 * ( np.log(2.0 * np.pi * self.var_)[None, :, :] + diff ** 2 / self.var_[None, :, :] ) # Sum over features, add log prior # Result shape: (n_samples, n_classes) return self.class_log_prior_[None, :] + log_density.sum(axis=2) def predict(self, X): """Predict class with highest joint log-likelihood.""" X = np.asarray(X, dtype=float) jll = self._joint_log_likelihood(X) return self.classes_[np.argmax(jll, axis=1)] def score(self, X, y): """Return classification accuracy.""" return float(np.mean(self.predict(X) == np.asarray(y)))

In [5]:

# Quick sanity check on demo data
gnb_scratch = GaussianNBScratch(var_smoothing=1e-9)
gnb_scratch.fit(X_demo, y_demo)

preds = gnb_scratch.predict(X_demo)
acc = gnb_scratch.score(X_demo, y_demo)
print(f'Training accuracy: {acc:.4f}')
assert acc > 0.90, f'Expected accuracy > 0.90, got {acc}'

# Check learned parameters are close to true values
print(f'\nLearned means:\n{gnb_scratch.mean_}')
print(f'True means: Class 0 = [1.0, 2.0], Class 1 = [3.0, 4.0]')
print(f'\nLearned log-priors: {gnb_scratch.class_log_prior_}')
print(f'Expected log-priors (balanced): [{np.log(0.5):.4f}, {np.log(0.5):.4f}]')

assert np.isclose(gnb_scratch.class_log_prior_[0], gnb_scratch.class_log_prior_[1], atol=1e-10)
print('\nSanity checks passed ✓')

# Quick sanity check on demo data gnb_scratch = GaussianNBScratch(var_smoothing=1e-9) gnb_scratch.fit(X_demo, y_demo) preds = gnb_scratch.predict(X_demo) acc = gnb_scratch.score(X_demo, y_demo) print(f'Training accuracy: {acc:.4f}') assert acc > 0.90, f'Expected accuracy > 0.90, got {acc}' # Check learned parameters are close to true values print(f'\nLearned means:\n{gnb_scratch.mean_}') print(f'True means: Class 0 = [1.0, 2.0], Class 1 = [3.0, 4.0]') print(f'\nLearned log-priors: {gnb_scratch.class_log_prior_}') print(f'Expected log-priors (balanced): [{np.log(0.5):.4f}, {np.log(0.5):.4f}]') assert np.isclose(gnb_scratch.class_log_prior_[0], gnb_scratch.class_log_prior_[1], atol=1e-10) print('\nSanity checks passed ✓')

Training accuracy: 0.9850

Learned means:
[[0.98820304 1.97230685]
 [3.05883416 3.98505686]]
True means: Class 0 = [1.0, 2.0], Class 1 = [3.0, 4.0]

Learned log-priors: [-0.69314718 -0.69314718]
Expected log-priors (balanced): [-0.6931, -0.6931]

Sanity checks passed ✓

Test on the Iris Dataset¶

The classic Iris dataset has 150 samples, 4 features, 3 classes — a perfect test for multi-class Gaussian NB.

In [6]:

from sklearn.datasets import load_iris
from ml_first_principles.data_utils import train_test_split

iris = load_iris()
X_iris, y_iris = iris.data, iris.target
X_train, X_test, y_train, y_test = train_test_split(X_iris, y_iris, test_size=0.3, random_state=SEED)

gnb_iris = GaussianNBScratch(var_smoothing=1e-9)
gnb_iris.fit(X_train, y_train)

train_acc = gnb_iris.score(X_train, y_train)
test_acc = gnb_iris.score(X_test, y_test)
print(f'Iris — Scratch GaussianNB')
print(f'  Train accuracy: {train_acc:.4f}')
print(f'  Test accuracy:  {test_acc:.4f}')
assert test_acc > 0.85, f'Expected test accuracy > 0.85, got {test_acc}'

from sklearn.datasets import load_iris from ml_first_principles.data_utils import train_test_split iris = load_iris() X_iris, y_iris = iris.data, iris.target X_train, X_test, y_train, y_test = train_test_split(X_iris, y_iris, test_size=0.3, random_state=SEED) gnb_iris = GaussianNBScratch(var_smoothing=1e-9) gnb_iris.fit(X_train, y_train) train_acc = gnb_iris.score(X_train, y_train) test_acc = gnb_iris.score(X_test, y_test) print(f'Iris — Scratch GaussianNB') print(f' Train accuracy: {train_acc:.4f}') print(f' Test accuracy: {test_acc:.4f}') assert test_acc > 0.85, f'Expected test accuracy > 0.85, got {test_acc}'

Iris — Scratch GaussianNB
  Train accuracy: 0.9714
  Test accuracy:  0.9333

5. Library Comparison¶

Compare against:

1. sklearn.naive_bayes.GaussianNB
2. ml_first_principles.probabilistic_models.GaussianNB (the src/ library)

In [7]:

from sklearn.naive_bayes import GaussianNB as SklearnGNB

sk_gnb = SklearnGNB(var_smoothing=1e-9)
sk_gnb.fit(X_train, y_train)

sk_train_acc = sk_gnb.score(X_train, y_train)
sk_test_acc = sk_gnb.score(X_test, y_test)
print(f'Iris — sklearn GaussianNB')
print(f'  Train accuracy: {sk_train_acc:.4f}')
print(f'  Test accuracy:  {sk_test_acc:.4f}')

# Match predictions
scratch_preds = gnb_iris.predict(X_test)
sk_preds = sk_gnb.predict(X_test)
pred_match = np.mean(scratch_preds == sk_preds)
print(f'\nPrediction match rate: {pred_match:.4f}')

from sklearn.naive_bayes import GaussianNB as SklearnGNB sk_gnb = SklearnGNB(var_smoothing=1e-9) sk_gnb.fit(X_train, y_train) sk_train_acc = sk_gnb.score(X_train, y_train) sk_test_acc = sk_gnb.score(X_test, y_test) print(f'Iris — sklearn GaussianNB') print(f' Train accuracy: {sk_train_acc:.4f}') print(f' Test accuracy: {sk_test_acc:.4f}') # Match predictions scratch_preds = gnb_iris.predict(X_test) sk_preds = sk_gnb.predict(X_test) pred_match = np.mean(scratch_preds == sk_preds) print(f'\nPrediction match rate: {pred_match:.4f}')

Iris — sklearn GaussianNB
  Train accuracy: 0.9714
  Test accuracy:  0.9333

Prediction match rate: 1.0000

In [8]:

# Parameter matching: means, variances, and priors
print('=== Parameter Comparison (Scratch vs sklearn) ===')
print(f'\nClass log-priors:')
print(f'  Scratch: {gnb_iris.class_log_prior_}')
print(f'  sklearn: {np.log(sk_gnb.class_prior_)}')
# sklearn's GaussianNB exposes class_prior_ (probabilities), not log-priors
assert np.allclose(gnb_iris.class_log_prior_, np.log(sk_gnb.class_prior_), atol=1e-10), \
    'Log-priors do not match'
print('  ✓ Match')

print(f'\nClass means (first 2 features):')
print(f'  Scratch:\n{gnb_iris.mean_[:, :2]}')
print(f'  sklearn:\n{sk_gnb.theta_[:, :2]}')
assert np.allclose(gnb_iris.mean_, sk_gnb.theta_, atol=1e-10), \
    'Means do not match'
print('  ✓ Match')

print(f'\nClass variances (first 2 features):')
print(f'  Scratch:\n{gnb_iris.var_[:, :2]}')
print(f'  sklearn:\n{sk_gnb.var_[:, :2]}')
assert np.allclose(gnb_iris.var_, sk_gnb.var_, atol=1e-6), \
    'Variances do not match'
print('  ✓ Match')

# Parameter matching: means, variances, and priors print('=== Parameter Comparison (Scratch vs sklearn) ===') print(f'\nClass log-priors:') print(f' Scratch: {gnb_iris.class_log_prior_}') print(f' sklearn: {np.log(sk_gnb.class_prior_)}') # sklearn's GaussianNB exposes class_prior_ (probabilities), not log-priors assert np.allclose(gnb_iris.class_log_prior_, np.log(sk_gnb.class_prior_), atol=1e-10), \ 'Log-priors do not match' print(' ✓ Match') print(f'\nClass means (first 2 features):') print(f' Scratch:\n{gnb_iris.mean_[:, :2]}') print(f' sklearn:\n{sk_gnb.theta_[:, :2]}') assert np.allclose(gnb_iris.mean_, sk_gnb.theta_, atol=1e-10), \ 'Means do not match' print(' ✓ Match') print(f'\nClass variances (first 2 features):') print(f' Scratch:\n{gnb_iris.var_[:, :2]}') print(f' sklearn:\n{sk_gnb.var_[:, :2]}') assert np.allclose(gnb_iris.var_, sk_gnb.var_, atol=1e-6), \ 'Variances do not match' print(' ✓ Match')

=== Parameter Comparison (Scratch vs sklearn) ===

Class log-priors:
  Scratch: [-1.07044141 -1.07044141 -1.15745279]
  sklearn: [-1.07044141 -1.07044141 -1.15745279]
  ✓ Match

Class means (first 2 features):
  Scratch:
[[5.01111111 3.43333333]
 [5.94166667 2.73333333]
 [6.53030303 2.97575758]]
  sklearn:
[[5.01111111 3.43333333]
 [5.94166667 2.73333333]
 [6.53030303 2.97575758]]
  ✓ Match

Class variances (first 2 features):
  Scratch:
[[0.12765432 0.15111111]
 [0.30798611 0.09      ]
 [0.3057484  0.1072911 ]]
  sklearn:
[[0.12765432 0.15111111]
 [0.30798611 0.09      ]
 [0.3057484  0.1072911 ]]
  ✓ Match

In [9]:

# Compare with the src/ library implementation
from ml_first_principles.probabilistic_models import GaussianNB as LibGNB

lib_gnb = LibGNB(var_smoothing=1e-9)
lib_gnb.fit(X_train, y_train)

lib_test_acc = lib_gnb.score(X_test, y_test)
lib_preds = lib_gnb.predict(X_test)

print(f'Iris — src/ library GaussianNB')
print(f'  Test accuracy: {lib_test_acc:.4f}')

lib_match = np.mean(scratch_preds == lib_preds)
print(f'  Prediction match (scratch vs lib): {lib_match:.4f}')

# Parameters should match
assert np.allclose(gnb_iris.mean_, lib_gnb.mean_, atol=1e-10), 'Means do not match lib'
assert np.allclose(gnb_iris.var_, lib_gnb.var_, atol=1e-6), 'Variances do not match lib'
assert np.allclose(gnb_iris.class_log_prior_, lib_gnb.class_log_prior_, atol=1e-10), \
    'Log-priors do not match lib'
print('  Parameters match ✓')

# Compare with the src/ library implementation from ml_first_principles.probabilistic_models import GaussianNB as LibGNB lib_gnb = LibGNB(var_smoothing=1e-9) lib_gnb.fit(X_train, y_train) lib_test_acc = lib_gnb.score(X_test, y_test) lib_preds = lib_gnb.predict(X_test) print(f'Iris — src/ library GaussianNB') print(f' Test accuracy: {lib_test_acc:.4f}') lib_match = np.mean(scratch_preds == lib_preds) print(f' Prediction match (scratch vs lib): {lib_match:.4f}') # Parameters should match assert np.allclose(gnb_iris.mean_, lib_gnb.mean_, atol=1e-10), 'Means do not match lib' assert np.allclose(gnb_iris.var_, lib_gnb.var_, atol=1e-6), 'Variances do not match lib' assert np.allclose(gnb_iris.class_log_prior_, lib_gnb.class_log_prior_, atol=1e-10), \ 'Log-priors do not match lib' print(' Parameters match ✓')

Iris — src/ library GaussianNB
  Test accuracy: 0.9333
  Prediction match (scratch vs lib): 1.0000
  Parameters match ✓

6. Experiments and Failures — VERIFY¶

6.1 NB Works Despite Independence Violation¶

The independence assumption is almost always violated in practice, yet NB often classifies correctly. Let's demonstrate this with correlated features.

In [10]:

# Generate data with CORRELATED features (violates independence)
n = 200

# Class 0: correlated Gaussian (rho=0.8)
rho = 0.8
cov0 = np.array([[1.0, rho], [rho, 1.0]])
X0_corr = rng.multivariate_normal([0.0, 0.0], cov0, size=n)

# Class 1: correlated Gaussian (rho=0.8)
cov1 = np.array([[1.0, rho], [rho, 1.0]])
X1_corr = rng.multivariate_normal([3.0, 3.0], cov1, size=n)

X_corr = np.vstack([X0_corr, X1_corr])
y_corr = np.array([0] * n + [1] * n)

# Shuffle
idx = rng.permutation(len(y_corr))
X_corr, y_corr = X_corr[idx], y_corr[idx]

# Split
X_train_c, X_test_c, y_train_c, y_test_c = train_test_split(
    X_corr, y_corr, test_size=0.3, random_state=SEED
)

gnb_corr = GaussianNBScratch(var_smoothing=1e-9)
gnb_corr.fit(X_train_c, y_train_c)
acc_corr = gnb_corr.score(X_test_c, y_test_c)

print(f'Accuracy with correlated features (rho={rho}): {acc_corr:.4f}')
print(f'Correlation in class 0: {np.corrcoef(X0_corr.T)[0, 1]:.3f}')
print(f'\nDespite features being correlated (violating the naive assumption),')
print(f'NB still achieves {acc_corr*100:.1f}% accuracy because classification')
print(f'only needs the correct ranking of posterior probabilities.')

# Generate data with CORRELATED features (violates independence) n = 200 # Class 0: correlated Gaussian (rho=0.8) rho = 0.8 cov0 = np.array([[1.0, rho], [rho, 1.0]]) X0_corr = rng.multivariate_normal([0.0, 0.0], cov0, size=n) # Class 1: correlated Gaussian (rho=0.8) cov1 = np.array([[1.0, rho], [rho, 1.0]]) X1_corr = rng.multivariate_normal([3.0, 3.0], cov1, size=n) X_corr = np.vstack([X0_corr, X1_corr]) y_corr = np.array([0] * n + [1] * n) # Shuffle idx = rng.permutation(len(y_corr)) X_corr, y_corr = X_corr[idx], y_corr[idx] # Split X_train_c, X_test_c, y_train_c, y_test_c = train_test_split( X_corr, y_corr, test_size=0.3, random_state=SEED ) gnb_corr = GaussianNBScratch(var_smoothing=1e-9) gnb_corr.fit(X_train_c, y_train_c) acc_corr = gnb_corr.score(X_test_c, y_test_c) print(f'Accuracy with correlated features (rho={rho}): {acc_corr:.4f}') print(f'Correlation in class 0: {np.corrcoef(X0_corr.T)[0, 1]:.3f}') print(f'\nDespite features being correlated (violating the naive assumption),') print(f'NB still achieves {acc_corr*100:.1f}% accuracy because classification') print(f'only needs the correct ranking of posterior probabilities.')

Accuracy with correlated features (rho=0.8): 0.9583
Correlation in class 0: 0.811

Despite features being correlated (violating the naive assumption),
NB still achieves 95.8% accuracy because classification
only needs the correct ranking of posterior probabilities.

6.2 Failure Case — Duplicate Correlated Features¶

When features are near-duplicates, NB double-counts evidence, making posteriors overconfident. The accuracy may stay okay, but the probabilities become poorly calibrated.

In [11]:

# Create data where feature 2 is a noisy copy of feature 1
n_fail = 300
X0_base = rng.normal(0, 1, size=(n_fail, 1))
X1_base = rng.normal(2, 1, size=(n_fail, 1))

# Duplicate feature with small noise (near-perfect copy)
noise_level = 0.01
X0_dup = np.hstack([X0_base, X0_base + rng.normal(0, noise_level, size=(n_fail, 1))])
X1_dup = np.hstack([X1_base, X1_base + rng.normal(0, noise_level, size=(n_fail, 1))])

X_dup = np.vstack([X0_dup, X1_dup])
y_dup = np.array([0] * n_fail + [1] * n_fail)
idx_dup = rng.permutation(len(y_dup))
X_dup, y_dup = X_dup[idx_dup], y_dup[idx_dup]

X_tr_dup, X_te_dup, y_tr_dup, y_te_dup = train_test_split(
    X_dup, y_dup, test_size=0.3, random_state=SEED
)

# Train on 1 feature vs 2 duplicate features
gnb_1feat = GaussianNBScratch(var_smoothing=1e-9)
gnb_1feat.fit(X_tr_dup[:, :1], y_tr_dup)

gnb_2feat = GaussianNBScratch(var_smoothing=1e-9)
gnb_2feat.fit(X_tr_dup, y_tr_dup)

acc_1f = gnb_1feat.score(X_te_dup[:, :1], y_te_dup)
acc_2f = gnb_2feat.score(X_te_dup, y_te_dup)

print(f'Accuracy with 1 feature:           {acc_1f:.4f}')
print(f'Accuracy with 2 duplicate features: {acc_2f:.4f}')

# Show overconfidence via log-likelihoods
jll_1f = gnb_1feat._joint_log_likelihood(X_te_dup[:5, :1])
jll_2f = gnb_2feat._joint_log_likelihood(X_te_dup[:5])

print(f'\nLog-likelihood gap (|class0 - class1|) for first 5 test points:')
print(f'  1 feature:  {np.abs(jll_1f[:, 0] - jll_1f[:, 1])}')
print(f'  2 features: {np.abs(jll_2f[:, 0] - jll_2f[:, 1])}')
print(f'\nDuplicate features roughly DOUBLE the log-likelihood gap,')
print(f'making the model overconfident without improving accuracy.')

# Create data where feature 2 is a noisy copy of feature 1 n_fail = 300 X0_base = rng.normal(0, 1, size=(n_fail, 1)) X1_base = rng.normal(2, 1, size=(n_fail, 1)) # Duplicate feature with small noise (near-perfect copy) noise_level = 0.01 X0_dup = np.hstack([X0_base, X0_base + rng.normal(0, noise_level, size=(n_fail, 1))]) X1_dup = np.hstack([X1_base, X1_base + rng.normal(0, noise_level, size=(n_fail, 1))]) X_dup = np.vstack([X0_dup, X1_dup]) y_dup = np.array([0] * n_fail + [1] * n_fail) idx_dup = rng.permutation(len(y_dup)) X_dup, y_dup = X_dup[idx_dup], y_dup[idx_dup] X_tr_dup, X_te_dup, y_tr_dup, y_te_dup = train_test_split( X_dup, y_dup, test_size=0.3, random_state=SEED ) # Train on 1 feature vs 2 duplicate features gnb_1feat = GaussianNBScratch(var_smoothing=1e-9) gnb_1feat.fit(X_tr_dup[:, :1], y_tr_dup) gnb_2feat = GaussianNBScratch(var_smoothing=1e-9) gnb_2feat.fit(X_tr_dup, y_tr_dup) acc_1f = gnb_1feat.score(X_te_dup[:, :1], y_te_dup) acc_2f = gnb_2feat.score(X_te_dup, y_te_dup) print(f'Accuracy with 1 feature: {acc_1f:.4f}') print(f'Accuracy with 2 duplicate features: {acc_2f:.4f}') # Show overconfidence via log-likelihoods jll_1f = gnb_1feat._joint_log_likelihood(X_te_dup[:5, :1]) jll_2f = gnb_2feat._joint_log_likelihood(X_te_dup[:5]) print(f'\nLog-likelihood gap (|class0 - class1|) for first 5 test points:') print(f' 1 feature: {np.abs(jll_1f[:, 0] - jll_1f[:, 1])}') print(f' 2 features: {np.abs(jll_2f[:, 0] - jll_2f[:, 1])}') print(f'\nDuplicate features roughly DOUBLE the log-likelihood gap,') print(f'making the model overconfident without improving accuracy.')

Accuracy with 1 feature:           0.7889
Accuracy with 2 duplicate features: 0.7889

Log-likelihood gap (|class0 - class1|) for first 5 test points:
  1 feature:  [0.10608085 2.89029082 2.03506952 1.55886487 1.85526616]
  2 features: [0.21845117 5.80786212 4.03792952 3.13102919 3.75942679]

Duplicate features roughly DOUBLE the log-likelihood gap,
making the model overconfident without improving accuracy.

In [12]:

# Visualize: decision boundary for correlated data
from ml_first_principles.visualization import plot_decision_boundary

fig = plot_decision_boundary(
    gnb_corr, X_corr, y_corr,
    title='Gaussian NB — Correlated Features (rho=0.8)'
)
plt.show()

print('The decision boundary is still reasonable despite the independence')
print('assumption being violated — NB is robust for classification.')

# Visualize: decision boundary for correlated data from ml_first_principles.visualization import plot_decision_boundary fig = plot_decision_boundary( gnb_corr, X_corr, y_corr, title='Gaussian NB — Correlated Features (rho=0.8)' ) plt.show() print('The decision boundary is still reasonable despite the independence') print('assumption being violated — NB is robust for classification.')

The decision boundary is still reasonable despite the independence
assumption being violated — NB is robust for classification.

6.3 Failure Case — Non-Gaussian Features¶

When a feature's true distribution is multimodal within a class, Gaussian NB assigns wrong density values.

In [13]:

# Class 0: bimodal distribution (mixture of two Gaussians)
n_bi = 200
X0_bi_a = rng.normal(-2, 0.5, size=(n_bi // 2, 1))
X0_bi_b = rng.normal(2, 0.5, size=(n_bi // 2, 1))
X0_bi = np.vstack([X0_bi_a, X0_bi_b])
X0_bi = np.hstack([X0_bi, rng.normal(0, 1, size=(n_bi, 1))])  # add a normal feature

# Class 1: unimodal at 0
X1_bi = rng.normal(0, 0.8, size=(n_bi, 2))

X_bi = np.vstack([X0_bi, X1_bi])
y_bi = np.array([0] * n_bi + [1] * n_bi)
idx_bi = rng.permutation(len(y_bi))
X_bi, y_bi = X_bi[idx_bi], y_bi[idx_bi]

X_tr_bi, X_te_bi, y_tr_bi, y_te_bi = train_test_split(
    X_bi, y_bi, test_size=0.3, random_state=SEED
)

gnb_bi = GaussianNBScratch(var_smoothing=1e-9)
gnb_bi.fit(X_tr_bi, y_tr_bi)
acc_bi = gnb_bi.score(X_te_bi, y_te_bi)

print(f'Accuracy on bimodal data: {acc_bi:.4f}')
print(f'\nGaussian NB struggles here because class 0 has a bimodal')
print(f'feature 1 distribution — a single Gaussian cannot capture it.')

fig, ax = plt.subplots(1, 1, figsize=(8, 4))
ax.hist(X0_bi[:, 0], bins=30, alpha=0.5, label='Class 0 (bimodal)', density=True)
ax.hist(X1_bi[:, 0], bins=30, alpha=0.5, label='Class 1 (unimodal)', density=True)
ax.set(title='Feature 1 Distribution — Bimodal vs Unimodal',
       xlabel='Feature 1', ylabel='Density')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Class 0: bimodal distribution (mixture of two Gaussians) n_bi = 200 X0_bi_a = rng.normal(-2, 0.5, size=(n_bi // 2, 1)) X0_bi_b = rng.normal(2, 0.5, size=(n_bi // 2, 1)) X0_bi = np.vstack([X0_bi_a, X0_bi_b]) X0_bi = np.hstack([X0_bi, rng.normal(0, 1, size=(n_bi, 1))]) # add a normal feature # Class 1: unimodal at 0 X1_bi = rng.normal(0, 0.8, size=(n_bi, 2)) X_bi = np.vstack([X0_bi, X1_bi]) y_bi = np.array([0] * n_bi + [1] * n_bi) idx_bi = rng.permutation(len(y_bi)) X_bi, y_bi = X_bi[idx_bi], y_bi[idx_bi] X_tr_bi, X_te_bi, y_tr_bi, y_te_bi = train_test_split( X_bi, y_bi, test_size=0.3, random_state=SEED ) gnb_bi = GaussianNBScratch(var_smoothing=1e-9) gnb_bi.fit(X_tr_bi, y_tr_bi) acc_bi = gnb_bi.score(X_te_bi, y_te_bi) print(f'Accuracy on bimodal data: {acc_bi:.4f}') print(f'\nGaussian NB struggles here because class 0 has a bimodal') print(f'feature 1 distribution — a single Gaussian cannot capture it.') fig, ax = plt.subplots(1, 1, figsize=(8, 4)) ax.hist(X0_bi[:, 0], bins=30, alpha=0.5, label='Class 0 (bimodal)', density=True) ax.hist(X1_bi[:, 0], bins=30, alpha=0.5, label='Class 1 (unimodal)', density=True) ax.set(title='Feature 1 Distribution — Bimodal vs Unimodal', xlabel='Feature 1', ylabel='Density') ax.legend() ax.grid(True, alpha=0.3) plt.tight_layout() plt.show()

Accuracy on bimodal data: 0.9333

Gaussian NB struggles here because class 0 has a bimodal
feature 1 distribution — a single Gaussian cannot capture it.

7. Connections¶

From	To	Link
Naive Bayes (generative)	Logistic Regression (discriminative)	Same posterior form, different parameter estimation
Class-conditional Gaussians	Gaussian Mixture Models / Clustering	Unsupervised version drops labels
Naive independence	Feature engineering	Decorrelate features to better match the assumption
Laplace smoothing	Bayesian priors	Smoothing = Dirichlet prior on multinomial parameters

Takeaway¶

Naive Bayes is the simplest probabilistic classifier:

• No optimization loop — closed-form parameter estimation.
• Scales linearly with data size and feature count.
• Strong baseline for text classification and high-dimensional data.
• Probability estimates are poor (due to the independence assumption), but the classification ranking is often correct.

When you need better calibrated probabilities, use logistic regression (the discriminative counterpart) or calibrate NB's outputs post-hoc.