32 - Naive Bayes: Bayes' theorem with a naive assumption

A doctor who receives a patient with fever, cough and a sore throat does not reason about one symptom at a time: she starts from how frequent each disease is in general, then updates her conviction at every clue. A spam filter does the same with the words of a message: "free" shifts the judgement towards spam, "meeting" towards a legitimate message, and the final decision is the sum of all the clues. Naive Bayes is this: a classifier that applies Bayes' theorem and, to make the calculations feasible, makes a deliberately simplified assumption (hence "naive") about the relationship between the clues.

It is one of the simplest algorithms in machine learning and for this very reason an excellent starting point: you will implement it by hand in a few lines, you will see that it coincides with the library, and you will understand in which cases that naive assumption is paid for dearly.

An intuition before the formulas

Imagine receiving a message that contains the word "free". How convinced are you that it is spam? It depends on two things. First, on how much spam arrives in general in your inbox (if 90% of the mail is spam, you start out already suspicious). Second, on how typical the word "free" is of spam compared with legitimate messages. Bayes' theorem is the rule that combines these two ingredients coherently. A numerical example, with figures invented only to illustrate the mechanism: 10% of the messages are spam; a filter that looks for "free" finds it in 95% of the spam and in 2% of the legitimate messages. If the filter fires, how likely is it to be spam?

spam and the filter fires:       0.10 × 0.95 = 0.095
legitimate and the filter fires: 0.90 × 0.02 = 0.018
P(spam | fires) = 0.095 / (0.095 + 0.018) = 0.84

The filter's 95% "sensitivity" does not mean a 95% probability of spam: the final probability is 84%, because legitimate messages are many and even a 2% of false alarms weighs. Getting this step wrong (confusing P(fires | spam) with P(spam | fires)) is one of the most common reasoning errors, and Bayes is exactly the formula that connects them.

Bayes' theorem, derived in three lines

The conditional probability of A given B is the fraction of the cases in which B happens in which A also happens:

P(A | B) = P(A ∩ B) / P(B)

The same joint event can be written in two ways, P(A ∩ B) = P(B | A) · P(A). Substituting into the previous formula gives Bayes' theorem:

P(A | B) = P(B | A) · P(A) / P(B)

Applied to classification, A is the class c (for example "spam") and B are the observed features x = (x₁, ..., x_n). The three quantities have precise names:

P(c | x)  =  P(x | c) · P(c)  /  P(x)
 posterior    likelihood · prior / evidence

• Prior P(c): how frequent the class is before seeing the message (the "base rate" of spam).
• Likelihood P(x | c): how typical it would be to see these features if the class were c.
• Posterior P(c | x): what we want, the probability of the class after seeing the data.

The three quantities have different roles and only two are estimated from the data: the prior is counted on the labels, the likelihood is modelled with the naive assumption, the posterior is what gets computed. The most frequent confusion is swapping likelihood and posterior, that is P(fires | spam) with P(spam | fires), and it is the error that the numerical example of the filter has just shown.

A notation before going on: argmax. Writing argmax_c f(c) means "the value of c for which f(c) is maximum", not the maximum itself. Here it is the class with the highest score.

To choose the class there is no need to compute P(x): it is the same for all the classes. It is enough to compare the numerators, so the decision rule is

ŷ = argmax_c  P(c) · P(x₁, ..., x_n | c)

The naive assumption

The problem is P(x₁, ..., x_n | c): estimating the joint probability of all the features together would require, with n binary features, about 2ⁿ parameters per class, and data you will never have. The "naive" assumption is to assume that, knowing the class, the features are conditionally independent:

P(x₁, ..., x_n | c) = P(x₁ | c) · P(x₂ | c) · ... · P(x_n | c)

From 2ⁿ parameters we go to about n per class. The final rule is:

ŷ = argmax_c  P(c) · Π_i P(x_i | c)

The assumption is almost always false: in a text, "prize" and "win" often appear together, in spam too. Yet the classifier works surprisingly well, because to classify the probabilities do not need to be right, it is enough that the winning class is the right one. When the assumption is too wrong, however, the resulting probabilities become unreliable: we will see it with an experiment.

Two practical tricks: logarithms and smoothing

