Naive Bayes

Two measurements can carry nearly the same information. A classifier that multiplies their class-conditional likelihoods can then count that information twice: the predicted class may stay useful while its reported probability becomes too extreme. Naive Bayes makes this tradeoff explicit. It models each class with simple feature distributions and combines them using a conditional-independence assumption. For text counts, the corresponding multinomial model instead assumes independent token draws given the class and document length.

We will fit both Gaussian measurements and multinomial word counts, compare the primitive calculations with scikit-learn, and inspect a seeded example with nearly duplicated measurements. Training uses class counts and feature summaries; prediction adds one contribution per feature or observed token count. These simple operations make it possible to trace a prediction back to its inputs.

The model is generative rather than discriminative. Logistic regression, which we studied in Logistic Regression, models $p(y \mid \mathbf{x})$ directly. Naive Bayes models $p(\mathbf{x} \mid y)$ and $p(y)$, then combines them with Bayes' theorem. This difference matters: a generative model can marginalize unobserved features under an appropriate observation model, but neither class implemented here silently accepts missing entries. Its probability estimates still depend on its distributional assumptions.

This chapter derives the Bayes decision rule, explains the conditional independence assumption, develops the Gaussian and multinomial event models, shows how smoothing and log-sum-exp keep computations stable, and examines what happens when independence fails. We will implement Gaussian and multinomial naive Bayes from NumPy primitives, compare them with scikit-learn, and study calibration and class imbalance.

The first advantage is speed. Training a Gaussian naive Bayes classifier requires only class counts, per-class feature means, and per-class feature variances. That is $O(nd)$ work for $n$ observations and $d$ features, with no iterative optimization. Multinomial naive Bayes is similarly cheap: it accumulates word counts per class and then normalizes them. Prediction is also fast because each class score is a sum of $d$ terms. For high-dimensional text data, where the number of features can be tens or hundreds of thousands, this matters.

The second advantage is behavior with small samples and many features. A full multivariate Gaussian model would need to estimate a $d \times d$ covariance matrix for each class, whose unregularized sample estimate is singular when a class has at most $d$ observations. Regularized full-covariance models remain possible. Naive Bayes replaces that matrix with $d$ separate variances per class. The number of parameters grows linearly in $d$, not quadratically. This is why naive Bayes often works as a strong baseline for document classification, spam filtering, and other sparse count problems.

The third advantage is interpretability. Each class score is a sum of contributions from the prior and from individual features. We can inspect which words or measurements push an observation toward a class. The model also handles multiclass problems without modification: we compute one score per class and normalize.

Finally, naive Bayes is a natural online learner. With count data, we can update class counts and feature counts as new observations arrive. With Gaussian data, we can update means and variances incrementally. This makes the model useful in streaming settings where retraining a more complex classifier would be expensive.

The central weakness is the conditional independence assumption. Naive Bayes treats features as independent once the class is known. When features remain dependent within a class, multiplying their marginal likelihoods can count overlapping evidence multiple times. A spam classifier may treat "free" and "win" as independent pieces of evidence. If those words remain dependent within spam, multiplying their marginal likelihoods discards that dependence. The result is a classifier that can rank observations well but produce probabilities that are too extreme.

The second weakness is the distributional assumption. Gaussian naive Bayes assumes each feature is normally distributed within each class. If a feature is skewed, multimodal, or bounded, a single Gaussian may be a poor description. Multinomial naive Bayes assumes that word counts follow a multinomial distribution with class-specific word probabilities. That is a reasonable approximation for some text data, but it ignores word order, grammar, and longer-range dependencies.

The third weakness is the zero-frequency problem. If a word never appears in the training documents of a class, the unsmoothed multinomial estimate for that word is zero. A test document containing that word then receives zero likelihood for that class, regardless of all other evidence. Smoothing fixes the numerical problem, but it does not make the model correct.

Class imbalance is another practical concern. The prior $p(y=k)$ can dominate the score when one class is much more frequent than another. A model that predicts the majority class most of the time can have high accuracy but poor recall for the minority class. We need to evaluate with class-specific metrics and consider prior adjustment or threshold tuning.

Finally, naive Bayes cannot capture feature interactions. If the class depends on a combination of features rather than on each feature separately, naive Bayes may fail even with unlimited data. This is a bias problem, not a variance problem. More data will not fix a model that cannot represent the true decision boundary.

This section derives the Bayes decision rule and the log score, then develops the Gaussian and multinomial event models that the implementation uses.

Every probabilistic classifier has to answer the same question: given a vector of measurements, which class is most plausible? Naive Bayes answers it by describing what data from each class look like, weighting those descriptions by how common each class is, and letting Bayes' theorem combine the two. The argument below builds that answer in stages: first the decision rule, then the obstacle that makes the rule impractical, then the assumption that removes the obstacle, and finally the numerical form a program evaluates.

