Word Embeddings & Representations

import numpy as np
from collections import defaultdict

class MultinomialNB:
    def __init__(self, alpha=1.0):  # alpha: Laplace smoothing
        self.alpha = alpha
        self.class_log_prior_ = {}   # log P(Y=k)
        self.feature_log_prob_ = {}  # log P(xⱼ=1|Y=k)

    def fit(self, X, y):  # X: (N, V) count matrix; y: (N,)
        classes, counts = np.unique(y, return_counts=True)
        N, V = X.shape
        for k, cnt in zip(classes, counts):
            self.class_log_prior_[k] = np.log(cnt / N)
            X_k = X[y == k]              # rows belonging to class k
            # Smoothed MLE: (count + α) / (total + α·V)
            feature_counts = X_k.sum(axis=0) + self.alpha  # (V,)
            self.feature_log_prob_[k] = np.log(feature_counts / feature_counts.sum())
        return self

    def _log_likelihood(self, X, k):
        return X @ self.feature_log_prob_[k]  # dot product of counts & log-probs

    def predict_log_proba(self, X):  # (N, K) log-posteriors (unnormalised)
        classes = sorted(self.class_log_prior_)
        log_posts = np.stack(
            [self.class_log_prior_[k] + self._log_likelihood(X, k) for k in classes], axis=1)
        return log_posts - log_posts.max(axis=1, keepdims=True)

    def predict(self, X):
        return np.array(sorted(self.class_log_prior_))[self.predict_log_proba(X).argmax(axis=1)]

# ---- Spam/ham demo ----
from sklearn.datasets import fetch_20newsgroups
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.model_selection import train_test_split

data = fetch_20newsgroups(subset='all', categories=['sci.med', 'rec.sport.hockey'])
vec  = CountVectorizer(max_features=5000); X_bow = vec.fit_transform(data.data).toarray()
X_tr, X_te, y_tr, y_te = train_test_split(X_bow, data.target, test_size=0.2, random_state=42)

nb = MultinomialNB(alpha=1.0); nb.fit(X_tr, y_tr)
print(f"Accuracy: {(nb.predict(X_te)==y_te).mean():.3f}")
# Accuracy: 0.978  (98% on 5000-dim bag-of-words)

import numpy as np

class GaussianNB:
    def fit(self, X, y):
        self.classes_ = np.unique(y)
        self.theta_  = {}  # per-class mean vectors
        self.sigma_  = {}  # per-class variance vectors
        self.prior_  = {}
        N = len(y)
        for k in self.classes_:
            Xk = X[y == k]
            self.theta_[k] = Xk.mean(axis=0)
            self.sigma_[k] = Xk.var(axis=0) + 1e-9  # var + epsilon for stability
            self.prior_[k] = len(Xk) / N

    def _log_likelihood(self, X, k):
        mu, var = self.theta_[k], self.sigma_[k]
        # log N(x; mu, var) = -0.5*log(2πvar) - (x-mu)²/(2var)
        return -0.5 * np.sum(np.log(2 * np.pi * var) + (X - mu)**2 / var, axis=1)

    def predict(self, X):
        log_posts = np.stack(
            [np.log(self.prior_[k]) + self._log_likelihood(X, k) for k in self.classes_], axis=1)
        return self.classes_[log_posts.argmax(axis=1)]