Logarithms. Multiplying a hundred probabilities smaller than 1 produces numbers so small that the computer rounds them to zero (underflow). Since the logarithm is increasing, the argmax does not change if you work on the logarithms, where products become sums:

ŷ = argmax_c  [ log P(c) + Σ_i log P(x_i | c) ]

Laplace smoothing. If a word has never appeared in the training spam, P(word | spam) = 0, and that single zero wipes out the whole product: one unknown term cancels the weight of all the other words. It is fixed by "pretending" to have seen every word α more times (usually α = 1): seen up close it is a prior, and the lesson on Bayesian statistics derives it as the mean of a Dirichlet posterior. For the multinomial model (word counts):

P(w | c) = ( count(w, c) + α ) / ( N_c + α · |V| )

where count(w, c) is how many times w appears in the messages of class c, N_c is the total number of words in that class and |V| the size of the vocabulary. The denominator adds α·|V| because we have added α to each of the |V| words: this way the probabilities still sum to 1.

A hand implementation on a mini spam filter

We use ten messages invented for the occasion, small enough to follow the numbers.

import numpy as np

spam = ["win a prize free now", "offer free click now",
        "prize win money immediately", "click here for winning money free",
        "offer exclusive prize free"]
ham  = ["meeting tomorrow at ten", "sending you the report of meeting",
        "lunch tomorrow with the team", "the report is ready for tomorrow",
        "we meet to lunch after that meeting"]

docs = spam + ham
y = np.array([1] * len(spam) + [0] * len(ham))        # 1 = spam
vocab = sorted({w for d in docs for w in d.split()})
idx = {w: i for i, w in enumerate(vocab)}

def conta(testo):
    v = np.zeros(len(vocab))
    for w in testo.split():
        if w in idx:                 # out-of-vocabulary words: ignored
            v[idx[w]] += 1
    return v

X = np.array([conta(d) for d in docs])
print("documents x words:", X.shape)

alpha = 1.0
prior = np.array([(y == c).mean() for c in (0, 1)])           # P(c)
cont = np.array([X[y == c].sum(axis=0) for c in (0, 1)])      # count(w, c)
verosim = (cont + alpha) / (cont.sum(axis=1, keepdims=True) + alpha * len(vocab))

def log_punteggi(testo):
    x = conta(testo)
    return np.log(prior) + x @ np.log(verosim).T      # log P(c) + sum_i x_i log P(w_i|c)

def posteriori(testo):
    s = log_punteggi(testo)
    e = np.exp(s - s.max())                           # stable softmax
    return e / e.sum()

for t in ["free prize tomorrow", "report of meeting", "click for winning money"]:
    p = posteriori(t)
    print("%-28s P(ham)=%.3f  P(spam)=%.3f" % (t, p[0], p[1]))

documents x words: (10, 32)
free prize tomorrow          P(ham)=0.133  P(spam)=0.867
report of meeting            P(ham)=0.949  P(spam)=0.051
click for winning money      P(ham)=0.038  P(spam)=0.962

What the code does, piece by piece. Vocabulary: the sorted list of all the words present in the training messages (32 words); idx assigns each one a column. conta turns a message into a vector of 32 numbers, one per word, with how many times it appears: it is the "bag of words" representation, which throws away the order and keeps only the counts. X is the matrix (10 messages × 32 words). prior holds P(c) for the two classes, both 0.5. cont sums the word counts separately per class, and verosim applies the formula with smoothing: the sum of each row is 1. log_punteggi computes log P(c) + Σ x_i log P(w_i|c) for the two classes; posteriori converts them back into normalised probabilities.

Note the line x @ np.log(verosim).T: multiplying by the vector of counts is equivalent to adding log P(w|c) once for each occurrence of the word. The final conversion from logarithmic scores to probabilities uses the trick of subtracting the maximum before exponentiating, to avoid overflow.

Take "free prize tomorrow": two words typical of spam and one typical of legitimate messages. The model assigns 0.867 to spam: not a simple 2 against 1 count, but the result of the learned frequencies, of the prior (5 messages out of 10 are spam) and of the smoothing.

Now the library, and a check that it coincides

from sklearn.naive_bayes import MultinomialNB