We have a feature vector $\mathbf{x} = (x_1, \ldots, x_d)^\top$ and a class label $y \in \{1, \ldots, K\}$. Under equal-cost zero-one loss, the Bayes decision rule returns a class with the largest posterior probability $p(y=k \mid \mathbf{x})$, using a specified tie convention. Unequal error costs can require a different decision. Bayes' theorem expresses this posterior as

$$p(y=k \mid \mathbf{x}) = \frac{p(\mathbf{x} \mid y=k)\,p(y=k)}{p(\mathbf{x})}.$$

Each of the three factors has a distinct job:

• $p(y=k)$ is the prior. It says how common class $k$ is before we look at any features at all.
• $p(\mathbf{x} \mid y=k)$ is the class-conditional likelihood. It says how plausible this particular measurement vector is if the class really is $k$. This is where the model's description of the data lives.
• $p(\mathbf{x})$ is the evidence. It is the marginal mass or density of the feature vector, with the class averaged out. It has no $k$ in it.

Only the first two factors carry information about which class we are in. The evidence can always be recovered afterwards by summing the numerator over classes:

$$p(\mathbf{x}) = \sum_{l=1}^{K} p(\mathbf{x} \mid y=l)\,p(y=l).$$

Because that denominator is one number shared by every class, it rescales all posteriors by the same amount. Rescaling cannot change which class is largest, so the decision rule needs to look only at the numerator:

$$\hat{y} = \arg\max_{k \in \{1,\ldots,K\}} p(y=k)\,p(\mathbf{x} \mid y=k).$$

Here $\arg\max$ means the value of $k$ that makes the expression as large as possible. We are not yet comparing calibrated posterior probabilities; we are only asking which class produces the largest unnormalized score. The normalized probabilities come back later, and the ranking established here survives that normalization unchanged.

At this point the decision rule is settled, and the real difficulty becomes visible. The likelihood $p(\mathbf{x} \mid y=k)$ is a joint density over $d$ dimensions. Describing it without any structure means modelling how every feature varies with every other feature inside each class. For an unrestricted discrete joint table with $m$ possible values per feature, there are $m^d$ cells. Other structured families have different costs: a full Gaussian covariance, for example, has a quadratic rather than exponential parameter count. Naive Bayes reduces the feature model to class-conditional marginals.

The naive Bayes assumption is

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

This says that, conditional on the class, the features are independent. It does not imply that they are independent overall: mixing classes can create marginal association even when features are independent within every class. Within spam, the factorization says that seeing "free" supplies no further information about whether "win" appears once the class is known. Whether that holds depends on the class-specific data-generating process; Bayes' theorem alone does not supply the assumption. But it is far cheaper to estimate: the joint likelihood collapses from one complicated $d$-dimensional object into a product of $d$ one-dimensional densities, and each of those can be estimated from a single column of the data.

Substituting the product into Bayes' theorem gives

$$p(y=k \mid \mathbf{x}) = \frac{\pi_k \prod_{j=1}^{d} p(x_j \mid y=k)}{\sum_{l=1}^{K} \pi_l \prod_{j=1}^{d} p(x_j \mid y=l)},$$

where $\pi_k = p(y=k)$ is the prior. The class score before normalization is

$$s_k(\mathbf{x}) = \log \pi_k + \sum_{j=1}^{d} \log p(x_j \mid y=k).$$

Taking logarithms is not cosmetic. Products of many small probability masses or densities can underflow floating-point arithmetic. A density can also exceed one, so very large density products can overflow. The logarithm is strictly increasing, so it does not change which class has the largest score. It converts the product into a sum, avoiding a direct product of small factors. The score still requires finite inputs and valid distribution parameters. It also exposes the structure of the model: each feature contributes one additive term, so a score can be decomposed into the prior plus a list of per-feature contributions that we can read off and interpret.

Once we have the scores, the normalized posterior is

$$p(y=k \mid \mathbf{x}) = \frac{\exp(s_k(\mathbf{x}))}{\sum_{l=1}^{K} \exp(s_l(\mathbf{x}))}.$$

For a batch of $n$ observations, we can store the scores in an $n \times K$ matrix $S$, where $S_{ik} = s_k(\mathbf{x}_i)$. Each column corresponds to a class, and each row corresponds to an observation. In code, we compute this matrix with broadcasting rather than writing nested loops.

One ingredient remains: the prior itself. It is usually estimated by its maximum likelihood value,

$$\hat{\pi}_k = \frac{n_k}{n},$$

where $n_k$ is the number of training observations in class $k$ and $n$ is the total number of observations. If a class is rare, we may want to smooth the prior or set it explicitly. With Laplace smoothing, the estimate becomes

$$\hat{\pi}_k = \frac{n_k + \alpha}{n + K\alpha}.$$

For $\alpha>0$, the extra pseudo-counts keep every prior strictly positive. This is the same idea we use for word probabilities later.

