Naive Bayes › Implementation in NumPy and scikit-learn
Machine Learning from ScratchAll chapters Naive Bayes › Implementation in NumPy and scikit-learnSection view · Sections in this chapter7 / 10 ↓
Reading tools Definitions, links & notes
Choose your expertise level to adjust how many terms are explained. Beginners see more tooltips, experts see fewer to maintain reading flow. Hover over underlined terms for instant definitions.
Article links
Make inline references clickable
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.
Gaussian Naive Bayes from scratch
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.
PythonIn [14]Expand
class GaussianNaiveBayes:
"""Gaussian naive Bayes classifier with variance smoothing."""
def __init__(self, var_smoothing=1e-9):
if not np.isfinite(var_smoothing) or var_smoothing < 0:
raise ValueError("var_smoothing must be finite and nonnegative.")
self.var_smoothing = var_smoothing
def fit(self, X, y):
X = np.asarray(X, dtype=float)
y = np.asarray(y)
if (
X.ndim != 2
or X.shape[0] == 0
or X.shape[1] == 0
or y.ndim != 1
or len(y) != len(X)
):
raise ValueError(
"Use a nonempty dense matrix and one label per row."
)
if not np.all(np.isfinite(X)):
raise ValueError(
"Features must be finite; handle missing values explicitly."
)
self.n_features_in_ = X.shape[1]
self.classes_, counts = np.unique(y, return_counts=True)
self.priors_ = counts / counts.sum()
self.means_ = np.zeros((len(self.classes_), X.shape[1]))
self.vars_ = np.zeros_like(self.means_)
for k, label in enumerate(self.classes_):
Xk = X[y == label]
self.means_[k] = Xk.mean(axis=0)
self.vars_[k] = Xk.var(axis=0)
epsilon = self.var_smoothing * X.var(axis=0).max()
self.vars_ += epsilon
if not np.all(np.isfinite(self.vars_)) or np.any(self.vars_ <= 0):
raise ValueError(
"Positive finite class variances are required; remove constant features or choose valid smoothing."
)
self.epsilon_ = epsilon
return self
def _joint_log_likelihood(self, X):
X = np.asarray(X, dtype=float)
if (
X.ndim != 2
or X.shape[1] != self.n_features_in_
or not np.all(np.isfinite(X))
):
raise ValueError(
"Prediction requires a finite matrix with the fitted feature count."
)
jll = np.empty((X.shape[0], len(self.classes_)))
for k in range(len(self.classes_)):
mu = self.means_[k]
var = self.vars_[k]
jll[:, k] = (
np.log(self.priors_[k])
- 0.5 * np.sum(np.log(2.0 * np.pi * var))
- 0.5 * np.sum((X - mu) ** 2 / var, axis=1)
)
if not np.all(np.isfinite(jll)):
raise ValueError(
"Computed scores are nonfinite; inspect feature magnitudes and variances."
)
return jll
def predict_log_proba(self, X):
jll = self._joint_log_likelihood(X)
centered = jll - np.max(jll, axis=1, keepdims=True)
return centered - logsumexp(centered, axis=1, keepdims=True)
def predict_proba(self, X):
return np.exp(self.predict_log_proba(X))
def predict(self, X):
return self.classes_[np.argmax(self._joint_log_likelihood(X), axis=1)]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 , and summing over axis=1 produces one score per observation.
Numerical comparison with scikit-learn
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.
PythonIn [15]Expand
# Generate a three-class problem with two redundant, correlated features.
X_num, y_num = make_classification(
n_samples=600,
n_features=6,
n_informative=4,
n_redundant=2,
n_classes=3,
n_clusters_per_class=1,
class_sep=1.0,
flip_y=0.05,
random_state=123,
)
X_train_num, X_test_num, y_train_num, y_test_num = train_test_split(
X_num, y_num, test_size=0.3, random_state=123, stratify=y_num
)
custom_gnb = GaussianNaiveBayes().fit(X_train_num, y_train_num)
sk_gnb = GaussianNB().fit(X_train_num, y_train_num)
# Predict and collect probabilities on the shared held-out split.
custom_pred = custom_gnb.predict(X_test_num)
sk_pred = sk_gnb.predict(X_test_num)
custom_proba = custom_gnb.predict_proba(X_test_num)
sk_proba = sk_gnb.predict_proba(X_test_num)ConsoleOut [16]Expand
model accuracy log loss from scratch 0.85 0.5086 scikit-learn 0.85 0.5086 Max absolute probability difference: 1.89e-15
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.
Multinomial Naive Bayes from scratch
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.
PythonIn [17]Expand
class MultinomialNaiveBayes:
"""Multinomial naive Bayes classifier with additive smoothing."""
def __init__(self, alpha=1.0, fit_prior=True):
if not np.isfinite(alpha) or alpha <= 0:
raise ValueError("This implementation requires a finite alpha > 0.")
self.alpha = alpha
self.fit_prior = fit_prior
def fit(self, X, y):
X = np.asarray(X, dtype=float)
y = np.asarray(y)
if (
X.ndim != 2
or X.shape[0] == 0
or X.shape[1] == 0
or y.ndim != 1
or len(y) != len(X)
):
raise ValueError(
"Use a nonempty dense matrix and one label per row."
)
if not np.all(np.isfinite(X)):
raise ValueError(
"Features must be finite; handle missing values explicitly."
)
self.n_features_in_ = X.shape[1]
self.classes_, counts = np.unique(y, return_counts=True)
if self.fit_prior:
self.priors_ = counts / counts.sum()
else:
self.priors_ = np.full(len(self.classes_), 1.0 / len(self.classes_))
if np.any(X < 0):
raise ValueError("Multinomial features must be nonnegative.")
self.feature_count_ = np.zeros((len(self.classes_), X.shape[1]))
for k, label in enumerate(self.classes_):
self.feature_count_[k] = X[y == label].sum(axis=0)
self.class_count_ = counts.astype(float)
self.token_count_ = self.feature_count_.sum(axis=1)
smoothed = self.feature_count_ + self.alpha
normalizers = smoothed.sum(axis=1, keepdims=True)
if not np.all(np.isfinite(smoothed)) or not np.all(
np.isfinite(normalizers)
):
raise ValueError(
"Aggregated counts or smoothing overflow; rescale input magnitudes."
)
self.feature_log_prob_ = np.log(smoothed) - np.log(normalizers)
return self
def _joint_log_likelihood(self, X):
X = np.asarray(X, dtype=float)
if (
X.ndim != 2
or X.shape[1] != self.n_features_in_
or not np.all(np.isfinite(X))
or np.any(X < 0)
):
raise ValueError(
"Prediction requires finite nonnegative features with the fitted width."
)
jll = X @ self.feature_log_prob_.T + np.log(self.priors_)
if not np.all(np.isfinite(jll)):
raise ValueError(
"Computed scores are nonfinite; inspect input magnitudes."
)
return jll
def predict_log_proba(self, X):
jll = self._joint_log_likelihood(X)
centered = jll - np.max(jll, axis=1, keepdims=True)
return centered - logsumexp(centered, axis=1, keepdims=True)
def predict_proba(self, X):
return np.exp(self.predict_log_proba(X))
def predict(self, X):
return self.classes_[np.argmax(self._joint_log_likelihood(X), axis=1)]The smoothing shows up in two places. self.feature_count_ + self.alpha is the numerator of the smoothed word probability. The row-wise normalization by normalizers plays the role of in the denominator. With finite positive smoothed counts and normalizers, the probabilities for each class sum to one; the guard rejects aggregate overflow. Storing feature_log_prob_ rather than the probabilities themselves means prediction never multiplies small numbers: X @ feature_log_prob_.T evaluates in a single matrix product.
We now use a tiny document-count example. The vocabulary has four words: "free", "win", "meeting", and "offer". The training set has three spam documents and three ham documents, represented by individual document rows; fitting then aggregates counts within each class. The first test document contains "free" twice and "win" once. The second test document adds one occurrence of "offer". The word "offer" never appears in the ham training counts, so without smoothing the ham likelihood would be zero.
PythonIn [18]Expand
vocab = ["free", "win", "meeting", "offer"]
X_counts = np.array(
[
[4, 5, 1, 3],
[5, 4, 2, 2],
[3, 6, 0, 4],
[1, 1, 8, 0],
[2, 0, 7, 0],
[1, 2, 9, 0],
],
dtype=float,
)
y_counts = np.array(["spam"] * 3 + ["ham"] * 3)
X_counts_test = np.array(
[
[2, 1, 0, 0],
[2, 1, 0, 1],
],
dtype=float,
)
# Fit both the from-scratch and scikit-learn multinomial models with Laplace smoothing.
custom_mnb = MultinomialNaiveBayes(alpha=1.0).fit(X_counts, y_counts)
sk_mnb = MultinomialNB(alpha=1.0).fit(X_counts, y_counts)ConsoleOut [19]Expand
Smoothed word probabilities (alpha = 1):
free win meeting offer
ham 0.1429 0.1143 0.7143 0.0286
spam 0.3023 0.3721 0.0930 0.2326
Predicted spam probabilities:
document P(spam) from scratch P(spam) scikit-learn
free free win 0.9358 0.9358
free free win offer 0.9916 0.9916The three spam rows total , so adding one per word gives over a total of . The ham rows total , giving over a total of after smoothing. The pseudocount of one gives "offer" a nonzero probability in ham even though it never appeared in the ham training counts. Without smoothing, that probability would be zero and the second test document would have zero ham likelihood.
The counts above produce two class-specific word distributions. The figure compares those distributions side by side.
PythonIn [20]Expand
spam_idx_theta = list(custom_mnb.classes_).index("spam")
ham_idx_theta = list(custom_mnb.classes_).index("ham")
theta_spam = np.exp(custom_mnb.feature_log_prob_[spam_idx_theta])
theta_ham = np.exp(custom_mnb.feature_log_prob_[ham_idx_theta])
word_positions = np.arange(len(vocab))VisualizationOut [21]Expand
Look at the two documents side by side. Adding a single "offer" token adds a log word-probability term for each class. Here the smoothed offer probability is larger in spam, so the extra token increases the spam score relative to ham. A zero count contributed no term before; multinomial naive Bayes does not separately penalize word absence.
Effect of smoothing
The smoothing parameter controls the trade-off between the observed counts and a uniform distribution. A very small approaches the unsmoothed maximum likelihood estimate. A large pulls the word probabilities toward and weakens the effect of the data. The table below shows the predicted spam probability for both test documents as varies.
PythonIn [22]Expand
# Sweep the smoothing parameter and record spam probabilities for both documents.
alphas = [0.01, 0.1, 1.0, 10.0]
alpha_rows = []
for alpha in alphas:
model = MultinomialNaiveBayes(alpha=alpha).fit(X_counts, y_counts)
spam_idx = list(model.classes_).index("spam")
probs = model.predict_proba(X_counts_test)[:, spam_idx]
alpha_rows.append(
{
"alpha": alpha,
"P(spam | no offer)": probs[0],
"P(spam | offer)": probs[1],
}
)
alpha_results = pd.DataFrame(alpha_rows)ConsoleOut [23]Expand
alpha P(spam | no offer) P(spam | offer) 0.01 0.9574 0.9999 0.10 0.9555 0.9994 1.00 0.9358 0.9916 10.00 0.7751 0.8548
With , the model is close to the raw counts. The word "offer" still has a tiny ham probability, so the document containing "offer" is strongly but not absolutely classified as spam. As increases, each class token distribution moves toward the uniform distribution; both converge to the same uniform limit. At , the four pseudo-counts total 40 per class, comparable with the observed class token totals of 39 and 31, and the predictions move toward the prior. In practice, we choose by cross-validation rather than by inspecting a single test document.
Over this grid, the document without "offer" moves from about 0.9574 to 0.7751, while the document with "offer" moves from about 0.9999 to 0.8548. The absolute probability change is therefore larger for the first document. The unseen ham token creates a zero-frequency problem without smoothing, but the normalized probability change also depends on the other word ratios and the starting probability.
The four-point table hides the shape of the change. A finer sweep shows how both documents respond as smoothing increases.
PythonIn [24]Expand
alpha_grid = np.logspace(-2, 1.2, 80)
alpha_curve_rows = []
for alpha_value in alpha_grid:
model_curve = MultinomialNaiveBayes(alpha=alpha_value).fit(
X_counts, y_counts
)
spam_idx_curve = list(model_curve.classes_).index("spam")
probs_curve = model_curve.predict_proba(X_counts_test)[:, spam_idx_curve]
alpha_curve_rows.append(
{
"alpha": alpha_value,
"without offer": probs_curve[0],
"with offer": probs_curve[1],
}
)
alpha_curve = pd.DataFrame(alpha_curve_rows)VisualizationOut [25]Expand
Key Parameters
The parameter alpha controls additive smoothing. The scikit-learn default alpha=1.0 is Laplace smoothing. Compare it with smaller values such as or inside validation. Larger values add more pseudo-count mass relative to the observed token counts, pulling the fitted probabilities further toward uniform; whether that improves the intended prediction needs validation. A value of zero reintroduces the zero-frequency problem when a token is absent from a class sample.
The parameter fit_prior determines whether the class priors are learned from the training labels. If it is True, the prior is the observed class frequency. If it is False, the prior is uniform. We can also pass an explicit class_prior array in scikit-learn when we have external information about class frequencies or when we want to correct for class imbalance. For example, in a rare-disease screening problem, we might set the prior to the known population prevalence rather than the prevalence observed in a case-control sample.
For Gaussian naive Bayes, the main parameter is var_smoothing. It adds a small fraction of the largest feature variance to every class-feature variance. The default is . Increasing it stabilizes the model when some features have zero or near-zero variance within a class. Decreasing it makes the model closer to the unsmoothed maximum likelihood estimate. The parameter priors can be used to set class priors explicitly, just as with the multinomial model.
A practical note about scaling: The unsmoothed Gaussian naive Bayes decision rule is invariant to separate invertible affine transformations of each feature. If we replace by , the mean and variance scale accordingly, and the log-density changes by a constant that is the same for every class. The decision is unchanged. However, var_smoothing is computed from the largest feature variance, so features with very different scales can receive uneven effective smoothing. Standardizing features, as discussed in Standardization, makes the smoothing parameter easier to interpret and can improve numerical behavior.
Key Methods
The fit method estimates the model parameters. For Gaussian naive Bayes, it computes class priors, per-class means, and per-class variances. For multinomial naive Bayes, it accumulates word counts and computes smoothed log word probabilities. Both methods return the fitted model, which allows chaining. Both primitive classifiers break an exact score tie by choosing the first class in the sorted np.unique class order.
The predict method returns the most likely class for each observation. It is equivalent to taking the argmax of the joint log-likelihood. The predict_proba method returns the normalized posterior probabilities. The predict_log_proba method returns the logarithms of those probabilities. For extremely small posterior probabilities, predict_log_proba can preserve distinctions across observations that probability outputs lose to underflow.
The scikit-learn implementations also provide partial_fit, which updates the model with a batch of observations without retraining from scratch. This is useful for streaming data and for datasets that do not fit in memory. Our from-scratch classes do not implement partial_fit, but the update equations are straightforward: for Gaussian naive Bayes, we update running means and variances; for multinomial naive Bayes, we add the new counts to the existing count matrix.
The scikit-learn score method returns the mean accuracy on the given data. Accuracy is a useful summary when the classes are balanced, but it can be misleading when they are not. We should pair it with precision, recall, F1, and the confusion matrix, as discussed in Classification Metrics and Decision Thresholds.
Comments
Add to the conversation
A question, a correction, or an example from your own work.
No comments yet. Be the first to share your thoughts!