nb = MultinomialNB(alpha=1.0).fit(X, y)
for t in ["free prize tomorrow", "report of meeting", "click for winning money"]:
    print("%-28s sklearn P(spam)=%.3f   by hand P(spam)=%.3f" % (
        t, nb.predict_proba(conta(t).reshape(1, -1))[0, 1], posteriori(t)[1]))

free prize tomorrow          sklearn P(spam)=0.867   by hand P(spam)=0.867
report of meeting            sklearn P(spam)=0.051   by hand P(spam)=0.051
click for winning money      sklearn P(spam)=0.962   by hand P(spam)=0.962

What the code does: MultinomialNB(alpha=1.0) is the same thing we wrote, with the same smoothing. fit(X, y) learns priors and likelihoods; predict_proba returns the probabilities for the two classes, and we take column 1 (spam). The reshape(1, -1) turns the vector into a one-row matrix, because the library always expects a matrix (n samples × d features).

The numbers coincide: the library does nothing more than what you have written, and this is worth remembering. When a model looks like a black box to you, it is often a simple idea with careful engineering.

word              clogged_extruder     unlevelled_bed
under_extrusion           6                     1
click                     4                     0
first_layer               1                     7
detachment                0                     5
temperature               3                     1
TOTAL                    14                    14

Continuous features: Gaussian Naive Bayes

With numerical features (heights, measurements) you do not count words: you assume that within each class each feature follows a Gaussian, whose mean μ and variance σ² you estimate from the data of that class:

P(x_i | c) = 1 / sqrt(2π σ²_ic) · exp( −(x_i − μ_ic)² / (2 σ²_ic) )

log P(x_i | c) = −½ [ log(2π σ²_ic) + (x_i − μ_ic)² / σ²_ic ]

Let us implement it in vectorised form on the Iris dataset (three species of flowers, four measurements), with a stratified split as in the previous lesson. The code writes the class from scratch in vectorised form: for each class it keeps the prior, the mean and the variance of each measurement, and evaluates the logarithm of the Gaussian on all the samples and all the classes in one go, then compares it with the scikit-learn one on the same predictions:

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

class GaussianNBDaZero:
    def fit(self, X, y):
        self.classi = np.unique(y)
        self.prior = np.array([(y == c).mean() for c in self.classi])
        self.mu = np.array([X[y == c].mean(axis=0) for c in self.classi])
        self.var = np.array([X[y == c].var(axis=0) for c in self.classi]) + 1e-9
        return self

    def log_congiunta(self, X):
        # X: (n, d)  mu, var: (k, d)  -> result (n, k)
        diff2 = (X[:, None, :] - self.mu[None, :, :]) ** 2
        ll = -0.5 * (np.log(2 * np.pi * self.var)[None] + diff2 / self.var[None]).sum(axis=2)
        return ll + np.log(self.prior)

    def predict(self, X):
        return self.classi[self.log_congiunta(X).argmax(axis=1)]

Xi, yi = load_iris(return_X_y=True)
Xi_tr, Xi_te, yi_tr, yi_te = train_test_split(Xi, yi, test_size=0.3, stratify=yi, random_state=0)

mio = GaussianNBDaZero().fit(Xi_tr, yi_tr)
lib = GaussianNB().fit(Xi_tr, yi_tr)
print("accuracy by hand:    %.3f" % (mio.predict(Xi_te) == yi_te).mean())
print("accuracy sklearn:    %.3f" % lib.score(Xi_te, yi_te))
print("identical predictions:", bool((mio.predict(Xi_te) == lib.predict(Xi_te)).all()))
print("shape of mu:", mio.mu.shape, "-> (classes, features)")

accuracy by hand:    0.978
accuracy sklearn:    0.978
identical predictions: True
shape of mu: (3, 4) -> (classes, features)

What the code does: fit computes, for each class, the prior, the mean of each feature and the variance of each feature (the + 1e-9 avoids divisions by zero when a variance is null). log_congiunta evaluates the formula of the logarithm of the Gaussian for each sample, class and feature, and sums them over the features. The (n, k) result holds one score per sample-class pair and predict takes the argmax per row. Iris has 3 classes and 4 features, so mu has shape (3, 4). The fact that the two implementations make the same 45 predictions out of 45 is the confirmation that the library hides nothing more.

