Probability and Statistics for Machine Learning

Phase 0 · Week 5–6 · ~20–25 hours

Every machine learning model is, at its core, a probabilistic model. Linear regression assumes Gaussian noise. Logistic regression models a Bernoulli distribution over class labels. Neural language models define a probability distribution over token sequences. When you understand this, the model design decisions that seem arbitrary — the choice of loss function, the form of regularization, the meaning of “softmax output” — become derived necessities rather than engineering folklore. This file builds that probabilistic vocabulary from first principles.

---

Why Probability Is the Language of ML (Not Just a Prerequisite)

The connection is exact, not approximate:

• Mean Squared Error is the negative log-likelihood under a Gaussian noise model

• Cross-entropy loss is the negative log-likelihood under a categorical distribution

• L2 regularization is a Gaussian prior on weights (MAP estimation, not “prevents overfitting”)

• Dropout is approximate variational inference

• Batch normalization can be understood as normalizing intermediate distributions

If you treat these as heuristics (“cross-entropy works well for classification”), you’ll use them correctly most of the time but fail to generalize. If you understand the probabilistic derivation, you can construct appropriate loss functions for novel problems and know when standard choices are wrong.

---

Core Concepts

1. Probability Foundations

The two rules everything derives from:

Sum rule:  P(X) = Σᵧ P(X, Y)
Product rule: P(X, Y) = P(X|Y) · P(Y)

From these two rules, you can derive Bayes’ theorem, marginal probabilities, conditional independence, and everything else. Don’t let statistics courses bury you in formulas before you internalize these two.

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

# Simulate the sum and product rules concretely
# Two fair dice
np.random.seed(42)
n_rolls = 100000
die1 = np.random.randint(1, 7, n_rolls)
die2 = np.random.randint(1, 7, n_rolls)

# Joint probability P(die1=3, die2=4)
joint = np.mean((die1 == 3) & (die2 == 4))
print(f"P(d1=3, d2=4) empirical: {joint:.4f}")
print(f"P(d1=3, d2=4) theoretical: {1/36:.4f}")

# Marginal probability P(die1=3) via sum rule
marginal = np.mean(die1 == 3)
print(f"P(d1=3) empirical: {marginal:.4f}")
print(f"P(d1=3) theoretical: {1/6:.4f}")

# Conditional probability P(sum=7 | die1=3)
sum_7_given_d1_3 = np.mean((die1 + die2 == 7)[die1 == 3])
print(f"P(sum=7 | d1=3) empirical: {sum_7_given_d1_3:.4f}")
print(f"P(sum=7 | d1=3) theoretical: {1/6:.4f}")

---

2. Key Probability Distributions

The distributions that appear most frequently in ML:

Gaussian (Normal) Distribution

X ~ N(μ, σ²) — the foundation of regression, noise models, and many priors.

mu, sigma = 0.0, 1.0
x = np.linspace(-4, 4, 1000)
pdf = stats.norm.pdf(x, mu, sigma)

# Why it's ubiquitous: Central Limit Theorem
# Sum of many independent RVs → Gaussian, regardless of their individual distributions
n_samples = 1000
sample_means = []
for _ in range(n_samples):
    # Sum of 30 Bernoulli(0.5) → approximately Gaussian
    samples = np.random.binomial(1, 0.5, 30)
    sample_means.append(np.mean(samples))

plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(x, pdf)
plt.title('N(0, 1) PDF')
plt.xlabel('x'); plt.ylabel('p(x)')

plt.subplot(1, 2, 2)
plt.hist(sample_means, bins=50, density=True)
# Overlay the theoretical Gaussian
mu_clt = 0.5
sigma_clt = np.sqrt(0.5 * 0.5 / 30)
x_clt = np.linspace(0.2, 0.8, 100)
plt.plot(x_clt, stats.norm.pdf(x_clt, mu_clt, sigma_clt), 'r-', lw=2)
plt.title('CLT: Bernoulli means → Gaussian')
plt.tight_layout()
plt.savefig('clt_demo.png', dpi=150)

ML connections for Gaussian:

• MSE loss derivation: if y = wᵀx + ε where ε ~ N(0, σ²), then maximizing likelihood = minimizing MSE

• Feature scaling: assumes roughly Gaussian inputs for gradient descent efficiency

• Bayesian linear regression: Gaussian prior on weights → L2 regularization

Bernoulli and Categorical Distributions

X ~ Bernoulli(p) for binary outcomes; X ~ Categorical(p₁,...,pₖ) for multi-class.

# Bernoulli: coin flips, binary classification labels
p = 0.3
bernoulli_pmf = {0: 1-p, 1: p}

# Categorical: dice rolls, multi-class classification labels
# The softmax output of a classifier IS a categorical distribution
logits = np.array([2.0, 1.0, 0.5, -0.5])
softmax = np.exp(logits) / np.sum(np.exp(logits))
print(f"Softmax (categorical probs): {softmax}")
print(f"Sum to 1: {softmax.sum():.6f}")