from sklearn.datasets import load_iris
Xi, yi = load_iris(return_X_y=True)
gnb = GaussianNB(); gnb.fit(Xi[:120], yi[:120)
print(f"Iris accuracy: {(gnb.predict(Xi[120:])==yi[120:]).mean():.3f}")
# Iris accuracy: 0.967

Initialise: π, μ₁…μₖ, Σ₁…Σₖ  (random or k-means init)

repeat until convergence:

  E-step  (compute responsibilities)
    for n = 1…N, k = 1…K:
      γₙₖ  ←  πₖ · N(xₙ; μₖ, Σₖ)  /  Σⱼ πⱼ · N(xₙ; μⱼ, Σⱼ)
      # γₙₖ = soft probability that xₙ belongs to component k

  M-step  (update parameters using weighted MLE)
    Nₖ  ←  Σₙ γₙₖ              # effective count for component k
    πₖ  ←  Nₖ / N               # new mixing weight
    μₖ  ←  Σₙ γₙₖ xₙ / Nₖ      # weighted mean
    Σₖ  ←  Σₙ γₙₖ (xₙ-μₖ)(xₙ-μₖ)ᵀ / Nₖ  # weighted covariance

until | log P(X; θᵗ) - log P(X; θᵗ⁻¹) | < ε

import numpy as np
from scipy.stats import multivariate_normal

class GaussianMixtureModel:
    def __init__(self, K=3, max_iter=100, tol=1e-6, random_state=0):
        self.K, self.max_iter, self.tol = K, max_iter, tol
        self.rng = np.random.default_rng(random_state)

    def fit(self, X):
        N, D = X.shape
        K   = self.K

        # --- Initialise with k-means++ style ---
        idx      = self.rng.choice(N, K, replace=False)
        self.mu  = X[idx].copy()
        self.sig = [np.eye(D) for _ in range(K)]
        self.pi  = np.ones(K) / K

        prev_ll = -np.inf
        for it in range(self.max_iter):
            # --- E-step: compute responsibilities ---
            log_r = np.zeros((N, K))
            for k in range(K):
                log_r[:, k] = np.log(self.pi[k] + 1e-300) + \
                              multivariate_normal.logpdf(X, self.mu[k], self.sig[k])
            # Normalise with log-sum-exp for stability
            log_r -= log_r.max(axis=1, keepdims=True)
            r     = np.exp(log_r)
            r    /= r.sum(axis=1, keepdims=True)

            # --- M-step: update parameters ---
            Nk = r.sum(axis=0)
            self.pi  = Nk / N
            self.mu  = (r.T @ X) / Nk[:, None]
            for k in range(K):
                diff = X - self.mu[k]
                self.sig[k] = (r[:, k:k+1] * diff).T @ diff / Nk[k] + 1e-6*np.eye(D)

            # --- Check convergence via log-likelihood ---
            ll = self._log_likelihood(X)
            if abs(ll - prev_ll) < self.tol: break
            prev_ll = ll
        return self

    def _log_likelihood(self, X):  # sum log P(xₙ; θ) over all n
        log_pdfs = np.stack([np.log(self.pi[k] + 1e-300) +
                             multivariate_normal.logpdf(X, self.mu[k], self.sig[k])
                             for k in range(self.K)], axis=1)
        return np.logaddexp.reduce(log_pdfs, axis=1).sum()

    def predict(self, X):  # hard cluster assignment
        r = np.stack([self.pi[k] * multivariate_normal.pdf(X, self.mu[k], self.sig[k])
                      for k in range(self.K)], axis=1)
        return r.argmax(axis=1)

    def sample(self, n):  # generate new data from the learned model
        z = self.rng.choice(self.K, size=n, p=self.pi)
        return np.stack([self.rng.multivariate_normal(self.mu[k], self.sig[k]) for k in z])

import numpy as np
from sklearn.datasets import make_blobs

X_blobs, _ = make_blobs(n_samples=500, centers=4, cluster_std=0.8, random_state=7)

def bic(gmm, X):  # Bayesian Information Criterion: lower = better
    n_params = gmm.K * (X.shape[1] + X.shape[1)**2 + 1)
    return -2 * gmm._log_likelihood(X) + n_params * np.log(len(X))

def aic(gmm, X):  # Akaike IC: penalises less than BIC
    n_params = gmm.K * (X.shape[1] + X.shape[1)**2 + 1)
    return -2 * gmm._log_likelihood(X) + 2 * n_params

print(f"{'K':>4} {'LL':>10} {'BIC':>10} {'AIC':>10}")
for K in range(1, 8):
    gmm = GaussianMixtureModel(K=K).fit(X_blobs)
    ll  = gmm._log_likelihood(X_blobs)
    print(f"{K:>4} {ll:>10.1f} {bic(gmm,X_blobs):>10.1f} {aic(gmm,X_blobs):>10.1f}")
# K=4 has the best (lowest) BIC and AIC -- correctly identifies 4 clusters

import numpy as np

# Prove EM never decreases log-likelihood
gmm_track = GaussianMixtureModel(K=3, max_iter=1)  # 1 step at a time
gmm_track.fit(X_blobs)
ll_history = []

# Run 50 EM iterations and record log-likelihood
for step in range(50):
    gmm_track.fit(X_blobs)  # 1 more EM step (single iteration)
    ll_history.append(gmm_track._log_likelihood(X_blobs))

# Verify monotonic increase
diffs = np.diff(ll_history)
print(f"Min LL increase per step: {diffs.min():.6f}")
print(f"All steps non-decreasing: {(diffs >= -1e-9).all()}")
# Min LL increase: -0.000000  (never decreases)
# All steps non-decreasing: True

Input: observations x₁…xₜ, HMM parameters (π, A, B)
Output: P(X | λ) = P(x₁,…,xₜ | λ)

# Initialisation
for k = 1…K:
  α₁[k]  ←  π[k] · B[k](x₁)           # P(z₁=k, x₁)

# Recursion (dynamic programming, O(K²T))
for t = 2…T:
  for k = 1…K:
    αₜ[k]  ←  B[k](xₜ) · Σⱼ αₜ₋₁[j] · A[j,k]

# Termination
P(X|λ)  ←  Σₖ αₜ[k]

Input: observations x₁…xₜ, HMM parameters (π, A, B)
Output: z* = argmax_Z P(Z, X | λ)

# Initialisation
for k = 1…K:
  δ₁[k]  ←  log π[k] + log B[k](x₁)
  ψ₁[k]  ←  0                    # no predecessor at t=1

# Recursion (same DP as Forward, but max instead of sum)
for t = 2…T, k = 1…K:
  δₜ[k]  ←  log B[k](xₜ) + max_j [ δₜ₋₁[j] + log A[j,k] ]
  ψₜ[k]  ←  argmax_j [ δₜ₋₁[j] + log A[j,k] ] # backpointer

# Backtrack from best final state
z*ₜ  ←  argmax_k δₜ[k]
for t = T-1 … 1:  z*ₜ ← ψₜ₊₁[z*ₜ₊₁]

import numpy as np

class HMM:
    def __init__(self, pi, A, B):
        # pi: (K,)  A: (K,K)  B: (K,V) emission probs (discrete)
        self.pi = np.log(pi + 1e-300)
        self.A  = np.log(A  + 1e-300)
        self.B  = np.log(B  + 1e-300)
        self.K  = len(pi)

    def forward(self, obs):  # log-domain for stability
        T = len(obs); K = self.K
        alpha = np.full((T, K), -np.inf)
        alpha[0] = self.pi + self.B[:, obs[0])
        for t in range(1, T):
            for k in range(K):
                alpha[t, k] = np.logaddexp.reduce(alpha[t-1] + self.A[:, k]) + self.B[k, obs[t]]
        # log P(obs | HMM) = logsumexp of final alpha row
        return np.logaddexp.reduce(alpha[-1])  # scalar log-likelihood

    def viterbi(self, obs):  # most likely state sequence
        T = len(obs); K = self.K
        delta = np.full((T, K), -np.inf)
        psi   = np.zeros((T, K), dtype=int)
        delta[0] = self.pi + self.B[:, obs[0])
        for t in range(1, T):
            trans = delta[t-1:, None] + self.A  # (K, K)
            psi[t]   = trans.argmax(axis=0)
            delta[t] = trans.max(axis=0) + self.B[:, obs[t]]
        path = np.zeros(T, dtype=int)
        path[-1] = delta[-1].argmax()
        for t in range(T-2, -1, -1):  path[t] = psi[t+1, path[t+1]]
        return path

# ---- Demo: 2-state HMM (fair/loaded die) ----
pi = np.array([0.5, 0.5])
A  = np.array([[0.7, 0.3], [0.4, 0.6]])
B  = np.array([[1/6]*6, [0.1, 0.1, 0.1, 0.1, 0.1, 0.5]])  # biased towards 6

hmm  = HMM(pi, A, B)
obs  = [5, 5, 5, 1, 2, 5, 5]  # lots of 6s (0-indexed as 5)
print(f"Log P(obs): {hmm.forward(obs):.3f}")
print(f"State seq:  {hmm.viterbi(obs)}")
# Log P(obs): -10.82
# State seq:  [1 1 1 0 0 1 1]  -- mostly state 1 (loaded die) as expected

import numpy as np
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV

# Fit KDE on normal data; flag low-density points as anomalies
np.random.seed(0)
X_normal  = np.random.randn(500, 2)
X_anomaly = np.random.uniform(-5, 5, (20, 2))
X_all     = np.vstack([X_normal, X_anomaly])

# Select bandwidth via cross-validation
params = {'bandwidth': np.logspace(-1, 1, 10)}
grid   = GridSearchCV(KernelDensity(kernel='gaussian'), params, cv=5)
grid.fit(X_normal)
best_bw = grid.best_params_['bandwidth']
print(f"Best bandwidth: {best_bw:.3f}")

kde = KernelDensity(kernel='gaussian', bandwidth=best_bw).fit(X_normal)

# Log-density scores; low score = anomaly
log_density = kde.score_samples(X_all)  # (520,)
threshold   = np.percentile(log_density[:500], 5)  # 5th percentile of normal
flagged     = log_density < threshold

print(f"Normal points flagged: {flagged[:500].sum()} / 500 ({flagged[:500].mean():.1%})")
print(f"Anomalies detected:    {flagged[500:].sum()} / 20  ({flagged[500:].mean():.1%})")
# Normal points flagged: 25 / 500  (5.0%  -- calibrated as expected)
# Anomalies detected:    18 / 20   (90.0% -- correctly flagging outliers)

import numpy as np
from sklearn.naive_bayes import MultinomialNB as SKLearnNB
from sklearn.linear_model import LogisticRegression
from sklearn.datasets import fetch_20newsgroups
from sklearn.feature_extraction.text import CountVectorizer
from sklearn.model_selection import train_test_split

ng = fetch_20newsgroups(subset='all', categories=['sci.med', 'rec.sport.hockey'])
vec = CountVectorizer(max_features=5000)
X   = vec.fit_transform(ng.data).toarray()
y   = ng.target
X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=0.3, random_state=0)