The score above is a sum of terms of the form $\log p(x_j \mid y=k)$. Nothing in the derivation so far has said what that quantity is. Specifying it is choosing an event model, and the choice is where the modelling happens: it decides what kind of data the classifier can describe well. The two choices developed here are Gaussian for continuous features and multinomial for count features.

For Gaussian naive Bayes, we assume that feature $j$ within class $k$ follows a normal distribution:

$$x_j \mid y=k \sim \mathcal{N}(\mu_{jk}, \sigma_{jk}^2).$$

where:

• $x_j$: the value of feature $j$ for a given observation
• $k$: the class label
• $\mu_{jk}$: the mean of feature $j$ within class $k$
• $\sigma_{jk}^2$: the variance of feature $j$ within class $k$
• $\mathcal{N}(\mu, \sigma^2)$: the normal distribution with mean $\mu$ and variance $\sigma^2$

The density (probability density function) is

$$p(x_j \mid y=k) = \frac{1}{\sqrt{2\pi\sigma_{jk}^{2}}} \exp\left(-\frac{(x_j-\mu_{jk})^2}{2\sigma_{jk}^2}\right).$$

where:

• $p(x_j \mid y=k)$: the conditional density of feature $j$ given class $k$
• $\mu_{jk}$: the mean of feature $j$ within class $k$
• $\sigma_{jk}^2$: the variance of feature $j$ within class $k$
• $\pi$: the constant $\pi \approx 3.14159$

The factor $1/\sqrt{2\pi\sigma_{jk}^2}$ normalizes the density so that it integrates to one. The exponent measures the squared distance between $x_j$ and the class mean, measured in units of the class standard deviation. A point far from the mean gets a small density. The variance $\sigma_{jk}^2$ controls how wide the distribution is: a large variance makes distant points less penalized.

When a class has observations with positive within-class variation for a feature, its Gaussian maximum likelihood estimates are the sample mean and the variance with divisor $n_k$:

$$\begin{aligned} \hat{\mu}_{jk} &= \frac{1}{n_k}\sum_{i:y_i=k} x_{ij}, \\ \hat{\sigma}_{jk}^2 &= \frac{1}{n_k}\sum_{i:y_i=k} (x_{ij} - \hat{\mu}_{jk})^2. \end{aligned}$$

where:

• $\hat{\mu}_{jk}$: the estimated mean of feature $j$ within class $k$
• $\hat{\sigma}_{jk}^2$: the estimated variance of feature $j$ within class $k$
• $n_k$: the number of training observations in class $k$
• $x_{ij}$: the value of feature $j$ for observation $i$
• $y_i$: the class label of observation $i$, with the sum restricted to observations satisfying $y_i=k$

With positive within-class variation, dividing by $n_k$ gives the variance MLE; for $n_k>1$, dividing by $n_k-1$ gives the unbiased sample variance. If every value of the feature in a class is identical, the formula gives zero, but no maximum is attained over strictly positive Gaussian variances: the likelihood grows without bound as the variance approaches zero. The positive floor used below keeps the fitted density away from that zero-variance boundary. The relative difference between the two variance estimates decreases as $n_k$ grows. This is the same normal maximum likelihood calculation we saw in Maximum Likelihood Estimation.

The implementation uses the log density:

$$\log p(x_j \mid y=k) = -\frac{1}{2}\log(2\pi\sigma_{jk}^2) - \frac{(x_j-\mu_{jk})^2}{2\sigma_{jk}^2}.$$

where:

• $\log p(x_j \mid y=k)$: the log density of feature $j$ given class $k$
• $\mu_{jk}$: the mean of feature $j$ within class $k$
• $\sigma_{jk}^2$: the variance of feature $j$ within class $k$
• $x_j$: the observed value of feature $j$

Let us derive this slowly. Start with the density:

$$p(x) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right).$$

Take the logarithm of the whole expression. The logarithm of a product is the sum of logarithms:

$$\log p(x) = \log\left(\frac{1}{\sqrt{2\pi\sigma^2}}\right) + \log\left(\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\right).$$

The first term is $\log(1) - \frac{1}{2}\log(2\pi\sigma^2) = -\frac{1}{2}\log(2\pi\sigma^2)$. The second term is just the exponent because $\log(\exp(z)) = z$. So we obtain the expression above. The constant term $-\frac{1}{2}\log(2\pi\sigma^2)$ depends on the class through $\sigma_{jk}^2$, so it is not the same for every class and must be included.

In practice, a class may have zero variance for a feature if all its training values are identical. Then the log density contains $\log(0)$ and the squared-distance term divides by zero. We fix this by adding a small positive constant to every variance:

$$\tilde{\sigma}_{jk}^2 = \hat{\sigma}_{jk}^2 + \epsilon, \qquad \epsilon = \gamma \max_{j} \operatorname{Var}(x_{1j}, \ldots, x_{nj}).$$

