Naive Bayes

What this chapter covers

Last chapter we drew one line through one feature and ignored the other seven. This chapter we stop ignoring them. Naive Bayes looks at every measurement a patient has, asks each one "how surprised are you to see this value in a diabetic versus a healthy person," and adds up the surprise. The class with the least surprise wins. That's the whole model, and it's built on a probability rule that's older than every other algorithm in this course by two centuries.

We build it by hand in NumPy — priors, per-class Gaussians, a log-posterior, an argmax — then let scikit-learn's GaussianNB build the same thing and watch the two land on the identical prediction for all 154 held-out patients. In between you get to see the model made literal: the fitted bell curves it carries around, and an animation where the probability visibly tips from "healthy" to "diabetic" as we walk a query point across the data.

The data is the same Pima Indians Diabetes set from week 1: 768 patients, eight measurements each, a label for diabetes. This time all eight columns are in play.

A bit of history

The rule at the center of this is Bayes' theorem, and it comes from a Presbyterian minister. Thomas Bayes worked out how to update a belief in light of evidence sometime in the 1740s, never published it, and died in 1761. His friend Richard Price found the essay in his papers and read it to the Royal Society in 1763. Pierre-Simon Laplace, working independently, put the theorem on solid footing a few years later and used it to do real science with it. So the math we're about to lean on predates computers, statistics as a discipline, and Darwin.

The "naive" part is much younger and much more pragmatic. Somewhere in the 1960s people building automatic document classifiers noticed that if you pretend every word in a document is independent of the others — obviously false, "machine" and "learning" travel together — the Bayes computation collapses into something you can actually run: multiply a pile of per-word probabilities and pick the biggest. It was wrong and it worked anyway. By the 1990s naive Bayes was the workhorse behind spam filters, and Paul Graham's 2002 essay "A Plan for Spam" put a naive-Bayes filter in front of a generation of programmers. It is still the first thing you reach for when you need a text classifier and a result by lunch.

The intuition

Forget the theorem for a second. You already do this. A patient walks in with a very high glucose reading and a high BMI. Your gut says diabetes, and the reason your gut says it is that you've seen the distributions: high glucose is common among diabetics and rare among healthy people, so a high reading is evidence, and two pieces of evidence pointing the same way are stronger than one.

Naive Bayes writes that gut feeling down as arithmetic. For each class it keeps a little statistical portrait of what "normal" looks like — the average glucose of a diabetic, how much it varies, the same for BMI, for age, for all eight features. A new patient gets scored against both portraits: how well do they fit the diabetic profile, how well the healthy one. The catch, the naive part, is that it scores each feature on its own and then just adds the scores, as if glucose and BMI carried completely separate information. They don't — they're correlated — but pretending they're separate is what makes the whole thing a few lines of NumPy instead of a covariance matrix you have to invert.

Here's the two features we'll watch it work on, glucose against BMI, each patient colored by whether they had diabetes:

The two clouds overlap heavily — this is not clean, separable data — but the purple diabetic cloud sits up and to the right. That lean is the whole signal, and Naive Bayes turns it into a pair of probabilities for every point.

The math

Start with Bayes' theorem. We want the probability of a class $y$ given a patient's measurements $x = (x_1, \dots, x_D)$:

Read it right to left. $P(y)$ is the prior — the class frequency before we look at the patient. $P(x \mid y)$ is the likelihood — how probable this exact set of measurements is if the patient belongs to class $y$. $P(x)$ is the same for every class, so it can't change which class wins; we drop it and keep a proportionality:

The hard term is $P(x \mid y)$ — the joint probability of all eight features at once, which in general you cannot estimate without a mountain of data. This is where "naive" earns its name. Assume the features are conditionally independent given the class, and the joint factors into a product of one-feature terms:

That is the entire trick. Each feature contributes its own likelihood, and we multiply them. For continuous features like glucose we need a shape for $P(x_j \mid y)$, and the default choice — the one that makes this Gaussian Naive Bayes — is a normal distribution fitted per feature per class:

where $\mu_{cj}$ and $\sigma_{cj}^{2}$ are just the mean and variance of feature $j$ among the training patients of class $c$. Multiplying eight small densities together underflows toward zero fast, so we take logs and turn the product into a sum:

The prediction is the class that maximizes this:

Every symbol here becomes one line of code, and the sum is the part to keep your eye on — it's where the eight pieces of evidence get combined.

What it's good at, what it isn't

The good side is speed and thrift. Training is a single pass to compute a mean and a variance per feature per class — no iteration, no gradient, nothing to converge. It needs very little data to estimate those handful of numbers, so it holds up when you have twenty examples per class instead of twenty thousand. On text, where every word is a feature and there are tens of thousands of them, it's still one of the fastest useful classifiers you can train, and it remains the baseline a spam filter or a topic tagger has to beat. Reach for it when you want a real answer immediately and a floor under whatever you build next.

The bad side is the naive assumption itself. Features are almost never independent — glucose and BMI move together, and Naive Bayes double-counts correlated evidence as if it were two separate witnesses when it's really one witness talking twice. That makes its probability estimates poorly calibrated: it's prone to shoving the posterior toward a confident 0.99 when the honest answer is 0.7. The ranking usually survives this (the right class still tends to come out on top, which is why accuracy holds up), but do not read a Naive Bayes probability as a real probability. And the Gaussian flavor specifically assumes each feature is bell-shaped within a class; on a hard-skewed feature that assumption is a lie and the fit suffers.

The data

The Pima set has eight features: pregnancies, glucose, blood pressure, skin thickness, insulin, BMI, a diabetes-pedigree score, and age. All eight go into the headline model. For the two-dimensional pictures below I use glucose and BMI, dropping the handful of rows where either is recorded as zero (the dataset's stand-in for "missing"), because a zero would drag a Gaussian's mean and variance somewhere untrue. The scatter above is that clean two-feature slice; the accuracy face-off at the end uses all eight features, raw.

Build it, one function at a time

Seven short functions, no library. This is the order I'd actually type them: the priors first, then the per-class statistics, then the Gaussian, then the posterior that ties them together.

The prior is the easiest thing in the whole course — what fraction of the training patients belong to each class:

def class_priors(y, classes):
    """P(y = c) for each class — just its share of the training rows.

    The prior is what you'd guess with no features at all: if 35% of
    patients are diabetic, the prior leans 65/35 before we look at a single
    measurement.
    """
    return np.array([(y == c).mean() for c in classes])

Then the model itself. For each class, slice out its rows and take the mean and variance of every column. That pair of arrays is the entire fitted model — there is nothing else stored:

def class_stats(X, y, classes):
    """Per-class, per-feature Gaussian mean and variance — the whole model.

    For each class we slice out its rows and take the column means and
    variances. means[c, j] and vars[c, j] are the center and spread of
    feature j *within* class c. That pair of numbers per feature per class is
    everything Naive Bayes ever stores.
    """
    means = np.array([X[y == c].mean(axis=0) for c in classes])   # (C, D)
    variances = np.array([X[y == c].var(axis=0) for c in classes])  # (C, D)
    return means, variances

Now the Gaussian, in log form because that's the only form we'll use it in. Given a value, a mean, and a variance, it returns the log-density, and it broadcasts so one call scores a whole matrix against a whole class at once:

def gaussian_logpdf(x, mean, var):
    """Log of the Gaussian density N(x; mean, var), elementwise.

    We never want the raw density — a product of densities underflows fast.
    We want its log, because in log-space the product of per-feature
    likelihoods becomes a sum we can add up safely.
    """
    return -0.5 * np.log(2.0 * np.pi * var) - 0.5 * (x - mean) ** 2 / var

Here's the heart of it. For each class, add the log prior to the summed per-feature log-likelihoods. That .sum(axis=1) is the naive independence assumption made concrete — eight separate log-likelihoods added up as if they were independent evidence:

def log_posterior(X, priors, means, variances):
    """Unnormalized log posterior log P(y=c) + sum_j log P(x_j | y=c).

    X is (N, D). For each class we add the log prior to the summed
    per-feature log-likelihoods — that sum IS the naive independence
    assumption, treating the features as if they carry separate evidence.
    Returns an (N, C) grid: a score per point per class.
    """
    N, C = X.shape[0], len(priors)
    out = np.zeros((N, C))
    for c in range(C):
        ll = gaussian_logpdf(X, means[c], variances[c]).sum(axis=1)  # (N,)
        out[:, c] = np.log(priors[c]) + ll
    return out

With a score per class in hand, prediction is an argmax across the class axis. We never divide by $P(x)$ because it's identical for every class and can't change the winner:

def predict(X, priors, means, variances, classes):
    """Label each row of X by the class with the highest log posterior.

    argmax over the (N, C) score grid. The normalizing constant P(x) is the
    same for every class, so it can't change which one wins — we skip it.
    """
    scores = log_posterior(X, priors, means, variances)
    return classes[np.argmax(scores, axis=1)]

Accuracy, same one-liner as always:

def accuracy(y_true, y_pred):
    """Fraction of predictions that match the truth."""
    return float((np.asarray(y_true) == np.asarray(y_pred)).mean())

And a thin wrapper that estimates the whole model at once. The only subtlety is a variance floor — a feature with zero variance inside a class would make the Gaussian divide by zero, so we add a tiny epsilon, exactly the way scikit-learn does, which is what lets our numbers match theirs to the decimal:

def fit(X, y, var_smoothing=1e-9):
    """Estimate the whole model: classes, priors, means, variances.

    The only wrinkle is var_smoothing. A feature that never varies inside a
    class has zero variance, and the Gaussian divides by it. So — exactly as
    scikit-learn does — we add a floor: var_smoothing times the largest
    feature variance in the data, added to every variance. It keeps the math
    finite and matches the library's numbers to the decimal.
    """
    classes = np.unique(y)
    priors = class_priors(y, classes)
    means, variances = class_stats(X, y, classes)
    epsilon = var_smoothing * X.var(axis=0).max()
    return classes, priors, means, variances + epsilon

One extra helper, used only for the captions in the animation below: turn the log scores into actual probabilities with a numerically safe softmax. The prediction never needs this — the argmax doesn't care about the normalizing constant — but it's nice to be able to say "P(diabetes) = 0.73" out loud:

def posterior_proba(X, priors, means, variances):
    """Turn the log-posterior scores into real probabilities per class.

    A softmax over the log scores (subtract the row max first — the standard
    log-sum-exp trick — so exp never overflows). Only needed for reporting,
    e.g. "P(diabetes) = 0.82" in a caption; the prediction never needs it.
    """
    scores = log_posterior(X, priors, means, variances)
    scores = scores - scores.max(axis=1, keepdims=True)
    p = np.exp(scores)
    return p / p.sum(axis=1, keepdims=True)

Watch it work

First, the model made visible. These are the fitted Gaussians — the actual bell curves the classifier stores, one per class per feature, for glucose and BMI. This is the model; there is nothing else under the hood:

The two glucose curves are clearly shifted — the diabetic bell sits well to the right, so glucose is a strong witness. The two BMI curves overlap much more, so BMI is a weaker witness that still nudges the vote. Naive Bayes listens to both and weights them automatically by how separated their curves are: a well-separated feature swings the log-likelihood hard, an overlapping one barely moves it.

Now watch the two curves fight over individual patients. Each frame drops one held-out patient onto the glucose-BMI plane. The crosses are the two class centers, the rings are each class Gaussian at one and two standard deviations, and the bar underneath is the posterior — the combined verdict. The point is colored by the predicted class and gets a red ring when the posterior tipped the wrong way. The caption reads off both probabilities.

Play it and watch the bar. Down in the low-glucose corner the posterior is pinned near 0.95 healthy — every feature agrees, no contest. As the query climbs toward high glucose and high BMI the bar slides right, and somewhere in the overlapping middle it crosses 0.5 and the point flips from blue to purple. The red rings cluster exactly there, in the seam where the two Gaussians have nearly equal density, because that's where the evidence is genuinely mixed and a probabilistic model is honest enough to be unsure. This is the thing to take away: the prediction isn't a hard line, it's a tug-of-war between two bells, and the boundary is wherever they pull equally.

The full implementation

The whole file, top to bottom — the seven functions the animation and the face-off both ran on, plus the data loaders:

"""Gaussian Naive Bayes, built from scratch.

Model every feature, for every class, as a 1-D Gaussian. To label a new
point, multiply each class prior by the per-feature likelihoods (the "naive"
assumption that the features are independent given the class), and pick the
class with the biggest product. We do the multiplying in log-space so a
handful of tiny densities never underflows to zero.

Pure NumPy — no ML library anywhere in this file. Every function below
appears in the chapter one step at a time (the `# region:` markers are what
the book's include directives pull in).
"""

import numpy as np
import pandas as pd


# region: class_priors
def class_priors(y, classes):
    """P(y = c) for each class — just its share of the training rows.

    The prior is what you'd guess with no features at all: if 35% of
    patients are diabetic, the prior leans 65/35 before we look at a single
    measurement.
    """
    return np.array([(y == c).mean() for c in classes])
# endregion


# region: class_stats
def class_stats(X, y, classes):
    """Per-class, per-feature Gaussian mean and variance — the whole model.

    For each class we slice out its rows and take the column means and
    variances. means[c, j] and vars[c, j] are the center and spread of
    feature j *within* class c. That pair of numbers per feature per class is
    everything Naive Bayes ever stores.
    """
    means = np.array([X[y == c].mean(axis=0) for c in classes])   # (C, D)
    variances = np.array([X[y == c].var(axis=0) for c in classes])  # (C, D)
    return means, variances
# endregion


# region: gaussian_logpdf
def gaussian_logpdf(x, mean, var):
    """Log of the Gaussian density N(x; mean, var), elementwise.

    We never want the raw density — a product of densities underflows fast.
    We want its log, because in log-space the product of per-feature
    likelihoods becomes a sum we can add up safely.
    """
    return -0.5 * np.log(2.0 * np.pi * var) - 0.5 * (x - mean) ** 2 / var
# endregion


# region: log_posterior
def log_posterior(X, priors, means, variances):
    """Unnormalized log posterior log P(y=c) + sum_j log P(x_j | y=c).

    X is (N, D). For each class we add the log prior to the summed
    per-feature log-likelihoods — that sum IS the naive independence
    assumption, treating the features as if they carry separate evidence.
    Returns an (N, C) grid: a score per point per class.
    """
    N, C = X.shape[0], len(priors)
    out = np.zeros((N, C))
    for c in range(C):
        ll = gaussian_logpdf(X, means[c], variances[c]).sum(axis=1)  # (N,)
        out[:, c] = np.log(priors[c]) + ll
    return out
# endregion


# region: predict
def predict(X, priors, means, variances, classes):
    """Label each row of X by the class with the highest log posterior.

    argmax over the (N, C) score grid. The normalizing constant P(x) is the
    same for every class, so it can't change which one wins — we skip it.
    """
    scores = log_posterior(X, priors, means, variances)
    return classes[np.argmax(scores, axis=1)]
# endregion


# region: accuracy
def accuracy(y_true, y_pred):
    """Fraction of predictions that match the truth."""
    return float((np.asarray(y_true) == np.asarray(y_pred)).mean())
# endregion


# region: fit
def fit(X, y, var_smoothing=1e-9):
    """Estimate the whole model: classes, priors, means, variances.

    The only wrinkle is var_smoothing. A feature that never varies inside a
    class has zero variance, and the Gaussian divides by it. So — exactly as
    scikit-learn does — we add a floor: var_smoothing times the largest
    feature variance in the data, added to every variance. It keeps the math
    finite and matches the library's numbers to the decimal.
    """
    classes = np.unique(y)
    priors = class_priors(y, classes)
    means, variances = class_stats(X, y, classes)
    epsilon = var_smoothing * X.var(axis=0).max()
    return classes, priors, means, variances + epsilon
# endregion


# region: posterior_proba
def posterior_proba(X, priors, means, variances):
    """Turn the log-posterior scores into real probabilities per class.

    A softmax over the log scores (subtract the row max first — the standard
    log-sum-exp trick — so exp never overflows). Only needed for reporting,
    e.g. "P(diabetes) = 0.82" in a caption; the prediction never needs it.
    """
    scores = log_posterior(X, priors, means, variances)
    scores = scores - scores.max(axis=1, keepdims=True)
    p = np.exp(scores)
    return p / p.sum(axis=1, keepdims=True)
# endregion


FEATURES = [
    "Pregnancies", "Glucose", "BloodPressure", "SkinThickness",
    "Insulin", "BMI", "DiabetesPedigreeFunction", "Age",
]


def load_data(path="../data/diabetes.csv"):
    """Pima Indians Diabetes dataset: 768 patients, 8 features, Outcome.

    Returns (X, y): X is (N, 8) float, y is (N,) int (0 = no diabetes,
    1 = diabetes).
    """
    df = pd.read_csv(path)
    X = df[FEATURES].to_numpy(float)
    y = df["Outcome"].to_numpy(int)
    return X, y


def load_two_features(path="../data/diabetes.csv", cols=("Glucose", "BMI")):
    """Two clean features for the 2-D concept views.

    Glucose and BMI both use 0 to mean "not recorded"; those rows would drag
    a Gaussian's mean and variance around, so we drop them here. The headline
    accuracy still uses all eight raw features — this cleaning is only so the
    bell curves and the animation read clearly.
    """
    df = pd.read_csv(path)
    keep = (df[list(cols)] > 0).all(axis=1)
    df = df[keep]
    X = df[list(cols)].to_numpy(float)
    y = df["Outcome"].to_numpy(int)
    return X, y, list(cols)

The library version

Nobody hand-rolls this in production. scikit-learn's GaussianNB is the same model — same priors, same per-class Gaussians, same log-posterior argmax, same variance floor — with the edge cases handled:

def gnb_sklearn(X_train, y_train, X_test):
    """Fit a Gaussian Naive Bayes classifier and predict the test rows.

    The defaults match our scratch version: one Gaussian per feature per
    class, priors read off the training class frequencies, var_smoothing=1e-9.
    """
    clf = GaussianNB()
    clf.fit(X_train, y_train)
    return clf.predict(X_test)

The only reason to look inside is to confirm it really is the same model. Fit it and it exposes exactly the three arrays our fit computes — class_prior_, theta_ (the means), and var_ (the variances, with the same epsilon already folded in):

def gnb_fit(X_train, y_train):
    """Return the fitted model so we can read its estimated parameters.

    clf.class_prior_ is the priors, clf.theta_ the per-class feature means,
    clf.var_ the per-class feature variances (with the smoothing floor
    already folded in). These are the exact numbers our fit() reproduces.
    """
    clf = GaussianNB()
    clf.fit(X_train, y_train)
    return clf

The trace generator asserts our priors, means, and variances match sklearn's to floating-point tolerance, and that the two agree on every single test prediction. When a from-scratch build reproduces the library's fitted parameters exactly, you know you understand the library — there's no hidden step it's doing that you aren't.

Scratch versus library

Eighty/twenty stratified split, seeded, fit on the 614 training patients, scored on the 154 held out. Both models, all eight features:

Both score 0.7727 on the held-out set, and they don't just tie on the average — they make the identical call on all 154 patients. That's the result you want from a from-scratch build: not "close to" the library but bit-for-bit the same model, which is the proof that the four equations up top are the whole algorithm and sklearn isn't hiding anything.

Two honest caveats about that 0.7727. First, the floor is high: guessing "no diabetes" for everyone scores about 65% here, because only about 35% of the patients are diabetic, so the eight-feature model is buying you roughly twelve points over the dumbest possible guess — real, but not miraculous. Second, and this is the Naive Bayes footnote you should never forget, that accuracy comes from the ranking being right, not from the probabilities being trustworthy. The model will hand you a confident-looking 0.9 that you should not bet money on. It sorts patients well; it doesn't price them well.

Takeaways

Naive Bayes is the classifier I reach for when I want an answer before I've finished my coffee. It trains in one pass, it's happy with a few dozen examples per class, and on text — where the features are words and there are tens of thousands of them — it is still the baseline everything else has to beat. If you're standing up a spam filter, a language detector, a first-pass topic tagger, start here and make the fancy model justify replacing it.

The independence assumption is wrong and that's the point. Glucose and BMI are correlated; the model pretends they aren't, double-counts the overlap, and still ranks patients well enough to beat the majority-class baseline by a dozen points. Wrong but useful is a real category in machine learning, and this is its cleanest example. What you don't get for free is calibration — the posterior it prints is a decent sort key and a bad probability, so if a downstream decision actually needs "how sure are we," you calibrate it or you move to a model that estimates probabilities honestly. That's the crack this opens onto the next stretch of the course: logistic regression, coming later, keeps the probabilistic framing but drops the naive assumption and gives you a number you can trust. Naive Bayes is where you learn to think in priors and likelihoods. The models after it are where you learn to stop being naive.