Part V — Generative Models · Lecture 10

Naive Bayes and the Generative–Discriminative Divide

Instead of modelling p(y|x) directly, model how the data was made and invert with Bayes' theorem. The 'naive' assumption is almost always false, the classifier is often excellent anyway, and understanding why teaches you more about learning than the algorithm itself.

By the end you can

Derive the generative decision rule from Bayes' theorem, and see why the evidence term drops out
Count the parameters the full joint needs, and what conditional independence buys
Derive the MLE for the class prior and the class-conditionals, constraint included
Derive Laplace smoothing as a MAP estimate, not a hack
Prove that Bernoulli naive Bayes is a linear classifier — and explain why it is still not logistic regression
State the bias–variance trade-off between generative and discriminative fitting, and see the crossover

Assumed

Phase A, Lecture 14: Bayes' theorem, MLE, MAP, conjugate priors
Phase B, Lecture 4: where loss functions come from (likelihood → loss)
Phase B, Lecture 5: logistic regression and the log-odds

1Two ways to build a classifier

Every classifier we have built so far attacked \(p(y\mid\vx)\) head-on. Logistic regression posits \(p(y=1\mid\vx)=\sigma(\vw^\top\vx+b)\) and fits \(\vw\) to make the observed labels likely given the observed inputs. It never says anything about where \(\vx\) came from. That is the discriminative strategy.

The generative strategy is the other one. Model the whole data-generating process — first how a class is chosen, then how a point is drawn given that class — and then invert it with Bayes' theorem to get the classification you actually wanted.

Intuition

A discriminative model learns where the border is. A generative model learns what each country looks like, and deduces the border.

Learning the border is a smaller job, so with enough data you do it better. Learning what each country looks like is a bigger job, so you need more data — but you have committed to a structure, and structure is worth a great deal when data is short. That single sentence is the whole trade-off, and the rest of this lecture makes it precise.

Generative classifier

Model the joint distribution by factorising it as

\[ p(\vx, y) \;=\; \underbrace{p(y)}_{\text{class prior}}\;\underbrace{p(\vx\mid y)}_{\text{class-conditional}}, \]

estimate both factors from data, and classify by the posterior \(p(y\mid\vx)\) obtained from Bayes' theorem. A discriminative classifier models \(p(y\mid\vx)\) directly and never estimates \(p(\vx)\).

From the joint model to the decision rule0 / 5

\[\hat y(\vx) \;=\; \argmax_{y\in\{1,\dots,K\}} \; p(y\mid\vx)\]

Start from what we actually want. Phase B Lecture 1 established that under 0–1 loss the risk-minimising prediction is the most probable class. So the target is the posterior, whatever route we take to it.

\[p(y\mid\vx) \;=\; \frac{p(\vx\mid y)\,p(y)}{p(\vx)}\]

Bayes' theorem, exactly as in Phase A Lecture 14. The numerator is the joint \(p(\vx,y)\); the denominator is the evidence \(p(\vx)=\sum_{y'}p(\vx\mid y')p(y')\), which normalises the posterior so it sums to one over classes.

\[\hat y(\vx) \;=\; \argmax_{y}\; \frac{p(\vx\mid y)\,p(y)}{p(\vx)} \;=\; \argmax_{y}\; p(\vx\mid y)\,p(y)\]

The evidence drops out. \(p(\vx)\) does not depend on \(y\), so it is a positive constant across the classes we are comparing. Dividing every candidate by the same positive number cannot change which one is largest. This is why generative classifiers never need to compute the hardest term — but note carefully that it is needed if you want calibrated probabilities rather than just a label.

\[\hat y(\vx) \;=\; \argmax_{y}\; \big[\,\log p(y) + \log p(\vx\mid y)\,\big]\]

Take logs. By Phase A Lecture 2's argmax theorem, \(\log\) is strictly increasing so it preserves the argmax. Products of many small probabilities underflow to zero in floating point; sums of logs do not. Every real implementation works in log space, and this is not optional — with 60 features you underflow immediately.

\[\text{needed: } \; \hat p(y) \;\text{ for } K-1 \text{ free numbers}, \quad \hat p(\vx\mid y) \;\text{ for each } y\]

What remains. The prior is easy: \(K-1\) free parameters. The class-conditional \(p(\vx\mid y)\) is a full distribution over the input space, one per class. That is where all the difficulty of the generative approach lives, and where the naive assumption comes in.

2Why the class-conditional is hard, and what 'naive' buys

Suppose each of the \(d\) features is binary. A general distribution over \(\{0,1\}^d\) assigns a probability to each of the \(2^d\) possible patterns, so it has \(2^d-1\) free parameters — per class. At \(d=30\) that is over a billion numbers per class, from a dataset that might have a thousand rows. There is no version of this that works.

The naive Bayes assumption

Conditional on the class, the features are independent:

\[ p(\vx\mid y) \;=\; \prod_{j=1}^{d} p(x_j \mid y). \]

Read the conditioning bar carefully. This does not say the features are independent. It says that once you know the class, knowing one feature tells you nothing more about another.

Common trap

"Height" and "vocabulary size" are strongly dependent in the general population — both increase with age. Conditional on age = 7 they are nearly independent. That is the shape of the assumption: the class is supposed to explain away all of the dependence.

In real problems it never does exactly. The words New and York co-occur far more than naive Bayes believes, in spam and in non-spam alike. The assumption is false. We use it anyway, and the interesting question — answered at the end of this lecture — is why that so often works.

The parameter count collapses0 / 4

\[\text{full joint:}\quad p(\vx\mid y)\ \text{ needs } \ (m^{d}-1)\ \text{ parameters per class}\]

Count the cells. With \(m\) possible values per feature there are \(m^d\) distinct patterns \(\vx\). A general distribution gives each a probability; they must sum to one, which removes one degree of freedom. Total across \(K\) classes: \(K(m^d-1)\), plus \(K-1\) for the prior.

\[\text{naive:}\quad \prod_{j=1}^{d} p(x_j\mid y)\ \text{ needs } \ d(m-1)\ \text{ parameters per class}\]