Here $\gamma$ is a small number such as $10^{-9}$, and $\operatorname{Var}(x_{1j}, \ldots, x_{nj})$ is the variance of feature $j$ over the entire training set. The maximum is taken across features. This is the convention used by scikit-learn's GaussianNB through its var_smoothing parameter. With $\gamma>0$ and at least one feature having positive global training variance, this gives a positive smoothing constant. If every global variance is zero, the formula also gives zero; our implementation rejects any remaining nonpositive class variance instead of claiming smoothing repaired it.

For multinomial naive Bayes, the data are nonnegative integer counts rather than continuous measurements. Conditional on the class and document length, tokens are independent categorical draws. The resulting count components are multinomial and are generally dependent when the length is fixed. The classifier below conditions on length, omitting a separate class-dependent length model. Each observation is a count vector $\mathbf{x} = (x_1, \ldots, x_V)$, where $x_j$ is the number of times word $j$ occurs in a document, and $V$ is the vocabulary size. The document length is $N = \sum_{j=1}^{V} x_j$. The multinomial probability of the document given class $k$ is

$$p(\mathbf{x} \mid y=k) = \frac{N!}{x_1! \cdots x_V!} \prod_{j=1}^{V} \theta_{jk}^{x_j}.$$

where:

• $\mathbf{x}$: the document's count vector, with $x_j$ the number of times word $j$ occurs
• $N$: the document length, $N = \sum_{j=1}^{V} x_j$
• $x_j!$: the factorial of the count of word $j$, used in the multinomial coefficient
• $\prod_{j=1}^{V} \theta_{jk}^{x_j}$: the probability of any particular ordering of the word tokens with those counts
• $\theta_{jk}$: the probability that a randomly chosen token from a document of class $k$ is word $j$
• $V$: the vocabulary size

The fraction is the multinomial coefficient. It counts the number of distinct orderings of the $N$ word tokens when there are $x_j$ tokens of each word. The product gives the probability of each token under the class-specific word distribution. The parameter $\theta_{jk}$ is the probability that a randomly chosen token from a document of class $k$ is word $j$, constrained so that the class-specific word probabilities sum to one.

The multinomial coefficient depends only on the document, not on the class. When we compare classes, it appears in every numerator and in the denominator, so it cancels. Therefore, for classification the likelihood is proportional to the product of class-specific word probabilities, with the proportionality factor independent of the class. This gives

$$p(\mathbf{x} \mid y=k) \propto \prod_{j=1}^{V} \theta_{jk}^{x_j}.$$

Taking logs gives

$$\log p(\mathbf{x} \mid y=k) = C(\mathbf{x}) + \sum_{j=1}^{V} x_j \log \theta_{jk},$$

where:

• $C(\mathbf{x})$: the log multinomial coefficient $\log\!\left(\frac{N!}{x_1! \cdots x_V!}\right)$, which depends only on the document and is the same for every class
• $V$: the vocabulary size
• $x_j$: the count of word $j$ in the document
• $\theta_{jk}$: the probability that a randomly chosen token from a document of class $k$ is word $j$

The constant is the same for all classes and does not affect the decision or the normalized posterior. The implementation below therefore stores log word probabilities and omits the class-independent multinomial coefficient from its class scores.

The maximum likelihood estimate of $\theta_{jk}$ is the fraction of word tokens in class $k$ that are word $j$:

$$\hat{\theta}_{jk} = \frac{N_{jk}}{N_k},$$

where:

• $N_{jk} = \sum_{i:y_i=k} x_{ij}$: the total count of word $j$ in class $k$
• $N_k = \sum_{j=1}^{V} N_{jk}$: the total number of word tokens in class $k$

The unsmoothed estimate requires a positive class token total $N_k$. With smoothing parameter $\alpha>0$, even a class with zero observed tokens receives a defined uniform word distribution. We add pseudo-counts:

$$\hat{\theta}_{jk} = \frac{N_{jk} + \alpha}{N_k + \alpha V}.$$

The numerator adds $\alpha$ to every word count. The denominator adds $\alpha V$ so that the probabilities still sum to one. When $\alpha=1$, this is Laplace smoothing. When $0 < \alpha < 1$, it is sometimes called Lidstone smoothing. When $\alpha=0$, we recover the unsmoothed maximum likelihood estimate, which can be zero and cause the zero-frequency problem.

A related model is the Bernoulli event model. It represents a document by binary presence indicators $x_j \in \{0,1\}$ rather than counts. The likelihood, over the fixed vocabulary of size $V$, is

$$p(\mathbf{x} \mid y=k) = \prod_{j=1}^{V} \theta_{jk}^{x_j}(1-\theta_{jk})^{1-x_j}.$$

Here $\theta_{jk}$ is the probability that word $j$ appears at least once in a document of class $k$. The Bernoulli model explicitly penalizes words that are absent, because each absent word contributes the factor $(1-\theta_{jk})$, while the multinomial model only looks at words that are present: a zero count contributes the factor $\theta_{jk}^{0}=1$ and therefore no evidence. The smoothed estimate is