# Numerically stable softmax (ALWAYS use this in practice)
def stable_softmax(logits):
    shifted = logits - np.max(logits)  # prevents overflow
    exp_shifted = np.exp(shifted)
    return exp_shifted / np.sum(exp_shifted)

print(f"Stable softmax: {stable_softmax(logits)}")

Dirichlet Distribution

The Dirichlet is a distribution over probability vectors — the conjugate prior for the Categorical.

# Dirichlet(α): generates probability vectors
# Useful in topic modeling (LDA), Bayesian statistics
alpha = np.array([2.0, 3.0, 1.0])  # concentration parameters
samples = np.random.dirichlet(alpha, size=5)
print("Dirichlet samples (each row sums to 1):")
for s in samples:
    print(f"  {s.round(3)}  sum={s.sum():.3f}")
# Higher α[i] → more probability mass typically near category i

---

3. Bayes’ Theorem — The Central Equation

P(hypothesis | data) = P(data | hypothesis) · P(hypothesis) / P(data)
                     = likelihood × prior / evidence

posterior ∝ likelihood × prior

# Medical test example: concrete Bayes computation
# Disease prevalence: 1%
# Test sensitivity (true positive rate): 95%
# Test specificity (true negative rate): 90%

p_disease = 0.01
p_no_disease = 1 - p_disease
p_positive_given_disease = 0.95
p_positive_given_no_disease = 0.10  # false positive rate

# P(positive) = total probability
p_positive = (p_positive_given_disease * p_disease + 
              p_positive_given_no_disease * p_no_disease)

# Bayes: P(disease | positive test)
p_disease_given_positive = (p_positive_given_disease * p_disease) / p_positive
print(f"P(disease | positive test): {p_disease_given_positive:.4f}")
# ~8.7% — counterintuitive! Rare disease + imperfect test → low PPV

ML framing of Bayes:

• P(y | x): the posterior — what we want to learn (conditional probability of label given input)

• P(x | y): the likelihood — how inputs are generated given a label (Naive Bayes uses this)

• P(y): the prior — class frequency in training data

• Logistic regression directly models the posterior. Naive Bayes models the likelihood. This is the fundamental difference.

---

4. Maximum Likelihood Estimation (MLE)

MLE answers: “Given data, what parameters make this data most probable?”

θ_MLE = argmax_θ P(data | θ) = argmax_θ Σᵢ log P(xᵢ | θ)

# MLE for a Gaussian: estimate μ and σ² from data
np.random.seed(42)
true_mu, true_sigma = 5.0, 2.0
data = np.random.normal(true_mu, true_sigma, 1000)

# MLE estimates (derivation: set gradient of log-likelihood to zero)
mu_mle = np.mean(data)           # ∂/∂μ log L = 0  →  μ̂ = mean
sigma2_mle = np.var(data)        # MLE: divide by N (BIASED!)
sigma2_unbiased = np.var(data, ddof=1)  # Unbiased: divide by N-1

print(f"True μ: {true_mu}, MLE μ̂: {mu_mle:.4f}")
print(f"True σ²: {true_sigma**2}, MLE σ̂²: {sigma2_mle:.4f} (biased)")
print(f"Unbiased σ̂²: {sigma2_unbiased:.4f}")

# WHY cross-entropy is the loss for classification (MLE derivation):
# If labels follow Categorical(softmax(Wx + b)), then
# log P(y=c | x) = log softmax_c(Wx + b)
# -Σᵢ log P(yᵢ | xᵢ) = Cross-Entropy Loss
# MLE ≡ Cross-Entropy Minimization — this is not a coincidence, it's a derivation

---

5. Maximum A Posteriori (MAP) Estimation

MAP adds a prior over parameters:

θ_MAP = argmax_θ P(θ | data) = argmax_θ [log P(data | θ) + log P(θ)]

# MAP for linear regression with Gaussian prior on weights:
# P(w) = N(0, 1/λ · I) → log P(w) ∝ -λ||w||²
# MAP = MLE + L2 regularization term

# L2-regularized linear regression IS MAP with Gaussian prior
# L1-regularized linear regression IS MAP with Laplace prior

# This is not a metaphor — it's exact algebraic equivalence

def map_linear_regression(X, y, lambda_reg=0.1):
    """
    MAP estimate = Ridge regression
    Closed form: w_MAP = (XᵀX + λI)⁻¹ Xᵀy
    Compare to MLE: w_MLE = (XᵀX)⁻¹ Xᵀy
    """
    n, d = X.shape
    w_map = np.linalg.solve(X.T @ X + lambda_reg * np.eye(d), X.T @ y)
    return w_map

np.random.seed(42)
X = np.random.randn(50, 3)
w_true = np.array([1.0, -2.0, 0.5])
y = X @ w_true + 0.5 * np.random.randn(50)

w_mle = np.linalg.lstsq(X, y, rcond=None)[0]
w_map = map_linear_regression(X, y, lambda_reg=0.1)

print(f"True weights: {w_true}")
print(f"MLE weights:  {w_mle.round(4)}")
print(f"MAP weights:  {w_map.round(4)}  (shrunk toward 0)")

---

6. KL Divergence — The Distance Between Distributions