The shape of the tensors is the part that is most often gotten wrong: X[:, None, :] has shape (n, 1, d), self.mu[None, :, :] has shape (1, k, d), broadcasting produces (n, k, d) and the sum over the last axis gives (n, k): one score for each (sample, class) pair.

The trap: correlated features and overconfident probabilities

When two features are almost the same information, the model counts it twice, and the evidence is inflated. To see it in an extreme way, we repeatedly duplicate a column of Iris and observe how the probability assigned to an ambiguous point changes:

punto = Xi_te[[np.argmax(lib.predict_proba(Xi_te).max(axis=1) < 0.9)]]   # a non-trivial point
print("probabilities with the 4 original features:", np.round(lib.predict_proba(punto)[0], 3))

for copie in (3, 10):
    Xd = np.hstack([Xi_tr] + [Xi_tr[:, [2]]] * copie)     # repeat column 2 (petal length)
    pd_ = np.hstack([punto] + [punto[:, [2]]] * copie)
    g = GaussianNB().fit(Xd, yi_tr)
    print("with %2d extra copies of column 2:" % copie, np.round(g.predict_proba(pd_)[0], 3))

probabilities with the 4 original features: [0.   0.78 0.22]
with  3 extra copies of column 2: [0.    0.981 0.019]
with 10 extra copies of column 2: [0. 1. 0.]

What the code does. The first line picks a "non-trivial" point of the test set: predict_proba(...).max(axis=1) < 0.9 is a vector of true/false (the model is not very sure), argmax takes the first true; the double square brackets keep the matrix shape. In the loop, np.hstack places next to the 4 original columns copie copies of column 2 (petal length), both in the training set and on the point; then it trains a new GaussianNB and prints the point's probabilities for the three species.

The point starts at (0, 0.78, 0.22). With three extra copies of the same column it becomes (0, 0.981, 0.019); with ten, (0, 1, 0). The duplicated features bring no new information, but the model treats them as independent clues and pushes the probabilities towards the extremes. The accuracy often does not change much (the winning class can stay the same), but the probabilities produced are no longer calibrated: a "99%" does not mean that in 99% of similar cases the prediction is right. If you need reliable probabilities (for example to rank by risk or to decide a threshold) you must calibrate them separately, or use another model.

The three variants, in brief

• MultinomialNB: features = counts (words, frequencies). The classic for texts.
• BernoulliNB: binary features (the word is there / is not there). It also penalises the absence of a word, unlike the multinomial.
• GaussianNB: continuous features, one Gaussian per feature and class.

The three variants are chosen according to the type of feature, not preference: counts for the multinomial, presence or absence for the Bernoulli, continuous measurements for the Gaussian. Choosing the wrong variant gives an error (negative values in the multinomial) or, worse, a model that runs but ignores the information that matters, as BernoulliNB shows on long texts.

Trade-offs

• Pros: immediate training (just count), little data is enough, it scales to thousands of features, easy to interpret and to update.
• Cons: the conditional independence assumption is rarely true; the probabilities are poorly calibrated; it does not capture interactions between features (the joint presence of two words means something different from the sum).
• When to use it: a quick baseline for text classification, little data, need for a very light model. When not: highly correlated features and probabilities that must be reliable.

The balance can be read in only one way: Naive Bayes is an excellent starting point and a poor estimator of probabilities. It is adopted when speed, simplicity and little data matter, it is abandoned when the features are highly correlated or when the probabilities must be reliable. The most important line is the last, because it says when to stop using it.

Real traps

• Data leakage in the vocabulary: if you build the vocabulary (or CountVectorizer/TfidfVectorizer) on the whole dataset before splitting, the test set has already influenced the training. The vectoriser goes inside a Pipeline and is fitted only on the train set.
• Zero without smoothing: with α = 0 a word never seen in a class zeroes the probability. Keep α > 0.
• Imbalanced classes: the prior P(c) weighs a lot. With 1% spam, the model tends to say "ham" until the evidence is overwhelming. Look at precision and recall, not accuracy (previous lesson).
• Shape mismatch: predict_proba wants a matrix (n, d); passing a single vector (d,) gives an error. For one sample you need reshape(1, -1).
• Wrong features for the wrong variant: giving negative values (for example standardised features) to MultinomialNB raises an error, because counts cannot be negative.