$$\hat{\theta}_{jk} = \frac{N_{jk} + \alpha}{n_k + 2\alpha},$$

where $N_{jk}$ is the number of documents in class $k$ that contain word $j$, and $n_k$ is the number of documents in class $k$. The denominator has $2\alpha$ because each word contributes two possible outcomes: present or absent.

Two of the formulas above can fail on real data, and both failures are worth understanding rather than patching over. Smoothing solves the statistical failure; log-sum-exp solves the numerical one. They are separate problems that happen to appear in the same chapter.

Smoothing addresses a statistical issue: with limited data, a word or feature value may be missing from a class sample by chance. The unsmoothed estimate sets its probability to zero, which is an extreme claim: it asserts that the event is impossible, based on the absence of evidence rather than on evidence of absence. It also creates a numerical problem, because $\log(0)$ is negative infinity and a single such term wipes out every other piece of evidence in the score. Adding pseudo-counts keeps every probability strictly positive and regularizes the estimate toward a uniform distribution.

For multinomial naive Bayes, the effect of $\alpha$ is easy to see. A larger $\alpha$ pulls every $\theta_{jk}$ toward $1/V$. A smaller $\alpha$ lets the observed counts dominate. If we set $\alpha$ too small, rare words still get very small probabilities and the model can be unstable. As $\alpha$ tends to infinity with the observed counts fixed, every class word distribution tends to uniform and predictions tend toward the prior. In practice, $\alpha=1$ is a reasonable default, and we can tune it with cross-validation. The same principle applies to the Gaussian variance smoothing parameter $\gamma$: it should be small enough not to distort real variance differences, but large enough to prevent division by zero.

Even with strictly positive probabilities, multiplying many numbers can underflow. Suppose a document has 500 words and each word has probability $0.01$. The product is $0.01^{500} = 10^{-1000}$, which is far below the smallest representable positive double. The base-10 logarithm of the product is $-1000$, so the natural logarithm is $-1000 \ln 10 \approx -2302.6$, which is a perfectly manageable floating-point number. This is why we compute scores in log space.

After computing the log scores $s_k(\mathbf{x})$, we need to normalize them into probabilities. A direct computation of $\exp(s_k)$ can overflow if some scores are large, or underflow if all scores are very negative. For finite class scores, log-sum-exp avoids exponentiating their common large scale. Let

$$m = \max_{k} s_k(\mathbf{x}).$$

Then

$$\log \sum_{k=1}^{K} \exp(s_k) = m + \log \sum_{k=1}^{K} \exp(s_k - m).$$

The subtraction makes the largest exponent zero, so the largest term in the sum is $\exp(0)=1$. The sum is therefore between $1$ and $K$, which is numerically safe. The normalized posterior is

$$p(y=k \mid \mathbf{x}) = \exp\left(s_k - \left[m + \log \sum_{l=1}^{K} \exp(s_l - m)\right]\right).$$

For finite scores this is mathematically identical to the original formula. Individual extremely small posterior probabilities may still round to zero, while their log probabilities remain available. The reordering protects the normalization without changing the model. In code, scipy.special.logsumexp implements this calculation.

For Gaussian naive Bayes, the same log-sum-exp step is used after we compute the joint log-likelihood for each class. Positive variances make the Gaussian log-density formula well defined. Extreme finite measurements can still overflow squared residuals; these teaching classes reject nonfinite computed scores. Variance smoothing and score normalization address different numerical problems.

With the formulas in place, it is worth asking what kind of model we have built and where its behaviour comes from. Three properties matter for practice: the model is generative, its decision boundary has a specific geometric shape, and its parameter count is linear in the dimension.

Naive Bayes is a generative classifier. It specifies a full probability model for the features and the label, then uses Bayes' theorem to classify. This is different from a discriminative classifier such as logistic regression, which models the conditional probability of the label directly. Generative models can be trained with unlabeled data in some settings, can handle missing features by marginalizing over them, and can generate synthetic observations. Discriminative models learn the label conditional or decision rule directly and can be more accurate when the generative assumptions are wrong; the comparison still depends on the model family and available training data.

For Gaussian naive Bayes, the decision boundary has a specific shape. The difference between the log scores of two classes $k$ and $l$ is

$$s_k(\mathbf{x}) - s_l(\mathbf{x}) = \log \frac{\pi_k}{\pi_l} - \frac{1}{2}\sum_{j=1}^{d} \log \frac{\sigma_{jk}^2}{\sigma_{jl}^2} - \frac{1}{2}\sum_{j=1}^{d} \left[ \frac{(x_j-\mu_{jk})^2}{\sigma_{jk}^2} - \frac{(x_j-\mu_{jl})^2}{\sigma_{jl}^2} \right].$$