print(f"{'N_train':>8} {'NB acc':>10} {'LR acc':>10}")
for n in [10, 30, 100, 300, 1000, len(X_tr)]:
    Xs = X_tr[:n]; ys = y_tr[:n]
    nb = SKLearnNB().fit(Xs, ys)
    lr = LogisticRegression(max_iter=1000).fit(Xs, ys)
    print(f"{n:>8} {nb.score(X_te,y_te):>10.3f} {lr.score(X_te,y_te):>10.3f}")
# N_train       NB acc     LR acc
#       10       0.697      0.612  <-- NB wins with tiny data
#       30       0.820      0.763
#      100       0.920      0.910
#      300       0.961      0.965  <-- LR overtakes around here
#     1000       0.972      0.985  <-- LR consistently better

Initialise: fit Naïve Bayes on labelled data only → θ₀

repeat until convergence:

  E-step: for each unlabelled xₙ:
    γₙₖ ← P(Y=k | xₙ; θ)
    # soft label assignment using current model

  M-step: refit NB using:
    Labelled data with hard labels
    Unlabelled data with soft labels γₙₖ
    # unlabelled data improves feature estimates

from sklearn.semi_supervised import SelfTrainingClassifier
from sklearn.naive_bayes import MultinomialNB as SKLearnNB
import numpy as np

# Simulate semi-supervised scenario: only 50 labelled, 1000 unlabelled
y_semi = y_tr.copy()
unlabelled_mask = np.ones(len(y_tr), dtype=bool)
unlabelled_mask[np.random.choice(len(y_tr), 50, replace=False)] = False
y_semi[unlabelled_mask] = -1  # -1 = unlabelled in sklearn

# Supervised only (50 labels)
sup_only = SKLearnNB().fit(X_tr[~unlabelled_mask], y_tr[~unlabelled_mask])
print(f"Supervised only  (50 labels): {sup_only.score(X_te, y_te):.3f}")

# Self-training: uses unlabelled data via iterative label propagation
st = SelfTrainingClassifier(SKLearnNB(), max_iter=10)
st.fit(X_tr, y_semi)  # y_semi has -1 for unlabelled
print(f"Semi-supervised (50+1000):    {st.score(X_te, y_te):.3f}")
# Supervised only  (50 labels): 0.812
# Semi-supervised (50+1000):    0.921  (unlabelled data helps significantly)