KL(P || Q) = Σₓ P(x) log(P(x) / Q(x))    (discrete)
           = ∫ p(x) log(p(x) / q(x)) dx   (continuous)

Properties:

• KL ≥ 0 always (Gibbs inequality)

• KL(P||Q) = 0 iff P = Q

• NOT symmetric: KL(P||Q) ≠ KL(Q||P)

def kl_divergence(p, q, eps=1e-10):
    """KL(P || Q) for discrete distributions."""
    p = np.array(p) + eps  # prevent log(0)
    q = np.array(q) + eps
    p = p / p.sum()  # normalize
    q = q / q.sum()
    return np.sum(p * np.log(p / q))

# Example: two categorical distributions over 4 classes
P = [0.4, 0.3, 0.2, 0.1]  # "true" distribution
Q = [0.25, 0.25, 0.25, 0.25]  # uniform "model" distribution

kl_pq = kl_divergence(P, Q)
kl_qp = kl_divergence(Q, P)

print(f"KL(P||Q): {kl_pq:.4f}")
print(f"KL(Q||P): {kl_qp:.4f}")
print(f"Not symmetric: {not np.isclose(kl_pq, kl_qp)}")

# Cross-entropy relationship:
# H(P, Q) = H(P) + KL(P || Q)
# Cross-entropy = entropy + KL divergence
# Minimizing cross-entropy ≡ minimizing KL divergence from model Q to true P

def cross_entropy(p, q):
    p, q = np.array(p) + 1e-10, np.array(q) + 1e-10
    p, q = p / p.sum(), q / q.sum()
    return -np.sum(p * np.log(q))

def entropy(p):
    p = np.array(p) + 1e-10
    p = p / p.sum()
    return -np.sum(p * np.log(p))

print(f"\nH(P, Q) = {cross_entropy(P, Q):.4f}")
print(f"H(P) + KL(P||Q) = {entropy(P) + kl_pq:.4f}")
print(f"These are equal (as expected): {np.isclose(cross_entropy(P, Q), entropy(P) + kl_pq)}")

---

What Most Learners Get Wrong

They memorize Bayes’ theorem but never use it to derive their loss function.

The single most powerful exercise in this entire phase: take a model you’ve already trained (logistic regression, linear regression), and derive its loss function from the negative log-likelihood of the probabilistic model it implies. Then ask: what probabilistic model does L2 regularization add? What about L1? This exercise converts probability from “a required topic” into “the design language of ML.”

Second gap: confusing P(y=1|x) with P(x|y=1). These are reversed, and the reversal matters enormously. Logistic regression models the former directly (discriminative). Naive Bayes models the latter (generative). The choice between these two approaches has practical consequences for data efficiency, class imbalance handling, and missing data.

Third gap: not distinguishing between biased and unbiased estimators. MLE for variance divides by N and is biased. Sample variance divides by N-1 and is unbiased. This is not a trivia fact — it’s the reason small datasets systematically underestimate uncertainty in Bayesian models.

---

Practice Problems with Acceptance Criteria

Problem 1 — MLE Derivation Derive the MLE estimator for the parameter λ of a Poisson distribution from scratch (pen and paper). Set ∂log L/∂λ = 0 and solve. Verify numerically by generating Poisson(3) samples and checking that your estimator recovers λ≈3.

• Acceptance criteria: Correct closed-form derivation (λ̂ = x̄). Numerical verification with 1000 samples showing error < 0.1.

Problem 2 — Bayes’ Theorem in ML Implement Naive Bayes from scratch for text classification (use 20Newsgroups dataset, binary: ‘rec.sport.hockey’ vs ‘sci.space’). Do not use sklearn’s GaussianNB or MultinomialNB.

• Acceptance criteria: Accuracy > 85% on test set. Explain in comments why you used log-probabilities instead of raw probabilities.

Problem 3 — KL and Cross-Entropy Relationship Generate data from a mixture of two Gaussians. Fit a single Gaussian using MLE. Compute the KL divergence from true distribution to fitted model. Then fit a mixture model and show the KL decreases.

• Acceptance criteria: Numerical KL divergence computed and printed for both fits. Plot shows true vs fitted distributions.

---

Resources

Resource	Focus	Time	Free?
Sheldon Ross “A First Course in Probability”	Rigorous foundations	~10 hrs (ch. 1–7)	❌ (~$50)
StatQuest with Josh Starmer (YouTube)	Visual ML-stats connections	~6 hrs (prob/stats playlist)	✅
Bishop PRML Chapter 1–2	Probability for ML (advanced)	~8 hrs	✅ (PDF free)
Khan Academy Statistics	Fill gaps in fundamentals	~3 hrs	✅

Community note on StatQuest (r/learnmachinelearning, 2024–2026): Cited as frequently as Andrew Ng for clarifying statistical concepts. The videos on Bayes, MLE, and distributions are described as “the first time X actually made sense” by multiple practitioners. 20 hours for the full ML statistics playlist; individual videos are 8–20 minutes each. Use it when a concept isn’t clicking from a textbook.

---

Return to README.md · Next: 04_information_theory.md