The squared terms make this a quadratic function of $\mathbf{x}$ when the class variances differ. The boundary $s_k(\mathbf{x}) = s_l(\mathbf{x})$ is therefore a quadratic surface. If the variances are shared across classes, so that $\sigma_{jk}^2 = \sigma_{jl}^2$ for every feature $j$, the quadratic terms cancel and the boundary becomes linear. This connects naive Bayes to linear and quadratic discriminant analysis. With a shared diagonal covariance matrix, Gaussian naive Bayes is a diagonal LDA. With class-specific diagonal covariance matrices, it is a diagonal QDA. We will study the full-covariance versions in the next chapter.

The parameter count is another important property. Gaussian naive Bayes has $K$ priors, $Kd$ means, and $Kd$ variances, for $O(Kd)$ parameters. Multinomial naive Bayes has $K$ priors and $KV$ word probabilities, for $O(KV)$ parameters. A full multivariate Gaussian with class-specific covariance matrices would have $O(Kd^2)$ parameters. The linear parameter count is what makes naive Bayes usable in high dimensions.

These counts include all $K$ stored prior probabilities, although their sum-to-one constraint leaves only $K-1$ free prior parameters. Multinomial word probabilities likewise have $V-1$ free parameters per class. The full-covariance count below includes only the distinct entries of each symmetric covariance matrix.

The asymptotic notation becomes clearer when plotted. The next figure compares linear and quadratic parameter growth over a range of feature dimensions.

For fixed feature dimension and positive class probabilities, correctly specified Gaussian or multinomial class models admit consistent parameter estimates under the usual sampling and nondegeneracy conditions. Regularization must become negligible: Gaussian variance additions should tend to zero, while multinomial pseudo-counts can remain fixed because their ratio to growing class token totals then vanishes. With these conditions, the induced rule approaches Bayes risk. A misspecified independence model need not recover the true posterior, although its decision rule can still be useful.

Finally, good classification does not imply good probability estimates. Naive Bayes often ranks observations correctly even when its posterior probabilities are overconfident. If we need calibrated probabilities for risk assessment or decision making, we must check calibration separately. We will return to this point in the practical implications.

The first figure shows decision regions of a two-dimensional Gaussian naive Bayes classifier. Each class gets its own mean and variance for each feature. The background color is the predicted probability of class 1. The decision boundary is where the two class scores tie. Because the class variances differ, the boundary can curve. The model is not using the correlation between the two features; it fits two separate one-dimensional Gaussians per class.

The reliability diagram is a diagnostic, not a proof. In this simulation, the two features are nearly collinear, so the independence assumption is badly violated. The class-conditional marginals are not guaranteed to be Gaussian either, so this single curve does not isolate the effect of dependence. The model can still separate the classes, but its probabilities are not trustworthy as estimates of risk. If we needed calibrated probabilities, we would apply a post-hoc calibration method such as Platt scaling or isotonic regression, which we will discuss in Part XII: Calibration and Prediction Uncertainty.

The formulas are compact, but compact formulas hide their arithmetic. Before writing any code, we work one small Gaussian naive Bayes problem by hand, in full, so that every symbol in the score has a concrete number attached to it. For the Gaussian parameter estimates in this worked example, subscripts give the class first and the feature second; the general Gaussian formulas above use feature first and class second. The training data have two features and two classes, with five observations per class:

The first thing to check is whether anything needs smoothing. Each class has five observations, no feature is constant within a class, and no variance is zero, so the unsmoothed maximum likelihood estimates are safe here and we can use them directly.

The priors are equal because each class has five observations: $\pi_0 = \pi_1 = 0.5$. We want to classify the new point $\mathbf{x} = (4.5, 4.0)$, which sits in the gap between the two clusters but is closer to class 1.

We begin with the sufficient statistics. For class 0, the feature means are

$$\hat{\mu}_{0,1} = \frac{1+2+3+1+2}{5} = 1.8, \qquad \hat{\mu}_{0,2} = \frac{2+3+2+1+2}{5} = 2.0.$$

The variances use division by $n_k=5$:

$$\hat{\sigma}_{0,1}^2 = \frac{(1-1.8)^2 + (2-1.8)^2 + (3-1.8)^2 + (1-1.8)^2 + (2-1.8)^2}{5} = \frac{2.8}{5} = 0.56,$$ $$\hat{\sigma}_{0,2}^2 = \frac{(2-2)^2 + (3-2)^2 + (2-2)^2 + (1-2)^2 + (2-2)^2}{5} = \frac{2.0}{5} = 0.40.$$

For class 1, the means are

$$\begin{aligned} \hat{\mu}_{1,1} &= \frac{6+7+8+7+8}{5} = \frac{36}{5} = 7.2, \\ \hat{\mu}_{1,2} &= \frac{5+6+7+6+5}{5} = \frac{29}{5} = 5.8. \end{aligned}$$

The variances use division by $n_k=5$ in the same way:

