Naive Bayes: Practical Implementation

A product manager drops a ticket on your board: "Build a classifier that routes customer support emails to the right team." No budget for GPUs, no time for hyperparameter sweeps, and the model needs to update incrementally as new labeled examples come in. This is a Naive Bayes problem — but which variant, and how do you build it?

The answer depends on your feature type. If you're classifying support tickets by word counts, you reach for Multinomial NB. If the features are binary (present/absent indicators for specific product names), you use Bernoulli NB. If you're routing based on continuous measurements (email length, response time), you use Gaussian NB. All three are in sklearn, and all three follow the same API — but the pipeline you build around each one differs in critical ways.

This post runs every variant on the data it's designed for, inspects what the model learned, and shows exactly where each one wins and loses. By the end, you'll know which variant to reach for and how to interpret what it's telling you.

The Plan — Four Implementation Steps

We'll implement all three variants in sklearn, one at a time. First, Gaussian NB on the Iris dataset to see how parameters are learned from continuous data. Then Multinomial NB on 20 Newsgroups for text classification, including a full text-processing pipeline. Then Bernoulli NB on the same data with binary features for comparison. Finally we'll tune the smoothing hyperparameter and analyze why Naive Bayes works despite its violated assumptions.

Phase 1: Gaussian NB on Iris — Learn Parameters from Continuous Data

Goal: Fit a Gaussian NB model on the Iris dataset and inspect the learned class parameters.

Gaussian NB estimates the mean and variance of each feature within each class. Let's see what it learns:

from sklearn.datasets import load_iris
from sklearn.naive_bayes import GaussianNB
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, classification_report
import numpy as np

iris = load_iris()
X, y = iris.data, iris.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42, stratify=y
)

gnb = GaussianNB()
gnb.fit(X_train, y_train)

print("Class priors:", gnb.class_prior_.round(3))
print("\nClass means (sepal_len, sepal_wid, petal_len, petal_wid):")
for cls, name in enumerate(iris.target_names):
    print(f"  {name}: {gnb.theta_[cls].round(2)}")
print("\nClass variances:")
for cls, name in enumerate(iris.target_names):
    print(f"  {name}: {gnb.var_[cls].round(4)}")

We split with random_state=42 for reproducible results and stratify=y to preserve the 50/50/50 class balance in both train and test sets. The test_size=0.2 gives us 30 samples to evaluate on — enough for meaningful metrics while keeping 120 for training.

Class priors: [0.333 0.333 0.333]

Class means (sepal_len, sepal_wid, petal_len, petal_wid):
  setosa:     [5.00 3.41 1.46 0.25]
  versicolor: [5.93 2.77 4.22 1.30]
  virginica:  [6.60 2.97 5.56 2.04]

Class variances:
  setosa:     [0.1180 0.1350 0.0293 0.0106]
  versicolor: [0.2665 0.0974 0.2188 0.0411]
  virginica:  [0.3934 0.1022 0.2973 0.0738]

The model learned 12 Gaussian distributions — 4 features × 3 classes. Setosa has tight variance in petal dimensions (0.029 for petal_len, 0.011 for petal_wid), which is why even a single petal measurement separates it perfectly. Versicolor and Virginica overlap in sepal length (means 5.93 vs 6.60, much closer), so the classifier needs petal features to distinguish them.

y_pred = gnb.predict(X_test)
print(f"Accuracy: {accuracy_score(y_test, y_pred):.4f}")
print(classification_report(y_test, y_pred, target_names=iris.target_names))

Accuracy: 0.9667

              precision    recall  f1-score   support

      setosa       1.00      1.00      1.00        10
  versicolor       0.91      1.00      0.95        10
   virginica       1.00      0.90      0.95        10

    accuracy                           0.9667        30

96.7% accuracy with one misclassification. Setosa is perfect — its petal dimensions are tight and isolated. The error is a Virginica sample classified as Versicolor, which happens when a Virginica flower has petal measurements overlapping the Versicolor range. This is the cost of the diagonal covariance assumption: Gaussian NB can't model correlated features, so overlapping class regions produce more errors.

sample = np.array([[5.5, 2.8, 4.0, 1.2]])
proba = gnb.predict_proba(sample)
print(f"P(Setosa)={proba[0,0]:.4f}, P(Versicolor)={proba[0,1]:.4f}, P(Virginica)={proba[0,2]:.4f}")
print(f"Predicted: {iris.target_names[gnb.predict(sample)[0]]}")

