Naive Bayes › Formula

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.