$$\begin{aligned} \hat{\sigma}_{1,1}^2 &= \frac{(6-7.2)^2 + (7-7.2)^2 + (8-7.2)^2 + (7-7.2)^2 + (8-7.2)^2}{5} \\ &= \frac{(-1.2)^2 + (-0.2)^2 + (0.8)^2 + (-0.2)^2 + (0.8)^2}{5} \\ &= \frac{1.44 + 0.04 + 0.64 + 0.04 + 0.64}{5} \\ &= \frac{2.8}{5} = 0.56, \end{aligned}$$ $$\begin{aligned} \hat{\sigma}_{1,2}^2 &= \frac{(5-5.8)^2 + (6-5.8)^2 + (7-5.8)^2 + (6-5.8)^2 + (5-5.8)^2}{5} \\ &= \frac{(-0.8)^2 + (0.2)^2 + (1.2)^2 + (0.2)^2 + (-0.8)^2}{5} \\ &= \frac{0.64 + 0.04 + 1.44 + 0.04 + 0.64}{5} \\ &= \frac{2.8}{5} = 0.56. \end{aligned}$$

where:

• $6, 7, 8, 7, 8$: the class 1 values of feature 1 used in the mean
• $5, 6, 7, 6, 5$: the class 1 values of feature 2 used in the mean
• $36$: the sum of the class 1 values of feature 1
• $29$: the sum of the class 1 values of feature 2
• $7.2$: the class 1 mean for feature 1, $\hat{\mu}_{1,1}$
• $5.8$: the class 1 mean for feature 2, $\hat{\mu}_{1,2}$
• $2.8$: the sum of squared deviations for feature 1 and, separately, for feature 2
• $5$: the number of class 1 observations, $n_k$

With the statistics in hand, we evaluate one log density at a time. Each density has the same two-part structure: a normalizing constant that depends only on the variance, and a negative penalty that grows with the squared distance from the mean. For class 0 and feature 1, with $x_1=4.5$:

$$\log p(4.5 \mid y=0) = -\frac{1}{2}\log(2\pi \cdot 0.56) - \frac{(4.5-1.8)^2}{2 \cdot 0.56}.$$

where:

• $4.5$: the value of feature 1 for the new point $\mathbf{x}=(4.5, 4.0)$
• $1.8$: the class 0 mean for feature 1, $\hat{\mu}_{0,1}$
• $0.56$: the class 0 variance for feature 1, $\hat{\sigma}_{0,1}^2$
• $-\log(2\pi \cdot 0.56)/2$: the log-density normalizing contribution from the variance
• $(4.5-1.8)^2/(2 \cdot 0.56)$: the squared distance from the mean measured in units of the variance

The first term is $-\frac{1}{2}\log(2\pi \cdot 0.56) \approx -0.6290$. The squared distance is $(2.7)^2 = 7.29$, and the denominator is $1.12$, so the second term is $-7.29/1.12 \approx -6.5089$. The squared-distance term dominates the normalizing constant here. The new point is far from the class 0 mean, and in this calculation the distance penalty ($6.5089$) is much larger than the variance normalizing constant ($0.6290$). Thus

$$\log p(4.5 \mid y=0) \approx -0.6290 - 6.5089 = -7.1379.$$

For class 0 and feature 2, with $x_2=4.0$:

$$\log p(4.0 \mid y=0) = -\frac{1}{2}\log(2\pi \cdot 0.40) - \frac{(4.0-2.0)^2}{2 \cdot 0.40}.$$

where:

• $4.0$: the value of feature 2 for the new point
• $2.0$: the class 0 mean for feature 2, $\hat{\mu}_{0,2}$
• $0.40$: the class 0 variance for feature 2, $\hat{\sigma}_{0,2}^2$

The first term is $-\frac{1}{2}\log(2\pi \cdot 0.40) \approx -0.4608$. The squared distance is $(4.0-2.0)^2 = 4.0$, and the denominator is $2 \cdot 0.40 = 0.8$, so the second term is $-4.0/0.8 = -5.0$. Thus

$$\log p(4.0 \mid y=0) \approx -0.4608 - 5.0 = -5.4608.$$

The total log score for class 0 is the log prior plus the two log densities, evaluated before rounding:

$$s_0(\mathbf{x}) = \log(0.5) + \log p(4.5 \mid y=0) + \log p(4.0 \mid y=0) \approx -13.2919.$$

Now we repeat the calculation for class 1. For feature 1 the arithmetic is symmetric, because the distance from $4.5$ to $7.2$ is also $2.7$ and the variance is the same:

$$\log p(4.5 \mid y=1) \approx -7.1379.$$

For class 1 and feature 2, with $x_2=4.0$ and mean $5.8$:

$$\log p(4.0 \mid y=1) = -\frac{1}{2}\log(2\pi \cdot 0.56) - \frac{(4.0-5.8)^2}{2 \cdot 0.56}.$$

where:

• $5.8$: the class 1 mean for feature 2, $\hat{\mu}_{1,2}$
• $0.56$: the class 1 variance for feature 2, $\hat{\sigma}_{1,2}^2$