Count the cells again, after factorising. Each factor \(p(x_j\mid y)\) is a distribution over \(m\) values: \(m-1\) free numbers. There are \(d\) of them and they are stored separately, so the counts add instead of multiplying. Exponential in \(d\) becomes linear in \(d\).

\[K(m^d-1) \;\longrightarrow\; Kd(m-1) \qquad\text{e.g. } d=30,\,m=2,\,K=2:\quad 2.1\times10^{9} \;\longrightarrow\; 60\]

The size of the win. Two billion parameters to sixty. This is not a small efficiency gain; it is the difference between an estimator that cannot be fit at all and one that can be fit from a few hundred examples.

\[\text{cost: } \ \mathcal{H}_{\text{naive}} \subsetneq \mathcal{H}_{\text{full}} \ \Longrightarrow \ \text{irreducible approximation error}\]

And the price. In Lecture 7's language we have made an enormous cut to the variance by imposing an enormous bias. The naive family cannot represent the true joint unless the true joint happens to factorise. Whatever error that costs, no amount of data removes it. Keep this in mind when the learning curves appear later.

Parameters needed: the full joint versus the naive factorisationInteractive

full joint K(mᵈ−1)naive Kd(m−1)n = 1000 training rows

values per feature m2

number of classes K2

read off at d =12

The vertical axis is log₁₀ of the parameter count, so a straight line means exponential growth. The full joint is a straight line; the naive model is the flat one near the bottom. The grey line marks a thousand training rows: above it you have more parameters than data points, which is Lecture 7's overfitting regime in its most extreme form.

3Fitting it: maximum likelihood, done properly

Take binary features, \(x_j\in\{0,1\}\), and write \(\theta_{jk}=p(x_j=1\mid y=k)\) and \(\pi_k=p(y=k)\). The parameters are \(\{\pi_k\}\) and \(\{\theta_{jk}\}\). We fit them by maximum likelihood — the same principle as Lecture 4, now applied to a joint rather than a conditional model.

MLE for the prior and the class-conditionals0 / 5

\[\ell(\pi,\theta) = \sum_{i=1}^{n}\Big[\log \pi_{y_i} + \sum_{j=1}^{d}\big(x_{ij}\log\theta_{jy_i} + (1-x_{ij})\log(1-\theta_{jy_i})\big)\Big]\]

Write the joint log-likelihood. Each row contributes \(\log p(\vx_i,y_i)=\log p(y_i)+\sum_j\log p(x_{ij}\mid y_i)\) — the first equality is the factorisation of the joint, the second is the naive assumption. The Bernoulli pmf \(\theta^{x}(1-\theta)^{1-x}\) becomes that bracketed sum when logged. Note this is the joint likelihood of \((\vx_i,y_i)\); logistic regression maximises the conditional likelihood of \(y_i\) alone. Remember this difference — it is the whole generative/discriminative split in one line.