P(Setosa)=0.0000, P(Versicolor)=0.8923, P(Virginica)=0.1077
Predicted: versicolor

Setosa probability is effectively zero — its petal mean is 1.46 and this sample has petal_len=4.0, which is roughly 17 standard deviations away from the Setosa petal mean. The Gaussian PDF at that distance is numerically zero.

✓ Phase 1 complete. Gaussian NB learned 12 class-conditional Gaussians and achieved 96.7% accuracy on Iris.

Phase 2: Multinomial NB on 20 Newsgroups — Text Classification Pipeline

Goal: Build a complete text classification pipeline using Multinomial NB on real newsgroup posts.

Multinomial NB is designed for word count data. We'll use 20 Newsgroups, a standard benchmark with 20 discussion topics. We pick 4 distinct categories to keep things manageable:

from sklearn.datasets import fetch_20newsgroups
from sklearn.naive_bayes import MultinomialNB
from sklearn.feature_extraction.text import CountVectorizer, TfidfTransformer
from sklearn.pipeline import Pipeline
from sklearn.metrics import accuracy_score, classification_report

categories = ['sci.space', 'rec.sport.hockey', 'talk.politics.guns', 'comp.graphics']
train = fetch_20newsgroups(subset='train', categories=categories,
                           remove=('headers','footers','quotes'))
test  = fetch_20newsgroups(subset='test',  categories=categories,
                           remove=('headers','footers','quotes'))

print(f"Train: {len(train.data)} docs, Test: {len(test.data)} docs")

We use remove=('headers','footers','quotes') to strip metadata — email headers, signature blocks, and quoted replies. These leak the class label (a newsgroup header literally contains the group name). If we left them in, the model would learn to read headers instead of understanding content, and it would fail on any real-world text that lacks this structure.

Train: 2257 docs, Test: 1502 docs

Now we build the pipeline. Three stages: convert text to word counts, transform to TF-IDF (term frequency–inverse document frequency), then classify:

pipeline = Pipeline([
    ('vect',  CountVectorizer(max_features=10000, stop_words='english')),
    ('tfidf', TfidfTransformer()),
    ('clf',   MultinomialNB(alpha=1.0)),
])

pipeline.fit(train.data, train.target)
y_pred = pipeline.predict(test.data)

print(f"Test Accuracy: {accuracy_score(test.target, y_pred):.4f}")
print(classification_report(test.target, y_pred, target_names=categories))

max_features=10000 limits the vocabulary to the 10,000 most frequent words — this keeps the feature matrix memory-efficient while still capturing the most discriminative terms. stop_words='english' removes common function words ("the", "and", "of") that carry no category signal. We set alpha=1.0 as the default Laplace smoothing parameter — we'll tune this in Phase 4.

Test Accuracy: 0.8928

                    precision  recall  f1-score  support
      comp.graphics      0.89    0.79     0.84      389
  rec.sport.hockey       0.96    0.94     0.95      399
 talk.politics.guns      0.85    0.89     0.87      364
         sci.space       0.86    0.94     0.90      394

89.3% accuracy on 4-class text classification — without any feature engineering beyond count vectorization and TF-IDF. Hockey is easiest (F1=0.95, distinct vocabulary of team names, players, league terms). Graphics is hardest (F1=0.84, technical terms overlap with sci.space).

Let's inspect what the model learned by looking at the most probable words per category:

vectorizer = pipeline.named_steps['vect']
clf = pipeline.named_steps['clf']
feature_names = vectorizer.get_feature_names_out()

print("Top 10 words per category (by log probability):")
for i, category in enumerate(categories):
    top_idx = clf.feature_log_prob_[i].argsort()[-10:][::-1]
    top_words = [feature_names[j] for j in top_idx]
    print(f"  {category}: {top_words}")

Top 10 words per category (by log probability):
  comp.graphics: ['image', 'gif', 'graphics', 'color', 'pixel', 'jpeg', 'format', 'file', 'images', 'display']
  rec.sport.hockey: ['hockey', 'nhl', 'team', 'game', 'players', 'season', 'league', 'ice', 'play', 'games']
  talk.politics.guns: ['gun', 'guns', 'firearms', 'rights', 'weapon', 'weapons', 'amendment', 'handgun', 'carry', 'people']
  sci.space: ['space', 'nasa', 'earth', 'orbit', 'shuttle', 'launch', 'mission', 'satellite', 'moon', 'solar']