The five traps are all preparation errors, not algorithm errors: the vocabulary built on the whole dataset, the zero without smoothing, the ignored prior, the wrong shape of the input, the wrong variant for the type of feature. Four give an explicit error or show up in the numbers, the first does not: it is the one to check by hand with a pipeline.

A test on real texts: 20 newsgroups

The mini spam filter with ten messages was didactic, but it says nothing about the behaviour of Naive Bayes on real texts. We use the classic "20 newsgroups" corpus: messages posted to discussion groups, each labelled with its group. We take four topics (sci.space, rec.autos, talk.politics.guns, comp.graphics), 2,317 messages for training and 1,543 for the test. On the first run scikit-learn downloads the corpus (a few seconds). The message contains, besides the text, headers, footers and quotes of previous messages: information that in reality you might not have, or that identifies the group without the model having to understand the text. Let us measure how much they weigh.

from sklearn.datasets import fetch_20newsgroups
from sklearn.feature_extraction.text import CountVectorizer, TfidfVectorizer
from sklearn.naive_bayes import MultinomialNB, BernoulliNB, ComplementNB
from sklearn.pipeline import make_pipeline
from sklearn.linear_model import LogisticRegression

categories = ["sci.space", "rec.autos", "talk.politics.guns", "comp.graphics"]

def load(subset, clean):
    extra = dict(remove=("headers", "footers", "quotes")) if clean else {}
    return fetch_20newsgroups(subset=subset, categories=categories, random_state=0, **extra)

for clean in (False, True):
    train, test = load("train", clean), load("test", clean)
    pipe = make_pipeline(CountVectorizer(min_df=2), MultinomialNB(alpha=1.0)).fit(train.data, train.target)
    print("metadata removed: %-5s   training examples %d, test %d   vocabulary %d words   test accuracy %.3f" % (
        clean, len(train.data), len(test.data), len(pipe[0].vocabulary_), pipe.score(test.data, test.target)))

metadata removed: False   training examples 2317, test 1543   vocabulary 19001 words   test accuracy 0.962
metadata removed: True    training examples 2317, test 1543   vocabulary 13134 words   test accuracy 0.866

What the code does. CountVectorizer turns each text into a vector of word counts (the "bag of words" representation of the lesson), discarding the words that appear in fewer than two documents (min_df=2). The vectoriser and the classifier sit in the same Pipeline: the vocabulary is learned only on the training set, and the test set is transformed with that vocabulary (no leakage from the vocabulary). The difference between the two lines is the instructive part: with headers, footers and quotes the test accuracy is 0.962; removing them it drops to 0.866. Almost ten points are not "ability to understand the text": they are metadata (e-mail addresses, organisation names, "Subject:" lines, group names in the signatures) that identify the group without any need to read. A model evaluated like this would give an inflated accuracy, and in production, where that metadata is not there, it would collapse: it is leakage in its most concrete form. From here on we use the clean version, which is the honest estimate.

Which variant, which smoothing?

train, test = load("train", True), load("test", True)

print("smoothing alpha  ->  accuracy (Multinomial Naive Bayes)")
for alpha in (0.001, 0.01, 0.1, 0.5, 1.0, 3.0):
    p = make_pipeline(CountVectorizer(min_df=2), MultinomialNB(alpha=alpha)).fit(train.data, train.target)
    print("   %-6g  ->  %.3f" % (alpha, p.score(test.data, test.target)))

print("\nvariants and comparison with a discriminative model")
candidates = {
    "MultinomialNB (counts)": make_pipeline(CountVectorizer(min_df=2), MultinomialNB(alpha=0.1)),
    "BernoulliNB (presence/absence)": make_pipeline(CountVectorizer(min_df=2, binary=True), BernoulliNB()),
    "ComplementNB": make_pipeline(CountVectorizer(min_df=2), ComplementNB()),
    "logistic regression on TF-IDF": make_pipeline(TfidfVectorizer(min_df=2), LogisticRegression(max_iter=2000)),
}
for name, p in candidates.items():
    p.fit(train.data, train.target)
    print("%-48s  accuracy %.3f" % (name, p.score(test.data, test.target)))

