Math Foundations for ML

The map:

  Linear algebra  →  how data and models are REPRESENTED (vectors, matrices)
  Calculus        →  how models LEARN (gradients = which way to nudge)
  Probability     →  how models express UNCERTAINTY (distributions, likelihood)
  Statistics      →  how you know it WORKED (sampling, significance, baselines)

a · b = a₁b₁ + a₂b₂ + ... + aₙbₙ

[2, 1, 0] · [3, 1, 4] = 6 + 1 + 0 = 7

GRADIENT DESCENT (all of deep learning's training, one line):

    weights ← weights − learning_rate × gradient

"Take a small step DOWNHILL, repeat a few million times."

  - learning_rate = step size: too big → you leap over valleys and
    diverge; too small → training takes geologic time
    (the knob you met in fine-tuning)
  - backpropagation = the bookkeeping algorithm that computes the
    gradient efficiently through stacked layers (chain rule, applied
    relentlessly — the "how" of neural-networks training)

from typing import List, Tuple

def compute_loss(data: List[Tuple[float, float]], w: float, b: float) -> float:
    # Mean Squared Error (MSE)
    total_error = sum((y - (w * x + b)) ** 2 for x, y in data)
    return total_error / len(data)

def step_gradient(
    data: List[Tuple[float, float]], 
    w_current: float, 
    b_current: float, 
    learning_rate: float
) -> Tuple[float, float]:
    w_gradient = 0.0
    b_gradient = 0.0
    N = len(data)
    
    # Calculate slopes (derivatives of loss with respect to w and b)
    for x, y in data:
        prediction = w_current * x + b_current
        error = y - prediction
        w_gradient += -2 * x * error
        b_gradient += -2 * error
        
    # Nudge parameters opposite to gradient direction
    new_w = w_current - (w_gradient / N) * learning_rate
    new_b = b_current - (b_gradient / N) * learning_rate
    return new_w, new_b

# 1. Target line: y = 2x + 1 (with minor noise)
points = [(1.0, 3.1), (2.0, 4.9), (3.0, 7.2), (4.0, 8.8), (5.0, 11.1)]

# 2. Start from random initial weights and set learning rate
w, b = 0.0, 0.0
lr = 0.01

print("Training Loop:")
for epoch in range(1, 101):
    w, b = step_gradient(points, w, b, lr)
    if epoch % 20 == 0:
        loss = compute_loss(points, w, b)
        print(f"Epoch {epoch:03d} | Loss: {loss:.4f} | w: {w:.3f} | b: {b:.3f}")

Bayes' Theorem:

   P(A | B) = P(B | A) × P(A)
              ─────────────────
                    P(B)

Read: "The probability of A given B equals
       the probability of B given A, times the prior probability of A,
       divided by the probability of B."

Prior:              P(spam) = 0.10          (10% of emails are spam)
Likelihood:         P(word "lottery" | spam) = 0.80   (80% of spam contains "lottery")
P(word "lottery"):  P("lottery") = 0.05    (5% of all emails contain "lottery")

Posterior:
P(spam | "lottery") = (0.80 × 0.10) / 0.05 = 0.08 / 0.05 = 1.6 → normalized: high!

→ Seeing "lottery" updates our belief from 10% to 80% spam probability.

# Naive Bayes from scratch — text spam classification
def naive_bayes_predict(word_in_email: bool, p_spam=0.10, p_word_given_spam=0.80, p_word_given_ham=0.01) -> float:
    """
    Apply Bayes theorem to compute P(spam | word_present).
    """
    p_word = (p_word_given_spam * p_spam) + (p_word_given_ham * (1 - p_spam))
    p_spam_given_word = (p_word_given_spam * p_spam) / p_word
    return p_spam_given_word

print(f"P(spam | 'lottery' in email): {naive_bayes_predict(True):.3f}")  # ~0.615
print(f"P(spam | no 'lottery'): close to prior 10%")

Entropy H(S) = -Σ pᵢ × log₂(pᵢ)

Where pᵢ = proportion of class i in the dataset S.

Examples:
  Pure node:  100% cats     → H = -(1.0 × log₂(1.0)) = 0.0   (no uncertainty!)
  Coin flip:  50% cats, 50% dogs → H = -(0.5 × log₂(0.5) + 0.5 × log₂(0.5)) = 1.0
  3 classes equal: H = -(3 × 0.333 × log₂(0.333)) = 1.585

Information Gain = H(parent) - Σ (|child| / |parent|) × H(child)

Decision trees pick the feature and threshold that maximizes Information Gain.
→ This is the mathematical reason trees split the way they do.

import math

def entropy(labels: list) -> float:
    """Calculate Shannon entropy of a list of class labels."""
    if not labels:
        return 0.0
    total = len(labels)
    counts = {}
    for label in labels:
        counts[label] = counts.get(label, 0) + 1
    return -sum((count/total) * math.log2(count/total) for count in counts.values())

def information_gain(parent: list, left_child: list, right_child: list) -> float:
    """How much entropy is reduced by this split?"""
    n = len(parent)
    weighted_child_entropy = (len(left_child)/n) * entropy(left_child) + \
                              (len(right_child)/n) * entropy(right_child)
    return entropy(parent) - weighted_child_entropy

# Example: split customers (churns=1, stays=0)
parent = [1, 1, 0, 0, 0, 0, 1, 0]            # mixed — H = 1.0
left   = [1, 1, 1]                             # mostly churns
right  = [0, 0, 0, 0, 0]                       # mostly stays

print(f"Parent entropy:     {entropy(parent):.4f}")       # 1.0
print(f"Left entropy:       {entropy(left):.4f}")         # 0.918
print(f"Right entropy:      {entropy(right):.4f}")        # 0.0
print(f"Information Gain:   {information_gain(parent, left, right):.4f}")  # 0.549

import numpy as np

# --- Normal Distribution ---
mu, sigma = 0, 1
samples = np.random.normal(mu, sigma, 10000)
print(f"Mean: {samples.mean():.3f}, Std: {samples.std():.3f}")  # ~0, ~1

# --- Binomial Distribution ---
# 100 email, each has P(spam)=0.1 — how many spam?
spam_counts = np.random.binomial(n=100, p=0.1, size=5000)
print(f"Expected spam: {spam_counts.mean():.1f} (should be ≈10)")

# --- Why weights are initialized from Normal? ---
# He initialization for ReLU networks: mean=0, std=sqrt(2/n_inputs)
n_inputs = 512
weights = np.random.normal(0, np.sqrt(2/n_inputs), (512, 256))
print(f"Weight std: {weights.std():.4f}")  # ~0.0625

CLT states:
  The distribution of SAMPLE MEANS approaches a Normal distribution
  as the sample size increases — regardless of the original distribution.

  If you take many samples of size n from ANY distribution:
    Mean of samples → N(μ, σ²/n)    (Normal with same mean, smaller variance)

Why it matters:
  - You can compute confidence intervals on any metric
  - You can run hypothesis tests even on non-normal data
  - Error bars on ML evaluation metrics are valid

import numpy as np
import random

# Original distribution: wildly skewed exponential (NOT normal)
population = np.random.exponential(scale=2.0, size=100_000)
print(f"Population mean: {population.mean():.3f}, Skew: very right-skewed")

# Take 10,000 samples of size 50 — compute each sample mean
sample_means = []
for _ in range(10_000):
    sample = np.random.choice(population, size=50, replace=False)
    sample_means.append(sample.mean())

sample_means = np.array(sample_means)
print(f"Sample means mean: {sample_means.mean():.3f}")  # ≈ 2.0 (population mean)
print(f"Sample means std:  {sample_means.std():.4f}")   # ≈ σ/√50
# → The distribution of sample means IS approximately Normal — that's CLT!

The Process:
  1. Null Hypothesis (H₀): "The new model is no better than the old model"
  2. Alternative (H₁): "The new model IS better"
  3. Collect data: run experiment, collect results
  4. Compute p-value: probability of observing this data IF H₀ is true
  5. Decision:
       p-value < 0.05 → reject H₀ → statistically significant improvement
       p-value ≥ 0.05 → fail to reject H₀ → not enough evidence

p-value intuition:
  p = 0.03 means: "If the model was actually not better,
  there's only a 3% chance we'd see this big of a difference by luck."
  → Low p → the improvement is unlikely to be random → ship it.

Type I Error (False Positive):  Conclude model is better when it isn't
                                  → Significance level α = 0.05 = allowed Type I rate

Type II Error (False Negative): Conclude model is no better when it actually is
                                  → Power = 1 - β = probability of detecting a real effect

The trade-off:
  Tighter threshold (α=0.01) → fewer Type I errors, more Type II errors
  Loose threshold (α=0.10) → more Type I errors, fewer Type II errors

from scipy import stats
import numpy as np

# --- T-test: comparing two model's accuracy scores across 5 runs ---
model_a_scores = [0.842, 0.851, 0.838, 0.845, 0.849]
model_b_scores = [0.861, 0.874, 0.869, 0.878, 0.865]

t_stat, p_value = stats.ttest_ind(model_a_scores, model_b_scores)
print(f"t-statistic: {t_stat:.3f}")
print(f"p-value: {p_value:.4f}")

if p_value < 0.05:
    print("✅ Model B is significantly better (p < 0.05)")
else:
    print("❌ No statistically significant difference — don't ship on hype")

Eigenvector intuition:
  Most data "stretches" in some preferred directions.
  For a dataset of heights and weights, the main stretch direction is
  "both height and weight increase together" (they're correlated).

  An eigenvector is that "preferred stretch direction."
  Its eigenvalue measures HOW MUCH the data stretches in that direction.

PCA (Principal Component Analysis):
  Step 1: Center the data (subtract mean)
  Step 2: Compute the covariance matrix
  Step 3: Find eigenvectors (the principal components)
  Step 4: Sort by eigenvalue — larger eigenvalue = more variance explained
  Step 5: Project data onto the top k eigenvectors → compressed representation

import numpy as np

# --- PCA from scratch ---
def pca_from_scratch(X: np.ndarray, n_components: int) -> tuple:
    """
    Principal Component Analysis without sklearn.
    X: (n_samples, n_features)
    Returns: (transformed data, explained_variance_ratio)
    """
    # 1. Center the data
    X_centered = X - X.mean(axis=0)
    
    # 2. Compute covariance matrix
    cov_matrix = np.cov(X_centered, rowvar=False)
    
    # 3. Eigendecomposition — find principal directions
    eigenvalues, eigenvectors = np.linalg.eigh(cov_matrix)
    
    # 4. Sort by eigenvalue (descending — largest variance first)
    sorted_idx = np.argsort(eigenvalues)[::-1]
    eigenvalues = eigenvalues[sorted_idx]
    eigenvectors = eigenvectors[:, sorted_idx]
    
    # 5. Select top k components
    components = eigenvectors[:, :n_components]
    
    # 6. Project data
    X_reduced = X_centered @ components
    
    # 7. Explained variance ratio
    explained_variance_ratio = eigenvalues[:n_components] / eigenvalues.sum()
    
    return X_reduced, explained_variance_ratio

# Create sample data: 2D clusters (like height/weight in 4D feature space)
np.random.seed(42)
X = np.random.randn(200, 4)  # 200 samples, 4 features
X[:, 1] = X[:, 0] * 0.8 + np.random.randn(200) * 0.2  # make feature 0 and 1 correlated

X_reduced, evr = pca_from_scratch(X, n_components=2)
print(f"Original shape: {X.shape}")          # (200, 4)
print(f"Reduced shape:  {X_reduced.shape}")  # (200, 2)
print(f"Variance explained by PC1: {evr[0]:.1%}")
print(f"Variance explained by PC2: {evr[1]:.1%}")
print(f"Total variance explained:  {evr.sum():.1%}")

# Verify using sklearn
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X_sk = pca.fit_transform(X)
print(f"\nsklearn PCA variance ratios: {pca.explained_variance_ratio_}")