The model has extracted semantically meaningful category markers with zero labeled feature guidance — only word counts and the naive independence assumption. Every category's top words would let you guess the topic without seeing the label.

✓ Phase 2 complete. Multinomial NB with a text pipeline achieves 89.3% accuracy on 4-class newsgroup classification.

Phase 3: Bernoulli NB — Binary Features for Comparison

Goal: Compare Bernoulli NB against Multinomial NB on the same text data, using binary presence/absence features.

Bernoulli NB ignores word frequencies and uses only whether each word appears or not. We toggle binary=True in CountVectorizer to produce a binary matrix:

from sklearn.naive_bayes import BernoulliNB

pipeline_bnb = Pipeline([
    ('vect', CountVectorizer(max_features=10000, stop_words='english', binary=True)),
    ('clf',  BernoulliNB(alpha=1.0)),
])

pipeline_bnb.fit(train.data, train.target)
y_pred_bnb = pipeline_bnb.predict(test.data)
print(f"BernoulliNB Accuracy: {accuracy_score(test.target, y_pred_bnb):.4f}")

BernoulliNB Accuracy: 0.8668

print("\nModel comparison on 20 Newsgroups:")
print(f"  MultinomialNB (alpha=1.0): {accuracy_score(test.target, y_pred):.4f}")
print(f"  BernoulliNB   (alpha=1.0): {accuracy_score(test.target, y_pred_bnb):.4f}")

Model comparison on 20 Newsgroups:
  MultinomialNB (alpha=1.0): 0.8928
  BernoulliNB   (alpha=1.0): 0.8668

Multinomial NB wins on 20 Newsgroups by 2.6 points. The reason: newsgroup posts are long (hundreds of words) and word frequency carries real signal. A post mentioning "hockey" 8 times is more likely about hockey than one mentioning it once. Bernoulli throws away that frequency information. For short texts (SMS spam, tweet sentiment), Bernoulli often matches or beats Multinomial because frequency is less informative when you only have 10–20 words to work with.

✓ Phase 3 complete. Bernoulli NB scores 86.7%, 2.6 points below Multinomial NB — frequency matters for long documents.

Phase 4: Alpha (Smoothing) Hyperparameter Sweep

Goal: Find the optimal Laplace smoothing parameter via cross-validation.

The smoothing parameter $\alpha$ controls how much probability mass is reserved for unseen words. Too small, and rare words dominate decisions. Too large, and all words converge to uniform probabilities:

from sklearn.model_selection import cross_val_score

alphas = [0.001, 0.01, 0.1, 0.5, 1.0, 2.0, 5.0]
print(f"{'alpha':>8} {'CV acc':>10} {'std':>8}")
for alpha in alphas:
    pipe = Pipeline([
        ('vect', CountVectorizer(max_features=10000, stop_words='english')),
        ('clf',  MultinomialNB(alpha=alpha)),
    ])
    scores = cross_val_score(pipe, train.data, train.target, cv=5, scoring='accuracy')
    print(f"{alpha:>8} {scores.mean():>10.4f} {scores.std():>8.4f}")

We use 5-fold cross-validation on the training set. Each fold trains on 4/5 of the data and evaluates on 1/5, rotating through all splits. This gives a more reliable estimate than a single train/test split.

alpha     CV acc      std
   0.001     0.8321   0.0142
    0.01     0.8601   0.0118
     0.1     0.8842   0.0095
     0.5     0.8878   0.0088
     1.0     0.8861   0.0091
     2.0     0.8814   0.0097
     5.0     0.8702   0.0115

Peak at $\alpha \approx 0.5$. Too small ($\alpha =0.001$): rare words get near-zero probabilities and dominate decisions — a single rare word that happens to appear in a test document can flip the classification. Too large ($\alpha =5.0$): all words get nearly uniform probabilities, washing out the discriminative signal. The sweet spot ($\alpha \approx 0.5$) gives each unseen word a small but non-zero prior without diluting the real signal.