\[\max_{\pi} \sum_k n_k \log\pi_k \quad\text{s.t.}\quad \sum_k \pi_k = 1, \qquad n_k=\#\{i: y_i=k\}\]

The prior separates. Collect the \(\log\pi_{y_i}\) terms: each class \(k\) appears \(n_k\) times. Nothing in this piece involves \(\theta\), so we can optimise it alone. But we cannot just set the derivative to zero — the \(\pi_k\) must sum to one, so this is a constrained problem and needs a multiplier.

\[\mathcal{L}=\sum_k n_k\log\pi_k - \lambda\Big(\sum_k\pi_k-1\Big),\qquad \frac{\partial\mathcal{L}}{\partial\pi_k}=\frac{n_k}{\pi_k}-\lambda=0 \;\Rightarrow\; \pi_k=\frac{n_k}{\lambda}\]

Lagrange multiplier. Summing the constraint over \(k\): \(1=\sum_k\pi_k=\frac{1}{\lambda}\sum_k n_k=\frac{n}{\lambda}\), so \(\lambda=n\). Hence \(\hat\pi_k=n_k/n\) — the class frequency, which is what you would have guessed, but now it is derived rather than assumed, and the same argument gives the multinomial MLE anywhere it appears.

\[\frac{\partial \ell}{\partial\theta_{jk}} = \sum_{i:y_i=k}\Big(\frac{x_{ij}}{\theta_{jk}}-\frac{1-x_{ij}}{1-\theta_{jk}}\Big) = 0 \;\;\Longrightarrow\;\; \hat\theta_{jk} = \frac{n_{jk}}{n_k},\quad n_{jk}=\#\{i: y_i=k,\,x_{ij}=1\}\]

Each conditional separates too. \(\theta_{jk}\) appears only in rows of class \(k\) and only in feature \(j\), so the \(Kd\) parameters decouple into \(Kd\) independent one-parameter problems. Clearing denominators: \(n_{jk}(1-\theta)=(n_k-n_{jk})\theta \Rightarrow \theta=n_{jk}/n_k\). Fitting naive Bayes is counting. There is no iterative optimisation anywhere, which is why it is orders of magnitude faster to fit than logistic regression.

\[\frac{\partial^2\ell}{\partial\theta_{jk}^2} = -\frac{n_{jk}}{\theta_{jk}^2}-\frac{n_k-n_{jk}}{(1-\theta_{jk})^2} \;<\;0\]

Confirm it is a maximum. Both terms are negative for \(\theta\in(0,1)\), so \(\ell\) is strictly concave in each \(\theta_{jk}\) and the stationary point is the unique maximum. Phase A Lecture 7's second-derivative test, doing real work — a stationary point of a likelihood is not automatically the MLE.

4The zero-count catastrophe, and the right way to fix it

Suppose the word quinoa never appears in any spam message in your training set. Then \(\hat\theta_{\text{quinoa},\text{spam}} = 0/n_{\text{spam}} = 0\), and \(\log 0 = -\infty\). Any message containing quinoa now has log-probability \(-\infty\) under the spam class, so it can never be classified as spam — no matter what the other five hundred words say.

One unseen word vetoes all the other evidence. This is not a rounding problem; it is a structural failure of maximum likelihood on sparse counts, and it happens constantly in text. The standard fix is called Laplace smoothing, and it is usually presented as "add one to every count". It deserves better than that, because it is exactly the MAP estimate from Lecture 6.

Laplace smoothing is a MAP estimate with a Beta prior0 / 5

\[\theta_{jk} \sim \mathrm{Beta}(\alpha,\beta), \qquad p(\theta) \propto \theta^{\alpha-1}(1-\theta)^{\beta-1}\]

Put a prior on the parameter. Beta is the natural choice on \([0,1]\): it is the conjugate prior for the Bernoulli, meaning the posterior is again Beta, which is what makes the algebra close in one line. \(\alpha=\beta=1\) is the uniform prior — every value of \(\theta\) equally plausible before seeing data.

\[p(\theta_{jk}\mid \text{data}) \;\propto\; \underbrace{\theta^{\,n_{jk}}(1-\theta)^{\,n_k-n_{jk}}}_{\text{likelihood}}\cdot\underbrace{\theta^{\alpha-1}(1-\theta)^{\beta-1}}_{\text{prior}} \;=\; \theta^{\,n_{jk}+\alpha-1}(1-\theta)^{\,n_k-n_{jk}+\beta-1}\]

Multiply and collect exponents. This is Bayes' theorem for the parameter, Phase A Lecture 14. Conjugacy is visible right here: the product has the same functional form as the prior, so the posterior is \(\mathrm{Beta}(n_{jk}+\alpha,\ n_k-n_{jk}+\beta)\). No integral was needed.

\[\hat\theta^{\text{MAP}}_{jk} \;=\; \argmax_\theta\ p(\theta\mid\text{data}) \;=\; \frac{n_{jk}+\alpha-1}{n_k+\alpha+\beta-2}\]

The mode of a Beta\((a,b)\) is \((a-1)/(a+b-2)\) for \(a,b>1\). Derive it if you have not: set \(\frac{d}{d\theta}[(a-1)\log\theta+(b-1)\log(1-\theta)]=\frac{a-1}{\theta}-\frac{b-1}{1-\theta}=0\), which gives \(\theta(a+b-2)=a-1\).

\[\alpha=\beta=1+\!\alpha_{\!s} \;\;\Longrightarrow\;\; \hat\theta_{jk}=\frac{n_{jk}+\alpha_{\!s}}{n_k+2\alpha_{\!s}}\]

And there is 'add-\(\alpha_s\)' smoothing. With \(\alpha_s=1\) it is literally add-one. So the hack is a MAP estimate with a symmetric Beta prior — exactly the same relationship ridge regression has to a Gaussian prior in Lecture 6. Two different smoothers, one idea: the prior contributes pseudo-counts, and pseudo-counts keep an estimate off the boundary.

\[n_k \to \infty \;\Longrightarrow\; \frac{n_{jk}+\alpha_{\!s}}{n_k+2\alpha_{\!s}} \to \frac{n_{jk}}{n_k}, \qquad n_k = 0 \;\Longrightarrow\; \hat\theta_{jk} = \tfrac12\]

Both limits are what you want. With plenty of data the prior washes out and you recover the MLE; with no data you fall back to the prior mean rather than to a degenerate 0 or 1. Smoothing is a shrinkage estimator, and \(\alpha_s\) is its regularisation strength — tune it by cross-validation like any other hyperparameter.

Check yourself

You increase the smoothing parameter \(\alpha_s\) from 1 to 100 on a text classifier with a 10 000-word vocabulary and 500 training documents. What happens to the fitted model?

Not quite. It affects every count. \(\frac{n_{jk}+\alpha_s}{n_k+2\alpha_s}\) differs from \(n_{jk}/n_k\) for all \(j,k\); the effect is merely largest when \(n_{jk}\) is small.

Correct. With \(n_k\approx250\) per class, \(\alpha_s=100\) contributes 200 pseudo-counts against 250 real ones — the prior is now comparable in weight to the data. Every \(\hat\theta_{jk}\) collapses toward \(\tfrac12\), so \(w_j=\log\frac{\theta_{j1}}{1-\theta_{j1}}-\log\frac{\theta_{j0}}{1-\theta_{j0}}\to 0\) and only the prior term \(\log(\pi_1/\pi_0)\) survives. This is underfitting by regularisation, exactly as with a huge \(\lambda\) in ridge.

Not quite. Regularisation improves test error only up to a point; past the optimum it degrades both. There is a best \(\alpha_s\) and it is found by cross-validation, not by making \(\alpha_s\) large. And training accuracy certainly changes — the parameters changed.

Not quite. Nothing breaks: \(\frac{n_{jk}+100}{n_k+200}\) is a perfectly valid probability in \((0,1)\) for any counts. The estimate is heavily biased, not invalid.

5Naive Bayes is a linear classifier

This is the result that surprises people, and it is worth deriving carefully because it reframes everything. With binary features, the naive Bayes decision function is exactly a linear function of \(\vx\) — the same hypothesis class as logistic regression, reached by a completely different route.

The Bernoulli naive Bayes log-odds are affine in x0 / 5

\[\log\frac{p(y=1\mid\vx)}{p(y=0\mid\vx)} \;=\; \log\frac{\pi_1}{\pi_0} \;+\; \sum_{j=1}^{d}\log\frac{p(x_j\mid y=1)}{p(x_j\mid y=0)}\]

Take the log-odds. The evidence \(p(\vx)\) cancels in the ratio — the same cancellation as in the decision rule, but now it is exact rather than merely irrelevant, so we keep a genuine probability. The product over \(j\) becomes a sum of logs.

\[\log\frac{p(x_j\mid y=1)}{p(x_j\mid y=0)} \;=\; x_j\log\frac{\theta_{j1}}{\theta_{j0}} \;+\; (1-x_j)\log\frac{1-\theta_{j1}}{1-\theta_{j0}}\]

Expand one term using the Bernoulli pmf. \(p(x_j\mid y=k)=\theta_{jk}^{x_j}(1-\theta_{jk})^{1-x_j}\); the log of the ratio is the difference of the logs. Because \(x_j\in\{0,1\}\), exactly one of the two terms is active for any given observation — but written this way the expression is algebraically affine in \(x_j\), which is the point.

\[= \; x_j\underbrace{\Big[\log\frac{\theta_{j1}}{1-\theta_{j1}} - \log\frac{\theta_{j0}}{1-\theta_{j0}}\Big]}_{\textstyle w_j} \;+\; \underbrace{\log\frac{1-\theta_{j1}}{1-\theta_{j0}}}_{\text{constant in } x_j}\]

Regroup: collect the \(x_j\) terms. Write \(x_j\log\frac{\theta_{j1}}{\theta_{j0}}+(1-x_j)\log\frac{1-\theta_{j1}}{1-\theta_{j0}}\), expand the second product and gather the coefficient of \(x_j\): \(\log\frac{\theta_{j1}}{\theta_{j0}}-\log\frac{1-\theta_{j1}}{1-\theta_{j0}} = \log\frac{\theta_{j1}}{1-\theta_{j1}}-\log\frac{\theta_{j0}}{1-\theta_{j0}}\). Each weight is a difference of log-odds: how much more this feature favours class 1 than class 0.

\[\boxed{\;\log\frac{p(y=1\mid\vx)}{p(y=0\mid\vx)} \;=\; \vw^\top\vx + b, \qquad b=\log\frac{\pi_1}{\pi_0}+\sum_{j=1}^{d}\log\frac{1-\theta_{j1}}{1-\theta_{j0}}\;}\]

Sum over \(j\). Every \(x_j\)-coefficient goes into \(\vw\); every constant goes into \(b\) along with the log prior ratio. So Bernoulli naive Bayes produces \(p(y=1\mid\vx)=\sigma(\vw^\top\vx+b)\) — the identical functional form as logistic regression.

\[\text{same } \mathcal{H}, \quad \text{different } \vw : \qquad \vw^{\text{NB}} \ \text{from counting}, \qquad \vw^{\text{LR}} \ \text{from maximising } \textstyle\sum_i \log p(y_i\mid\vx_i)\]

The punchline. Same hypothesis class, different selection rule within it. Naive Bayes picks \(\vw\) by estimating each coordinate separately from counts; logistic regression picks \(\vw\) by jointly optimising the conditional likelihood, which lets correlated features negotiate their weights against each other. Naive Bayes cannot do that — by construction it never lets features talk. This single difference explains everything in the next section.

Where this shows up in ML

Two practical consequences of linearity that people often miss:

1. Naive Bayes cannot learn XOR — or any problem whose boundary is not a hyperplane. Its famed "surprising accuracy" is entirely within the linear class. Nobody advertises this because the problems it is used on (text, high-dimensional sparse counts) are mostly linearly separable anyway.

2. Its probabilities are badly calibrated even when its labels are right. The argmax survives a lot of abuse; the numeric posterior does not. Correlated features each contribute a full independent log-odds term, so the evidence gets counted several times over and the posterior is driven to 0 or 1. If you need a probability rather than a label, calibrate the output (Platt scaling, isotonic regression) or use a discriminative model.

Why naive Bayes is overconfident: duplicated evidence is counted twiceInteractive

naive Bayes posteriorcorrect posteriordecision threshold ½

p(x=1 | y=1)0.75

p(x=1 | y=0)0.45

number of identical copies of the feature5

One informative binary feature, observed to be 1, copied c times into the design matrix. The copies carry no new information, so the correct posterior is flat. Naive Bayes multiplies the same likelihood ratio c times and races to certainty. Real features are rarely exact copies — but they are routinely correlated, which is the same disease in milder form.

6The real trade-off: bias against variance, once more

We now have two ways to choose \(\vw\) in the same linear family. Which is better? The honest answer is it depends on \(n\), and the dependence is a clean instance of Lecture 7's decomposition.

Generative versus discriminative (Ng & Jordan, 2001)

For a pair of models in the same hypothesis class fitted generatively (naive Bayes) versus discriminatively (logistic regression):

Logistic regression has the lower asymptotic error. As \(n\to\infty\) it converges to the best linear classifier for the true distribution; naive Bayes converges to something else, because it is optimising the wrong objective when the independence assumption is false.
Naive Bayes converges faster. It approaches its (worse) limit at a rate of order \(\log d / n\), while logistic regression needs order \(d/n\) — that is, roughly \(O(\log d)\) versus \(O(d)\) training examples to get close.

So the two learning curves cross. Naive Bayes wins on small \(n\); logistic regression wins once \(n\) is large enough, and how large depends on \(d\) and on how badly the independence assumption fails.

Note

The convergence rates above are the paper's asymptotic result and we are quoting rather than deriving them — the proof is a uniform-convergence argument of the kind Phase C covers. Marked as deferred, not skipped. The crossing itself is not deferred: you will reproduce it numerically in the code block below, and the mechanism (fewer effective parameters, biased limit) follows from what we have already derived.

Intuition

Why does the estimator with the false assumption ever win? Because \(d\) separate one-parameter counting problems is a far easier statistical task than one \(d\)-parameter joint optimisation. Each \(\hat\theta_{jk}\) is estimated from all \(n_k\) rows of its class, so its variance is \(O(1/n_k)\) regardless of \(d\). The logistic weights are coupled through the Hessian and must sort out \(d\) correlated directions from the same \(n\) rows.

Naive Bayes is a high-bias, very-low-variance estimator. When \(n\) is small, variance dominates the error, and the biased estimator wins on total error even though it is aiming at the wrong target. That is Lecture 7's central lesson in its purest form.

7Continuous features: Gaussian naive Bayes, LDA, QDA

For real-valued features the natural class-conditional is a Gaussian. Three models sit on a spectrum defined entirely by what you allow the covariance to be:

Gaussian naive Bayes: \(p(\vx\mid y=k)=\mathcal{N}(\vm_k, \mathrm{diag}(\sigma_{k1}^2,\dots,\sigma_{kd}^2))\). Diagonal covariance is the conditional independence assumption for Gaussians — a Gaussian has independent coordinates precisely when its covariance is diagonal. \(2d\) parameters per class.
LDA (linear discriminant analysis): \(p(\vx\mid y=k)=\mathcal{N}(\vm_k,\mS)\) with a shared full covariance. Correlations are modelled, but assumed the same in every class. \(Kd + d(d{+}1)/2\) parameters.
QDA (quadratic discriminant analysis): \(p(\vx\mid y=k)=\mathcal{N}(\vm_k,\mS_k)\), a separate full covariance per class. \(K\big(d+d(d{+}1)/2\big)\) parameters — the most flexible and the hungriest.

Why shared covariance makes the boundary linear0 / 4

\[\log p(\vx\mid y=k) = -\tfrac12(\vx-\vm_k)^\top\mS_k^{-1}(\vx-\vm_k) - \tfrac12\log\det\mS_k - \tfrac{d}{2}\log 2\pi\]

The multivariate Gaussian log-density from Phase A Lecture 14. The last term is the same for every class and will cancel in any comparison, so drop it.

\[\log\frac{p(y=1\mid\vx)}{p(y=0\mid\vx)} = -\tfrac12(\vx-\vm_1)^\top\mS_1^{-1}(\vx-\vm_1) + \tfrac12(\vx-\vm_0)^\top\mS_0^{-1}(\vx-\vm_0) - \tfrac12\log\frac{\det\mS_1}{\det\mS_0} + \log\frac{\pi_1}{\pi_0}\]

Form the log-odds. In general this contains \(\vx^\top(\mS_0^{-1}-\mS_1^{-1})\vx\), a genuine quadratic form: the boundary is a conic — ellipse, hyperbola, or pair of lines. That is QDA, and where its name comes from.

\[\mS_1=\mS_0=\mS \;\Longrightarrow\; -\tfrac12\vx^\top\mS^{-1}\vx + \tfrac12\vx^\top\mS^{-1}\vx = 0\]

Impose the shared covariance and the quadratic terms cancel exactly. Expand both quadratics: each produces the identical \(-\tfrac12\vx^\top\mS^{-1}\vx\), with opposite signs. The \(\log\det\) term cancels too, since the determinants are now equal.

\[\log\frac{p(y=1\mid\vx)}{p(y=0\mid\vx)} = \underbrace{\big[\mS^{-1}(\vm_1-\vm_0)\big]^\top}_{\vw^\top}\vx \;\underbrace{-\;\tfrac12\big(\vm_1^\top\mS^{-1}\vm_1-\vm_0^\top\mS^{-1}\vm_0\big)+\log\tfrac{\pi_1}{\pi_0}}_{b}\]

What survives is affine. The cross terms \(\vx^\top\mS^{-1}\vm_k\) are linear in \(\vx\); \(\mS^{-1}\) is symmetric so \(\vx^\top\mS^{-1}\vm_k=\vm_k^\top\mS^{-1}\vx\), and they combine into a single \(\vw\). Hence linear discriminant analysis. Note \(\vw=\mS^{-1}(\vm_1-\vm_0)\), not \(\vm_1-\vm_0\): the inverse covariance rotates the mean-difference toward directions where the within-class noise is small. Gaussian naive Bayes, using a diagonal covariance, cannot perform that rotation at all — and the code below shows exactly how much that costs.

Diagonal covariance versus shared full covariance: what the naive assumption throws awayInteractive

class 0class 1Gaussian NB boundaryLDA boundary

within-class correlation ρ0.85

direction of the mean shift (degrees)0

separation ‖m₁ − m₀‖2.0

Both classes share one covariance with correlation ρ. At ρ = 0 the two boundaries coincide — the naive assumption is true and costs nothing. Turn ρ up and they separate, except when the mean shift happens to lie along a principal axis of the covariance (0° or 90° at ρ ≠ 0 puts it along an eigenvector, where diagonal is accidentally right). Maximum damage is at 0° with ρ large: LDA rotates its weight vector toward the low-noise direction and Gaussian NB cannot.

Practical engineering note

When to actually reach for naive Bayes in 2026. It is not a museum piece, but its niche is specific:

Text and other very high-dimensional sparse counts, where \(d\gg n\) and the counting estimator's variance advantage is decisive. Multinomial NB on TF-IDF is still a respectable baseline.
As the baseline you must beat. It fits in one pass over the data with no hyperparameter search beyond \(\alpha_s\). If your gradient-boosted ensemble cannot beat it, something is wrong with your pipeline, and you have learnt that in seconds rather than hours.
Streaming and online settings, since the sufficient statistics are just counts: update them incrementally and the model is always current. No refitting.
Missing features are free. If \(x_j\) is unobserved, drop its factor from the product — the remaining terms are still a valid likelihood. Almost no other model handles missingness this cleanly.

Do not reach for it when you need calibrated probabilities, when the boundary is not linear, or when \(n\) is comfortably larger than \(d\) and a discriminative model can be fitted.

8Implement it and watch the crossover

Naive Bayes from scratch: smoothing, linearity, and the learning-curve crossingPython · NumPy

import numpy as np
np.set_printoptions(precision=3, suppress=True)

# ============ PART 1: Bernoulli naive Bayes, and why smoothing is not cosmetic
K, d = 2, 60
rgp   = np.random.default_rng(0)
pw    = np.stack([rgp.beta(0.4, 3.0, size=d), rgp.beta(0.4, 3.0, size=d)])
prior = np.array([0.6, 0.4])

def sample_indep(n, rg):
    """Features ARE conditionally independent here: the naive model is correct."""
    y = rg.choice(K, size=n, p=prior)
    X = (rg.random((n, d)) < pw[y]).astype(float)
    return X, y

def fit_bnb(X, y, alpha=1.0):
    pi    = np.array([(y == k).mean() for k in range(K)])
    theta = np.stack([(X[y == k].sum(0) + alpha) / ((y == k).sum() + 2*alpha)
                      for k in range(K)])
    return np.log(pi), theta          # fitting is pure counting -- no optimiser

def score_bnb(X, logpi, theta):
    with np.errstate(divide='ignore', invalid='ignore'):
        return X @ np.log(theta).T + (1 - X) @ np.log(1 - theta).T + logpi

rg = np.random.default_rng(1)
Xtr, ytr = sample_indep(200, rg)
Xte, yte = sample_indep(4000, rg)

print("--- SMOOTHING IS NOT COSMETIC ---")
for alpha in (0.0, 1.0):
    lp, th = fit_bnb(Xtr, ytr, alpha)
    S    = score_bnb(Xte, lp, th)
    bad  = int((~np.isfinite(S)).any(1).sum())
    acc  = float((np.nan_to_num(S, neginf=-1e300).argmax(1) == yte).mean())
    degen = int((th == 0).sum() + (th == 1).sum())
    print(f"  alpha={alpha}:  degenerate theta estimates {degen:3d}"
          f"   test points with a -inf score {bad:4d}/{len(yte)}"
          f"   accuracy {acc:.4f}")

# ============ PART 2: the log-odds really are affine in x
lp, th = fit_bnb(Xtr, ytr, 1.0)
w = np.log(th[1]/(1-th[1])) - np.log(th[0]/(1-th[0]))       # the derivation, coded
b = np.log((1-th[1])/(1-th[0])).sum() + lp[1] - lp[0]
gap = score_bnb(Xte, lp, th)[:, 1] - score_bnb(Xte, lp, th)[:, 0]
print()
print("--- NAIVE BAYES IS A LINEAR CLASSIFIER ---")
print(f"  max |log-odds - (Xw + b)| over {len(yte)} points: "
      f"{np.abs(gap - (Xte@w + b)).max():.2e}")
print("  Same hypothesis class as logistic regression. Only the fitting rule differs.")

# ============ PART 3: break the assumption, and race the two fitting rules
d2, B = 40, 4                                  # 10 blocks of 4 correlated features
rgq = np.random.default_rng(0)
pz  = np.stack([rgq.beta(1.5, 1.5, size=d2//B), rgq.beta(1.5, 1.5, size=d2//B)])
pf  = np.stack([rgq.beta(1.0, 1.5, size=d2),    rgq.beta(1.0, 1.5, size=d2)])

def sample_dep(n, rg):
    # Within a block, features are copies of a latent bit 88% of the time.
    # So they are NOT conditionally independent -- naive Bayes is now misspecified.
    y  = rg.choice(K, size=n, p=np.array([0.5, 0.5]))
    z  = (rg.random((n, d2//B)) < pz[y]).astype(float)
    zx = np.repeat(z, B, axis=1)
    ind = (rg.random((n, d2)) < pf[y]).astype(float)
    return np.where(rg.random((n, d2)) < 0.88, zx, ind), y

def fit_logistic(X, y, lam=1e-2, T=40):
    """Newton / IRLS -- Lecture 05's update, with a small ridge for stability."""
    n, dd = X.shape
    A = np.hstack([X, np.ones((n, 1))]); wv = np.zeros(dd+1)
    R = lam*np.eye(dd+1); R[-1, -1] = 0.0        # never penalise the intercept
    for _ in range(T):
        p  = 1/(1 + np.exp(-A@wv))
        s  = np.clip(p*(1-p), 1e-9, None)
        g  = A.T@(p - y) + R@wv
        H  = (A*s[:, None]).T@A + R + 1e-9*np.eye(dd+1)
        st = np.linalg.solve(H, g); wv -= st
        if np.abs(st).max() < 1e-10: break
    return wv

XteD, yteD = sample_dep(20000, np.random.default_rng(99))
print()
print("--- LEARNING CURVES, ASSUMPTION VIOLATED (mean test accuracy, 25 repeats) ---")
print("   n_train   naive Bayes   logistic   winner")
for n in (10, 20, 40, 80, 200, 600, 2000, 8000):
    an, al = [], []
    for t in range(25):
        rgt = np.random.default_rng(1000 + t)
        Xa, ya = sample_dep(n, rgt)
        if len(np.unique(ya)) < 2: continue
        lpn, thn = fit_bnb(Xa, ya, 1.0)
        an.append(float((score_bnb(XteD, lpn, thn).argmax(1) == yteD).mean()))
        wv = fit_logistic(Xa, ya.astype(float))
        al.append(float(((XteD@wv[:-1] + wv[-1] > 0).astype(int) == yteD).mean()))
    a, l = np.mean(an), np.mean(al)
    print(f"   {n:6d}      {a:.4f}      {l:.4f}   "
          f"{'naive Bayes' if a > l else 'logistic'}")
print("  The curves cross. NB plateaus early at a biased limit; logistic keeps improving.")

# ============ PART 4: continuous features -- diagonal vs full covariance
mu = np.array([[-1.0, 0.0], [1.0, 0.0]])       # mean shift along x1
C  = np.array([[3.0, 2.7], [2.7, 3.0]])        # within-class correlation 0.9
Lc = np.linalg.cholesky(C)

def sample_gauss(n, rg):
    y = rg.integers(0, 2, n)
    return mu[y] + rg.normal(size=(n, 2)) @ Lc.T, y

Xg, yg   = sample_gauss(400,   np.random.default_rng(4))
Xgt, ygt = sample_gauss(20000, np.random.default_rng(7))

means = np.stack([Xg[yg == k].mean(0) for k in range(2)])
varsd = np.stack([Xg[yg == k].var(0)  for k in range(2)])          # GNB: diagonal
Sh    = sum(np.cov((Xg[yg == k] - means[k]).T, bias=True)*(yg == k).sum()
            for k in range(2))/len(yg)                              # LDA: shared full
Sk    = [np.cov(Xg[yg == k].T, bias=True) for k in range(2)]        # QDA: per class

def acc_gnb(X, y):
    s = np.stack([-0.5*(((X-means[k])**2)/varsd[k] + np.log(2*np.pi*varsd[k])).sum(1)
                  for k in range(2)], 1)
    return float((s.argmax(1) == y).mean())
def lin_rule(m, S):
    Si = np.linalg.inv(S)
    return Si@(m[1]-m[0]), -0.5*(m[1]@Si@m[1] - m[0]@Si@m[0])       # the derivation
def acc_lda(X, y):
    wv, bb = lin_rule(means, Sh)
    return float(((X@wv + bb > 0).astype(int) == y).mean())
def acc_qda(X, y):
    s = []
    for k in range(2):
        Si = np.linalg.inv(Sk[k]); D = X - means[k]
        s.append(-0.5*np.einsum('ij,jk,ik->i', D, Si, D) - 0.5*np.log(np.linalg.det(Sk[k])))
    return float((np.stack(s, 1).argmax(1) == y).mean())

wB, bB = lin_rule(mu, C)                        # the true Bayes-optimal rule
print()
print("--- CONTINUOUS FEATURES: what the diagonal assumption costs ---")
print(f"  within-class correlation {C[0,1]/C[0,0]:.2f}, mean shift along x1 only")
print(f"  Gaussian NB (diagonal)  {acc_gnb(Xgt, ygt):.4f}")
print(f"  LDA (shared full)       {acc_lda(Xgt, ygt):.4f}")
print(f"  QDA (per-class full)    {acc_qda(Xgt, ygt):.4f}")
print(f"  Bayes optimal           {float(((Xgt@wB + bB > 0).astype(int) == ygt).mean()):.4f}")
wN = (means[1]-means[0])/varsd.mean(0)
wL, _ = lin_rule(means, Sh)
cosang = wN@wL/(np.linalg.norm(wN)*np.linalg.norm(wL))
print(f"  angle between the GNB and LDA weight vectors: "
      f"{np.degrees(np.arccos(np.clip(cosang,-1,1))):.1f} degrees")
print("  Same data, same Gaussian family. Only the covariance structure differs.")

9Exercises

★

A spam filter uses two binary features: contains the word "free" and contains the word "offer". From 1000 training emails (400 spam, 600 ham):

"free" appears in 320 spam, 60 ham
"offer" appears in 240 spam, 30 ham

Fit a Bernoulli naive Bayes model with Laplace smoothing \(\alpha_s=1\).
Classify an email containing both words. Give the posterior probability of spam.
Compute \(\vw\) and \(b\) from the linearity derivation and check the log-odds agree.

1. Parameters. \(\hat\pi_1=0.4,\ \hat\pi_0=0.6\).

\[\theta_{\text{free},1}=\tfrac{320+1}{400+2}=\tfrac{321}{402}=0.7985,\qquad \theta_{\text{free},0}=\tfrac{60+1}{600+2}=\tfrac{61}{602}=0.1013,\] \[\theta_{\text{offer},1}=\tfrac{240+1}{402}=0.5995,\qquad \theta_{\text{offer},0}=\tfrac{30+1}{602}=0.0515.\]

2. Both words present, so \(\vx=(1,1)\):

\[\log\frac{p(1\mid\vx)}{p(0\mid\vx)} = \log\tfrac{0.4}{0.6} + \log\tfrac{0.7985}{0.1013} + \log\tfrac{0.5995}{0.0515}\] \[= -0.4055 + 2.0646 + 2.4545 = 4.1137.\]

So \(p(\text{spam}\mid\vx)=\sigma(4.1137)=0.9839\).

3. \(w_{\text{free}}=\log\frac{0.7985}{0.2015}-\log\frac{0.1013}{0.8987}=1.3769-(-2.1829)=3.5598\); \(w_{\text{offer}}=\log\frac{0.5995}{0.4005}-\log\frac{0.0515}{0.9485}=0.4034-(-2.9133)=3.3167\).

\(b=\log\frac{0.2015}{0.8987}+\log\frac{0.4005}{0.9485}+\log\frac{0.4}{0.6} =-1.4952-0.8622-0.4055=-2.7628\).

Check: \(\vw^\top(1,1)+b = 3.5598+3.3167-2.7628 = 4.1137\). Identical to part 2, as the derivation guarantees. ∎

Notice how large those weights are relative to the prior term: two words have swamped a prior that favoured ham 3:2. With correlated words — and "free" and "offer" certainly are correlated — this is precisely the double-counting the overconfidence figure showed.

★★

Prove that if the naive Bayes assumption is exactly true and \(n\to\infty\), then Bernoulli naive Bayes and logistic regression converge to the same classifier. Then explain, in one or two sentences, why they nonetheless differ at finite \(n\) — and which you would rather have when \(n=30\) and \(d=500\).

Proof sketch. If the assumption holds, the true conditional is

\[p(y=1\mid\vx)=\sigma\Big(\log\tfrac{\pi_1}{\pi_0}+\textstyle\sum_j\log\tfrac{p(x_j\mid1)}{p(x_j\mid0)}\Big) =\sigma(\vw^{*\top}\vx+b^*)\]

by the linearity derivation, with \(\vw^*,b^*\) built from the true \(\theta_{jk},\pi_k\). So the true conditional lies inside the logistic family — the model is correctly specified for both estimators.

By the law of large numbers \(\hat\theta_{jk}\to\theta_{jk}\) and \(\hat\pi_k\to\pi_k\) almost surely, so naive Bayes' \((\hat\vw,\hat b)\to(\vw^*,b^*)\). Logistic regression maximises the conditional log-likelihood, whose population maximiser is the true conditional whenever that conditional is in the family (the expected log-likelihood is maximised at the truth — a KL-divergence argument, Lecture 4); under the usual identifiability and regularity conditions the MLE is consistent, so it too converges to \((\vw^*,b^*)\). Same limit. ∎

At finite \(n\). They target the same point but along different paths: naive Bayes estimates \(2d+1\) parameters from independent one-dimensional counting problems, each using all \(n_k\) rows of its class; logistic regression solves one coupled \(d\)-dimensional optimisation. The counting estimator has far lower variance, so it gets closer sooner.

At \(n=30,\ d=500\): naive Bayes, without hesitation. Unregularised logistic regression is not even identifiable here (with \(d>n\) the data are separable, the likelihood has no finite maximiser, and the weights diverge — Lecture 5's separation pathology). Naive Bayes just counts, and with smoothing the counts are never degenerate.

★★

Take one binary feature \(x\) and duplicate it \(c\) times, so the design is \((x,x,\dots,x)\). Let \(r=\log\frac{p(x=1\mid y=1)}{p(x=1\mid y=0)}\) be the true log-likelihood ratio for a single observation of \(x=1\), and assume equal priors.

What is the correct posterior log-odds after observing \(x=1\)?
What does naive Bayes report?
Fit logistic regression on the same duplicated design with a ridge penalty \(\lambda\|\vw\|^2\). What does it do to the \(c\) weights, and what is the resulting log-odds? Why is this the right answer?

1. The \(c\) columns are the same random variable, so they carry the evidence of one observation. The correct log-odds is \(r\), independent of \(c\).

2. Naive Bayes treats the columns as conditionally independent, so it multiplies the likelihood ratio \(c\) times: reported log-odds \(=cr\). The odds are overstated by a factor \(e^{(c-1)r}\), which for \(r=0.5,\ c=8\) is \(e^{3.5}\approx33\). With \(c\) large the posterior is numerically 0 or 1 whatever the data says.

3. The likelihood depends on \(\vw\) only through \(\sum_j w_j x = x\sum_j w_j\), so it is a function of \(s=\sum_j w_j\) alone. Fix \(s\); the penalty \(\lambda\sum_j w_j^2\) is then minimised subject to \(\sum_j w_j=s\), and by symmetry (or Cauchy–Schwarz, or a one-line Lagrange multiplier) the minimiser is \(w_j=s/c\) for every \(j\), giving penalty \(\lambda s^2/c\).

So ridge splits the weight equally across the duplicates and the log-odds is \(\sum_j w_j x = s\) — a single number fitted to the actual conditional likelihood, not \(c\) copies of it. Discriminative fitting lets the correlated features negotiate; naive Bayes structurally cannot, because it never looks at more than one feature at a time.

Aside worth noticing: the effective penalty on \(s\) is \(\lambda/c\), so duplicating a feature quietly weakens its regularisation. Redundant columns are not free even for a discriminative model — they are just far less damaging.

★★★

Modify the code block above. Multinomial naive Bayes replaces the Bernoulli class-conditional with a multinomial over word counts: \(p(\vx\mid y=k)\propto\prod_j \theta_{jk}^{x_j}\) with \(\sum_j\theta_{jk}=1\).

Derive the MLE for \(\theta_{jk}\) under the simplex constraint (Lagrange multiplier, as in the lecture) and state its smoothed version.
Implement it and compare against Bernoulli NB on count data.
Is multinomial NB also a linear classifier in \(\vx\)? Prove or disprove.

1. The class-\(k\) log-likelihood contribution is \(\sum_{i:y_i=k}\sum_j x_{ij}\log\theta_{jk} = \sum_j N_{jk}\log\theta_{jk}\) where \(N_{jk}=\sum_{i:y_i=k}x_{ij}\) is the total count of word \(j\) in class \(k\). Maximise subject to \(\sum_j\theta_{jk}=1\):

\[\mathcal{L}=\sum_j N_{jk}\log\theta_{jk}-\lambda\Big(\sum_j\theta_{jk}-1\Big),\qquad \frac{N_{jk}}{\theta_{jk}}=\lambda \;\Rightarrow\; \theta_{jk}=\frac{N_{jk}}{\lambda},\]

and summing the constraint gives \(\lambda=\sum_{j'}N_{j'k}=N_{\cdot k}\). Hence \(\hat\theta_{jk}=N_{jk}/N_{\cdot k}\) — exactly the same Lagrange argument as for the class prior, which is the point of the exercise. Smoothed: \(\hat\theta_{jk}=\frac{N_{jk}+\alpha_s}{N_{\cdot k}+\alpha_s d}\). Note the denominator: \(d\) pseudo-counts total, one per vocabulary word, not \(2\alpha_s\) as in the Bernoulli case.

2. Sample counts with rg.poisson or a multinomial per document and compare. Expect multinomial NB to win on count data and Bernoulli NB to win when only presence/absence matters — each is the correct likelihood for its own data type, which is Lecture 4's whole message.

3. Yes, it is linear. \(\log\frac{p(y=1\mid\vx)}{p(y=0\mid\vx)}=\log\frac{\pi_1}{\pi_0}+\sum_j x_j\log\frac{\theta_{j1}}{\theta_{j0}} = \vw^\top\vx+b\) with \(w_j=\log(\theta_{j1}/\theta_{j0})\) and \(b=\log(\pi_1/\pi_0)\). The multinomial coefficient \(\frac{(\sum_j x_j)!}{\prod_j x_j!}\) depends only on \(\vx\), not on \(k\), so it cancels in the ratio — which is why it can be dropped from the model entirely. Note this \(b\) has no \(\sum_j\log\frac{1-\theta_{j1}}{1-\theta_{j0}}\) term: the Bernoulli model pays for the words that are absent, the multinomial model does not model absence at all. That is a real modelling difference, not a notational one.

Where this goes next

You now have both halves of the classical supervised picture in the linear world: discriminative (Lectures 2, 5, 6) and generative (this one). The remaining two lectures leave linearity behind by two very different routes — Lecture 11 partitions the input space with axis-aligned cuts and combines many weak partitioners, and Lecture 12 keeps a linear boundary but changes the space it lives in.

Previous9 · k-means, Gaussian Mixtures, and EM Next11 · Decision Trees, Random Forests, and Gradient Boosting