smoothing alpha  ->  accuracy (Multinomial Naive Bayes)
   0.001   ->  0.865
   0.01    ->  0.867
   0.1     ->  0.865
   0.5     ->  0.866
   1       ->  0.866
   3       ->  0.854

variants and comparison with a discriminative model
MultinomialNB (counts)                            accuracy 0.865
BernoulliNB (presence/absence)                    accuracy 0.684
ComplementNB                                      accuracy 0.870
logistic regression on TF-IDF                     accuracy 0.839

What the code does. The first loop varies the smoothing α over five orders of magnitude: the accuracy stays between 0.865 and 0.867 for α from 0.001 to 1 and drops to 0.854 only with α = 3. Here the value of α matters little, because the corpus is large and the vocabulary very rich; the effect of smoothing is felt with little data or with rare words.

The second comparison is more interesting. BernoulliNB, which models the presence or absence of each word, gets only 0.684: in long texts almost all the words of the vocabulary are absent, and the model, which also penalises absence, is dominated by thousands of uninformative "not there"s. ComplementNB (a variant designed for texts and imbalanced classes, which estimates the parameters using the documents of all the other classes) scores 0.870, just above the multinomial (0.865). Logistic regression on TF-IDF scores 0.839: lower, in this case.

No crushing victory for one approach: on short texts with a lot of vocabulary and few examples, Naive Bayes is a baseline that is hard to beat and it trains in a fraction of a second. Differences of half a point (0.865 against 0.870) should not be commented on: they are below the sampling noise of a test of 1,543 messages (about ±0.9 points).

What did the model learn? The most indicative words

An advantage of Naive Bayes is readability: the parameters are the log-probabilities P(word | class), and they can be inspected. For each class we look for the words that distinguish it most from the others, that is with the maximum difference between the log-probability in the class and the mean in the others.

pipe = make_pipeline(CountVectorizer(min_df=2), MultinomialNB(alpha=0.1)).fit(train.data, train.target)
vectorizer, nb = pipe[0], pipe[1]
vocabulary = np.array(vectorizer.get_feature_names_out())
log_prob = nb.feature_log_prob_                                       # (classes, words): log P(word | class)

for k, name in enumerate(train.target_names):
    others = np.delete(log_prob, k, axis=0).mean(axis=0)
    top = np.argsort(-(log_prob[k] - others))[:8]
    print("%-20s %s" % (name, ", ".join(vocabulary[top])))

comp.graphics        jpeg, tiff, 3d, animation, visualization, polygon, vga, viewer
rec.autos            tires, brake, toyota, honda, sho, brakes, mustang, car
sci.space            lunar, spacecraft, orbital, orbit, moon, orbiter, titan, launches
talk.politics.guns   firearms, firearm, militia, amendment, handgun, assault, nra, rkba

What the code does. feature_log_prob_ is the matrix (4 classes × words of the vocabulary) of the logarithms of P(word | class), with smoothing. For each class k we subtract the mean of the other three and take the eight words with the highest difference. The list is perfectly readable: for graphics "jpeg", "tiff", "3d", "polygon" appear; for cars "tires", "brake", "toyota", "mustang"; for space "lunar", "spacecraft", "orbit", "moon"; for guns "firearms", "militia", "amendment", "nra". This transparency is a real practical advantage: you see at once whether the model relies on sensible words or on artefacts (proper names, codes, words that depend on how the data were collected). If among the top words you found e-mail addresses or university names, you would know you have a leakage, as in the previous comparison.

The best-known flaw: the probabilities are too sure

We have already seen with the duplicated features that the independence assumption inflates the probabilities. Let us measure it on real text with a reliability table: we group the test predictions by predicted maximum probability and compare, in each group, the mean confidence of the model with the fraction of predictions that are actually correct. A calibrated model has the two columns equal.

proba = pipe.predict_proba(test.data)
confidence = proba.max(axis=1)
correct = proba.argmax(axis=1) == test.target