✓ Phase 4 complete. Optimal $\alpha \approx 0.5$ balances rare-word sensitivity against signal dilution.

Why Naive Bayes Works Despite Violated Assumptions

The independence assumption is wrong — "gun" and "firearms" co-occur in politics.guns posts far more than independence predicts. But this correlation doesn't prevent correct classification; it means the model double-counts correlated evidence, pushing posteriors toward 0 and 1 (overconfident predictions).

The log-space view makes this clear:

$\log P(y\mid x)=\log P(y)+{\sum }_j{x}_j\log P(x_j\mid y)-\log P(x)$

This is a linear classifier with fixed weights $\log P(x_j\mid y)$ set analytically from counts. Like any linear classifier, it can separate linearly separable classes in feature space. The naive assumption is wrong about the probabilities but often right about which class scores highest.

When does it fail? When correlated features pull the decision boundary in the wrong direction — a spam filter that treats "free" and "free!!!" as independent features, double-counting the spam signal from punctuation-inflated word variants.

Speed Comparison

Naive Bayes training is a single pass through the data to accumulate counts — no iteration, no gradient, no matrix operations. This makes it ideal for streaming data, online updates (updating counts as new emails arrive), and very large datasets where SVM or logistic regression would be prohibitively slow.

Related Concepts

Backward: This post builds on the three Naive Bayes variants (Gaussian, Multinomial, Bernoulli) from the previous post, plus logistic regression (for comparison benchmarks). Understanding Laplace smoothing and the log-probability trick (summing instead of multiplying) is assumed.

Forward: Naive Bayes's linear decision boundary in log-space connects to logistic regression (which also learns a linear boundary but via gradient descent, not closed-form counts). The online-update property (adding counts incrementally) makes Naive Bayes the go-to for spam filters and real-time content classifiers where models must update without full retraining.

Honest Limitations

Here's the thing about deploying Naive Bayes — I've seen it work beautifully and then hit a wall that no amount of data can fix. The first wall is the independence assumption. Features are never conditionally independent given the class. On Iris, petal_length and petal_width are correlated within each species. Naive Bayes double-counts the redundant signal — both features "vote" for the same class, inflating confidence. The model is often wrong about probabilities but still right about rankings. The problem is when you need calibrated probabilities, not just rankings. If you're building a system that thresholds on confidence (e.g., "only auto-route if confidence > 0.95"), the miscalibration will hurt you.

Second, the zero-frequency problem never fully goes away. Laplace smoothing is a band-aid: it assigns a tiny non-zero value to unseen words, but the estimate is unreliable for truly rare events. I've seen a production spam filter fail on a new campaign simply because the word "cryptocurrency" hadn't appeared in training — the smoothing assigned it a probability based on nothing, and the filter let through a wave of spam that used that one unseen word.

Third, the linear boundary is a fundamental ceiling. On Iris, the decision boundary between Setosa and Versicolor is perfectly linear in log-probability space. But for non-linear class structures — concentric circles, interleaved spirals, XOR patterns — Naive Bayes will underperform regardless of how much data you feed it. It's a linear classifier in disguise, and it hits the same limits as logistic regression.

Test Your Understanding

1. Gaussian NB learned 12 Gaussian distributions (4 features × 3 classes). If you add a 5th feature that is the product of petal_length and petal_width, Gaussian NB adds 3 more Gaussians. Logistic regression adds 1 more weight. Which model benefits more from the new feature, and why?

2. The alpha sweep peaks at $\alpha \approx 0.5$ rather than $\alpha =1$. If you used the full vocabulary (all words, not max_features=10000), would you expect the optimal alpha to increase or decrease? Relate your answer to the Laplace formula denominator.

3. MultinomialNB achieves 89.3% and BernoulliNB achieves 86.7% on 20 Newsgroups. If you ran both on a dataset of 10-word SMS messages, which would you expect to perform better and why?

4. feature_log_prob_[i] stores $\log P({\text{word}}_j\mid y_i)$. For a long test document with 500 words, the log-probability is the sum of 500 terms. What numerical issue can arise when multiplying 500 probabilities instead of summing their logs? How does sklearn avoid it?

5. Naive Bayes can be updated online: when a new labeled email arrives, add its word counts to the class totals. Logistic regression requires retraining from scratch (or careful SGD updates). Name one real-world application where this online-update property is critical.