The squared distance is $(-1.8)^2 = 3.24$, so the second term is $-3.24/1.12 \approx -2.8929$. The first term is $-\frac{1}{2}\log(2\pi \cdot 0.56) \approx -0.6290$. Thus

$$\log p(4.0 \mid y=1) \approx -0.6290 - 2.8929 = -3.5219.$$

This is the decisive term. Feature 1 contributes identically to both classes, so the entire difference between the two scores comes from feature 2, where the new point is much closer to the class 1 mean. The total log score for class 1, again evaluated before rounding, is

$$s_1(\mathbf{x}) = \log(0.5) + \log p(4.5 \mid y=1) + \log p(4.0 \mid y=1) \approx -11.3530.$$

Class 1 has the larger score by $s_1 - s_0 = -11.3530 - (-13.2919) = 1.9389$. To convert the scores into a posterior probability, we use the logistic form for two classes:

$$p(y=1 \mid \mathbf{x}) = \frac{\exp(s_1)}{\exp(s_0)+\exp(s_1)} = \frac{1}{1+\exp(s_0-s_1)}.$$

where:

• $s_0$: the unnormalized log score for class 0, computed above as $s_0(\mathbf{x}) \approx -13.2919$
• $s_1$: the unnormalized log score for class 1, computed above as $s_1(\mathbf{x}) \approx -11.3530$
• $s_0-s_1$: the score difference, which is the only quantity that matters after normalization

Substituting $s_0-s_1 = -13.2919 - (-11.3530) = -1.9389$ gives

$$p(y=1 \mid \mathbf{x}) = \frac{1}{1+\exp(-1.9389)} \approx 0.874.$$

So we classify the point as class 1 with a predicted probability of about $87.4\%$. Notice how the log-sum-exp normalization works here in a two-class setting: the difference of scores is the only quantity that matters, and the sigmoid turns that difference into a probability without ever computing the raw likelihoods.

Within class 0, the sample correlation between the two features is about $0.42$. Within class 1, it is about $0.43$. These sample correlations are a diagnostic to investigate; five observations per class cannot establish population conditional dependence. The fitted model nevertheless multiplies the two marginal likelihoods. It favors class 1, while the calibration of its probability remains a separate question.

We can verify the arithmetic with a short computation:

We now implement both event models from primitives. The goal is to see exactly where the formulas enter the code, then compare with scikit-learn to confirm that the implementation is correct. We use NumPy for the linear algebra and scikit-learn only as a reference.

The class below stores the class priors, per-class feature means, and per-class feature variances. The fit method computes the maximum likelihood estimates. The _joint_log_likelihood method evaluates the Gaussian log density for every class and every observation. The predict_log_proba method applies the log-sum-exp normalization.

The teaching classes accept finite dense feature matrices with matching labels. Use nonmissing class labels of one consistent string or numeric type; these small classes do not validate missing or mixed-type labels. The Gaussian class requires positive fitted variances; the multinomial class requires nonnegative values and alpha>0. Integer counts match the multinomial sampling interpretation, while nonnegative fractional weights use its score as an approximation.

Each line of _joint_log_likelihood maps directly onto a piece of the log-score formula. np.log(self.priors_[k]) is the log prior. The term -0.5 * np.sum(np.log(2.0 * np.pi * var)) collects the normalizing constants across features. The term -0.5 * np.sum((X - mu) ** 2 / var, axis=1) is the summed squared-distance penalty, broadcast over every observation at once.

The fit method uses np.var, which divides by the number of observations. That matches the maximum likelihood estimate. The smoothing constant is a fraction of the largest feature variance in the full training set, which matches scikit-learn's convention. The _joint_log_likelihood method uses broadcasting: (X - mu) ** 2 / var has shape $n \times d$, and summing over axis=1 produces one score per observation.

We generate a three-class problem with six features, two of which are redundant. Redundant features are linear combinations of informative features, so they introduce correlation. This is a useful stress test for the independence assumption. We fit both our implementation and scikit-learn's GaussianNB on the same training split and compare accuracy, log loss, and probabilities.

The two implementations agree to floating-point precision. The accuracy and log loss are identical up to rounding. This confirms that the from-scratch implementation follows the same Gaussian event model as scikit-learn. The log loss is a proper scoring rule: it rewards confident correct predictions and penalizes confident wrong ones. It is often more informative than accuracy when we care about probability quality.

The maximum absolute probability difference is worth reading literally. It is a numerical check, not a performance measure: a value near machine epsilon says that two independently written programs computed the same numbers. If it were large, the discrepancy would point to a difference in how means, variances, or smoothing constants were computed, and we would go back to the formulas to find which convention differed.

For count data, we accumulate word counts per class, add the smoothing parameter, and normalize. The _joint_log_likelihood method is a matrix multiplication between the document count matrix and the transposed matrix of log word probabilities. This is one of the reasons multinomial naive Bayes is so fast for text.