print("predicted confidence    predictions   mean confidence   fraction correct")
edges = [0, 0.5, 0.7, 0.9, 0.99, 0.999, 1.0001]
for lo, hi in zip(edges[:-1], edges[1:]):
    in_bin = (confidence >= lo) & (confidence < hi)
    print("   [%.3f, %.3f)      %8d        %.3f              %.3f" % (lo, min(hi, 1.0), in_bin.sum(), confidence[in_bin].mean(), correct[in_bin].mean()))

document = vectorizer.transform([test.data[0]])
n_tokens = int(document.sum())
log_likelihood = (document.toarray()[0] * log_prob[test.target[0]]).sum()
print("\nfirst test document: %d words; log-likelihood %.1f; as a probability, in floating point: %s" % (n_tokens, log_likelihood, np.exp(log_likelihood)))

predicted confidence    predictions   mean confidence   fraction correct
   [0.000, 0.500)            77        0.300              0.338
   [0.500, 0.700)            39        0.589              0.436
   [0.700, 0.900)            58        0.816              0.517
   [0.900, 0.990)            80        0.958              0.600
   [0.990, 0.999)            68        0.995              0.676
   [0.999, 1.000)          1221        1.000              0.956

first test document: 165 words; log-likelihood -1132.1; as a probability, in floating point: 0.0

What the code does. The 1,543 predictions are grouped into six intervals of maximum confidence. The result is a textbook example of overconfidence. In the interval 0.99-0.999, where the model says it is 99.5% sure, it is right in only 67.6% of the cases. In the 0.9-0.99 one (mean confidence 0.958) the accuracy is 0.60. Only in the extreme interval, above 0.999, which contains the bulk of the predictions (1,221 out of 1,543), is the fraction correct 0.956: high but still below the 100% the model declares.

The problem is the very independence assumption: the words of a text are strongly correlated (if "orbit" is there it is very likely that "spacecraft" is there too), but Naive Bayes counts them as independent clues and multiplies the likelihoods, accumulating excess evidence.

Practical consequence: do not use the probabilities of Naive Bayes as risk estimates (for example to rank or to apply cost thresholds as in the lesson on metrics). If you need reliable probabilities, you recalibrate (CalibratedClassifierCV) or switch to a discriminative model such as logistic regression.

The last line also shows why we work with logarithms: a message of 165 words has a log-likelihood of −1132, and exp(−1132) in floating point is exactly 0.0, because below about 10⁻³²⁴ the number is not representable. Without logarithms, all the classes would have probability zero and the classifier would not know how to choose.

When Naive Bayes makes sense today

After these experiments the assessment can be precise. Naive Bayes is an excellent baseline for text classification: it trains with a single pass over the data (just count), it holds up well with little material, and it can be read. It is also a good choice when you must update the model incrementally (you add counts, without retraining) or when computing resources are minimal. It is not the right choice, instead, when reliable probabilities are needed (see the reliability table), when the features are highly correlated and their interaction is decisive (short texts where word order matters, or tabular data with complex relationships), or when you are chasing the last point of accuracy. The recommended workflow is simple: start from Naive Bayes as a reference, and adopt a more sophisticated model only if validation, with several seeds and an estimate of the uncertainty, shows an advantage that justifies the extra complexity.

Three exercises

Exercise 1. Try TfidfVectorizer with MultinomialNB in place of CountVectorizer: TF-IDF reduces the weight of very common words; does the accuracy change? Exercise 2. Add stop_words="english" to the vectoriser (it removes "the", "and"...): and the list of indicative words? Exercise 3. Apply CalibratedClassifierCV(pipe, method="isotonic", cv=3) and redo the reliability table: do the two columns get closer?

The code exercises of this lesson are in a single file, esercizi.py (about 15 KB). Three exercises on a new domain, the diagnosis of 3D printer faults from workshop logs: Bayes with the base rate, a multinomial Naive Bayes written from scratch with parametric smoothing, and the measurement of how much a duplicated column inflates the model's confidence. It requires Python 3.12+ and NumPy 2.3. Run it with python3 esercizi.py: it also runs with the parts to complete left empty and tells you which ones are missing.

Other lessons in the series that link to this one: Probability: Reasoning About Uncertainty, Random Forest: why a hundred trees err less than one, Clustering: K-Means and DBSCAN without labels.

In the next lesson we change philosophy: instead of estimating probabilities, we will look for the boundary that separates the classes with the widest possible margin.