Math for Machine Learning — Handbook

import numpy as np
def gaussian(x, mu, sigma):
    z = (x - mu) / sigma
    return np.exp(-0.5 * z**2) / (sigma * np.sqrt(2*np.pi))

x = np.linspace(-5, 5, 200)
pdf = gaussian(x, mu=0.0, sigma=1.0)
# fit mu, sigma from data by maximum likelihood = sample mean & std
data = np.random.normal(2.0, 1.5, size=10_000)
mu_hat, sigma_hat = data.mean(), data.std(ddof=0)
print(mu_hat, sigma_hat)        # ~2.0, ~1.5

# Bayes' rule for a medical test: P(disease | positive)
p_d   = 0.01          # prior prevalence
sens  = 0.99          # P(+ | disease)
spec  = 0.95          # P(- | no disease)
p_pos = sens*p_d + (1-spec)*(1-p_d)      # total prob of a positive
post  = sens*p_d / p_pos                  # posterior
print(round(post, 3))   # ~0.167 -- low despite a "99% accurate" test

import numpy as np
# Expectation and variance as probability-weighted sums
x = np.array([1, 2, 3, 4, 5, 6])         # a die
p = np.full(6, 1/6)
E  = (x * p).sum()                        # mean = 3.5
Var = ((x - E)**2 * p).sum()              # 2.9167
# Monte-Carlo check
rolls = np.random.randint(1, 7, size=1_000_000)
print(E, rolls.mean(), Var, rolls.var())

import numpy as np
# CLT: means of samples from ANY distribution become Gaussian
pop = np.random.exponential(scale=1.0, size=1_000_000)   # skewed
means = [np.random.choice(pop, size=n).mean()
         for _ in range(5000)]
import numpy as np; means = np.array(means)
# standard error shrinks like 1/sqrt(n)
for n in (1, 4, 16, 64):
    se = pop.std() / np.sqrt(n)
    print(n, round(se, 3))

import numpy as np
# Conditional probability from a joint table P(X, Y)
joint = np.array([[0.10, 0.20],
                  [0.30, 0.40]])          # rows X, cols Y
p_x = joint.sum(axis=1, keepdims=True)    # marginal P(X)
p_y_given_x = joint / p_x                  # rows sum to 1
print(p_y_given_x)
# independence test: joint == outer(P(X), P(Y)) ?
indep = np.allclose(joint, p_x * joint.sum(0))
print(indep)

import numpy as np
a = np.array([2.0, 1.0]); b = np.array([1.0, 3.0])
dot   = a @ b                              # 5.0
cos   = dot / (np.linalg.norm(a)*np.linalg.norm(b))
angle = np.degrees(np.arccos(cos))
proj  = (dot / (b @ b)) * b                 # projection of a onto b
print(dot, round(cos,3), round(angle,1))

import numpy as np
# A matrix is a linear map: it sends the unit square to a parallelogram
A = np.array([[2.0, 1.0],
              [0.0, 1.5]])
basis = np.eye(2)                          # columns are i-hat, j-hat
print(A @ basis)                            # where the basis vectors land
print("area scale =", np.linalg.det(A))     # determinant = area factor

import numpy as np
A = np.array([[2.0, 1.0],
              [1.0, 2.0]])
vals, vecs = np.linalg.eig(A)               # A v = lambda v
# verify
for i in range(2):
    print(np.allclose(A @ vecs[:, i], vals[i] * vecs[:, i]))
print("eigenvalues:", vals)                 # 3 and 1

import numpy as np
# Quadratic form x^T A x; positive-definite A => bowl with a unique min
A = np.array([[3.0, 1.0],
              [1.0, 2.0]])
def q(x): return x @ A @ x
print(np.all(np.linalg.eigvalsh(A) > 0))    # positive definite?
# gradient is 2 A x  ->  minimized at x = 0
x = np.array([1.0, -1.0]); print(q(x), 2*A@x)

import numpy as np
# Least-squares = orthogonal projection of y onto column space of X
X = np.array([[1, 0.], [1, 1.], [1, 2.], [1, 3.]])
y = np.array([1., 3., 4., 6.])
beta = np.linalg.solve(X.T @ X, X.T @ y)    # normal equations
yhat = X @ beta                              # projection
resid = y - yhat
print(beta, "resid . X =", resid @ X)        # residual is orthogonal to X

import numpy as np
# A function is convex if the chord lies above the curve everywhere
f = lambda x: x**2 + 3
def convex_check(f, a, b, t=0.5):
    mid = t*a + (1-t)*b
    return f(mid) <= t*f(a) + (1-t)*f(b)
print(all(convex_check(f, *np.random.randn(2)) for _ in range(1000)))

import numpy as np
# Gradient descent on a 2-D bowl  f(x) = x^T A x
A = np.array([[3., 0.], [0., 1.]])
grad = lambda x: 2 * A @ x
x, lr = np.array([4.0, 4.0]), 0.1
for _ in range(50):
    x = x - lr * grad(x)
print(x)            # -> [0, 0]; ill-conditioned A zig-zags

import numpy as np
# L2 regularization shrinks weights toward zero (ridge)
X = np.random.randn(50, 5); w_true = np.array([3,0,0,1,0.])
y = X @ w_true + 0.5*np.random.randn(50)
for lam in (0, 1, 10):
    w = np.linalg.solve(X.T@X + lam*np.eye(5), X.T@y)
    print(lam, np.round(w, 2))    # larger lam -> smaller weights

import numpy as np
# The learning rate: too small crawls, too large diverges
grad = lambda x: 2*x
for lr in (0.01, 0.5, 1.1):
    x, hist = 5.0, []
    for _ in range(20):
        x -= lr*grad(x); hist.append(x)
    print(f"lr={lr}: x->{x:.2e}")   # 1.1 explodes

import numpy as np
# An estimator is a function of data; study its sampling distribution
true_mu = 5.0
ests = [np.random.normal(true_mu, 2, size=30).mean() for _ in range(10000)]
ests = np.array(ests)
print("bias =", ests.mean() - true_mu)       # ~0 (unbiased)
print("std error =", ests.std())              # ~2/sqrt(30)

import numpy as np
# Sample variance: divide by n is biased; n-1 (Bessel) is unbiased
true_var = 4.0
biased, unbiased = [], []
for _ in range(20000):
    x = np.random.normal(0, 2, size=5)
    biased.append(x.var(ddof=0)); unbiased.append(x.var(ddof=1))
print(np.mean(biased), np.mean(unbiased), true_var)

import numpy as np
# 95% confidence interval for the mean (known idea: estimate +/- 1.96 SE)
x = np.random.normal(10, 3, size=40)
m, se = x.mean(), x.std(ddof=1)/np.sqrt(len(x))
lo, hi = m - 1.96*se, m + 1.96*se
print(round(lo,2), round(hi,2))   # ~95% of such intervals cover the truth

import numpy as np
from scipy import stats
# Two-sample test: is the mean different between groups?
a = np.random.normal(0.0, 1, 50)
b = np.random.normal(0.4, 1, 50)
t, p = stats.ttest_ind(a, b)
print(f"t={t:.2f}, p={p:.3f}", "reject H0" if p < 0.05 else "keep H0")

import numpy as np
# MLE vs MAP for a coin: MAP adds a Beta(a,b) prior (pseudo-counts)
heads, n = 7, 10
mle = heads / n                              # 0.70
a, b = 2, 2                                   # weak prior toward 0.5
mAP = (heads + a - 1) / (n + a + b - 2)       # 0.643
print(mle, round(mAP, 3))

import numpy as np
# Closed-form linear regression via the normal equations
X = np.c_[np.ones(100), np.random.randn(100, 2)]
beta_true = np.array([1.0, 2.0, -1.0])
y = X @ beta_true + 0.3*np.random.randn(100)
beta = np.linalg.lstsq(X, y, rcond=None)[0]
print(np.round(beta, 2))

import numpy as np
# Logistic regression by gradient descent on the log-loss
def sigmoid(z): return 1/(1+np.exp(-z))
X = np.c_[np.ones(200), np.random.randn(200, 2)]
y = (X[:,1] + X[:,2] > 0).astype(float)
w = np.zeros(3)
for _ in range(500):
    p = sigmoid(X @ w)
    w -= 0.1 * X.T @ (p - y) / len(y)         # gradient of log-loss
print(np.round(w, 2))

import numpy as np
# Lloyd's algorithm: assign to nearest centroid, then recompute centroids
X = np.vstack([np.random.randn(100,2),
               np.random.randn(100,2)+5])
k = 2; C = X[np.random.choice(len(X), k, replace=False)]
for _ in range(10):
    d = ((X[:,None,:]-C[None])**2).sum(-1)     # distances to centroids
    lab = d.argmin(1)
    C = np.array([X[lab==j].mean(0) for j in range(k)])
print(np.round(C, 2))

import numpy as np
# PCA = eigenvectors of the covariance (directions of max variance)
X = np.random.randn(500, 2) @ np.array([[3,1],[0,1.]])
Xc = X - X.mean(0)
cov = np.cov(Xc.T)
vals, vecs = np.linalg.eigh(cov)
order = vals.argsort()[::-1]
print("explained var ratio:", np.round(vals[order]/vals.sum(), 3))
Z = Xc @ vecs[:, order][:, :1]                 # project to top component

import numpy as np
# Margin = distance from the separating line to the nearest points
w, b = np.array([1.0, -1.0]), 0.0
margin = lambda x: (w @ x + b) / np.linalg.norm(w)
pts = np.array([[2, 0.], [0, 2.], [3, 1.]])
print(np.round([margin(p) for p in pts], 3))   # signed distances

import numpy as np
# Hinge loss: zero once a point is correct AND past the margin
def hinge(y, score): return np.maximum(0, 1 - y*score)
y = np.array([1, 1, -1, -1.])
score = np.array([2.0, 0.3, -1.5, 0.2])
print(hinge(y, score))     # [0, 0.7, 0, 1.2]

import numpy as np
# The kernel trick: inner products in feature space without the map
def phi(x): return np.array([x[0]**2, np.sqrt(2)*x[0]*x[1], x[1]**2])
a, b = np.array([1.,2.]), np.array([3.,1.])
print(phi(a) @ phi(b))           # explicit feature dot product
print((a @ b)**2)                # = polynomial kernel k(a,b)=(a.b)^2

import numpy as np
# RBF (Gaussian) kernel: similarity decays with distance
def rbf(A, B, gamma=0.5):
    d2 = ((A[:,None,:]-B[None])**2).sum(-1)
    return np.exp(-gamma * d2)
X = np.random.randn(5, 2)
K = rbf(X, X)                    # 5x5 Gram matrix, 1 on the diagonal
print(np.round(K, 2))

import numpy as np
# Best split = the threshold that most reduces impurity (Gini)
def gini(y): 
    p = np.bincount(y, minlength=2)/len(y); return 1 - (p**2).sum()
x = np.array([1,2,3,4,5,6.]); y = np.array([0,0,0,1,1,1])
best = min(((x[:-1]+x[1:])/2),
           key=lambda t: (gini(y[x<=t])*np.mean(x<=t) +
                          gini(y[x>t])*np.mean(x>t)))
print("best threshold:", best)   # ~3.5

from sklearn.tree import DecisionTreeClassifier
import numpy as np
X = np.random.randn(300, 2); y = (X[:,0]*X[:,1] > 0).astype(int)
tree = DecisionTreeClassifier(max_depth=3).fit(X, y)
print("train acc:", tree.score(X, y))
# deeper trees fit more -- and overfit; depth is the complexity knob

import numpy as np
# Bagging: average many trees trained on bootstrap samples -> less variance
from sklearn.tree import DecisionTreeRegressor
X = np.sort(np.random.rand(80,1),0); y = np.sin(6*X).ravel()+0.1*np.random.randn(80)
preds = []
for _ in range(50):
    idx = np.random.randint(0, len(X), len(X))     # bootstrap
    preds.append(DecisionTreeRegressor(max_depth=4).fit(X[idx], y[idx]).predict(X))
ensemble = np.mean(preds, 0)                          # smoother than any one tree

import numpy as np
# Gradient boosting: each tree fits the residual of the last
from sklearn.tree import DecisionTreeRegressor
X = np.sort(np.random.rand(120,1),0); y = np.sin(6*X).ravel()
pred = np.zeros_like(y); lr = 0.3
for _ in range(30):
    resid = y - pred
    pred += lr * DecisionTreeRegressor(max_depth=2).fit(X, resid).predict(X)
print("MSE:", np.mean((y-pred)**2))

import numpy as np
# Naive Bayes: assume features independent given the class
from sklearn.naive_bayes import GaussianNB
X = np.vstack([np.random.randn(100,2), np.random.randn(100,2)+3])
y = np.r_[np.zeros(100), np.ones(100)]
print("acc:", GaussianNB().fit(X, y).score(X, y))
# posterior  P(class|x)  ~  P(class) * prod_i P(x_i|class)

import numpy as np
from sklearn.mixture import GaussianMixture
X = np.vstack([np.random.randn(200,2), np.random.randn(200,2)+4])
gmm = GaussianMixture(n_components=2).fit(X)    # fit by EM
print("means:\n", np.round(gmm.means_, 2))
resp = gmm.predict_proba(X[:3])                  # soft cluster memberships
print(np.round(resp, 2))

import numpy as np
# Markov chain: next state depends only on the current state
P = np.array([[0.9, 0.1],
              [0.5, 0.5]])           # transition matrix
# stationary distribution = left eigenvector with eigenvalue 1
vals, vecs = np.linalg.eig(P.T)
pi = np.real(vecs[:, np.isclose(vals, 1)].ravel()); pi /= pi.sum()
print(np.round(pi, 3))

# Beta-Binomial conjugacy: posterior stays Beta, just add counts
a, b = 2, 2                     # prior Beta(2,2)
heads, tails = 8, 2
a_post, b_post = a + heads, b + tails
mean = a_post / (a_post + b_post)
print(a_post, b_post, round(mean, 3))   # closed-form Bayesian update

import numpy as np
# SVD factorizes ANY matrix: A = U S V^T; truncate for low-rank approx
A = np.random.randn(20, 10) @ np.random.randn(10, 8)   # rank <= 8
U, s, Vt = np.linalg.svd(A, full_matrices=False)
k = 3
A_k = U[:,:k] * s[:k] @ Vt[:k]            # best rank-3 approximation
print("rel error:", np.linalg.norm(A-A_k)/np.linalg.norm(A))

import numpy as np
# Whitening: rotate + scale so features have unit variance, no correlation
X = np.random.randn(1000,2) @ np.array([[2,1.],[0,1]])
Xc = X - X.mean(0)
cov = np.cov(Xc.T); vals, vecs = np.linalg.eigh(cov)
W = vecs @ np.diag(1/np.sqrt(vals)) @ vecs.T     # whitening matrix
Xw = Xc @ W
print(np.round(np.cov(Xw.T), 2))                  # ~identity

import numpy as np
from sklearn.manifold import TSNE
# t-SNE: preserve local neighborhoods for 2-D visualization
X = np.vstack([np.random.randn(100,10),
               np.random.randn(100,10)+5])
emb = TSNE(n_components=2, perplexity=30, init='pca').fit_transform(X)
print(emb.shape)        # (200, 2) -- two visible clusters

import numpy as np
# Random projection: a random matrix nearly preserves distances (JL lemma)
X = np.random.randn(500, 1000)
R = np.random.randn(1000, 50) / np.sqrt(50)
Z = X @ R                                  # 1000-D -> 50-D
# pairwise distances are roughly preserved
d_hi = np.linalg.norm(X[0]-X[1]); d_lo = np.linalg.norm(Z[0]-Z[1])
print(round(d_hi,1), round(d_lo,1))

import numpy as np
from sklearn.metrics import confusion_matrix
y_true = np.array([1,1,0,0,1,0,1,0])
y_pred = np.array([1,0,0,0,1,1,1,0])
tn, fp, fn, tp = confusion_matrix(y_true, y_pred).ravel()
precision = tp/(tp+fp); recall = tp/(tp+fn)
print(f"P={precision:.2f} R={recall:.2f} F1={2*precision*recall/(precision+recall):.2f}")

import numpy as np
from sklearn.metrics import roc_auc_score, roc_curve
y = np.r_[np.zeros(100), np.ones(100)]
scores = np.r_[np.random.rand(100), np.random.rand(100)+0.5]
fpr, tpr, thr = roc_curve(y, scores)
print("AUC =", round(roc_auc_score(y, scores), 3))  # ranking quality

import numpy as np
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LogisticRegression
X = np.random.randn(200, 4); y = (X[:,0] > 0).astype(int)
scores = cross_val_score(LogisticRegression(), X, y, cv=5)
print("fold scores:", np.round(scores, 3), "mean:", round(scores.mean(),3))

import numpy as np
# Calibration: among samples predicted p=0.7, ~70% should be positive
p = np.random.rand(10000)
y = (np.random.rand(10000) < p).astype(int)        # perfectly calibrated
bins = np.linspace(0,1,11)
for lo, hi in zip(bins[:-1], bins[1:]):
    m = (p>=lo)&(p<hi)
    if m.sum(): print(f"pred~{(lo+hi)/2:.1f}  actual={y[m].mean():.2f}")

import numpy as np
# User-based collaborative filtering via cosine similarity
R = np.array([[5,4,0,0],[5,5,0,0],[0,0,4,5],[0,0,5,4.]])  # users x items
def cos(u, v):
    m = (u>0)&(v>0)
    return 0 if m.sum()==0 else u[m]@v[m]/(np.linalg.norm(u[m])*np.linalg.norm(v[m]))
sim = np.array([[cos(R[i],R[j]) for j in range(4)] for i in range(4)])
print(np.round(sim, 2))

import numpy as np
# Matrix factorization: R ~ U V^T learned by gradient descent
R = np.array([[5,3,0,1],[4,0,0,1],[1,1,0,5],[0,0,5,4.]])
mask = R>0; U = np.random.rand(4,2); V = np.random.rand(4,2)
for _ in range(2000):
    E = mask*(R - U@V.T)
    U += 0.01*(E@V); V += 0.01*(E.T@U)
print(np.round(U@V.T, 1))           # fills in the blanks

import numpy as np
# Implicit feedback: clicks/plays as confidence, not ratings
plays = np.array([[0,3,0,12],[5,0,0,0],[0,0,8,1.]])
conf = 1 + 40*plays                  # more plays -> more confident "like"
pref = (plays > 0).astype(float)     # binary preference target
print(conf)      # weighting used by ALS for implicit data

import numpy as np
# Cold start: no history -> fall back to content features / popularity
item_feats = np.random.randn(5, 8)          # new items have features
user_profile = np.random.randn(8)            # from items they liked
scores = item_feats @ user_profile           # content-based ranking
print("recommend item:", scores.argmax())

import numpy as np
# Standardize: subtract mean, divide by std (fit on TRAIN only)
Xtr = np.random.randn(100, 3)*[10,1,0.1] + [5,0,-2]
mu, sd = Xtr.mean(0), Xtr.std(0)
Xtr_s = (Xtr - mu)/sd
Xte_s = (np.random.randn(20,3)*[10,1,0.1] + [5,0,-2] - mu)/sd   # reuse mu,sd
print(np.round(Xtr_s.mean(0),2), np.round(Xtr_s.std(0),2))

import numpy as np, pandas as pd
# One-hot encode a categorical column
df = pd.DataFrame({'color': ['red','blue','red','green']})
oh = pd.get_dummies(df['color'])
print(oh.values)        # each category -> its own 0/1 column

import numpy as np
from sklearn.preprocessing import PolynomialFeatures
# Polynomial features let a linear model fit curves
X = np.linspace(-2,2,5).reshape(-1,1)
Xp = PolynomialFeatures(degree=3, include_bias=False).fit_transform(X)
print(np.round(Xp, 2))   # columns: x, x^2, x^3

import numpy as np, pandas as pd
# Binning turns a continuous feature into ordinal buckets
age = np.array([5, 17, 22, 40, 67, 90])
bins = [0, 18, 35, 65, 120]
labels = ['child','young','adult','senior']
print(pd.cut(age, bins=bins, labels=labels).tolist())

import numpy as np
# Autocorrelation: correlation of a series with a lagged copy of itself
x = np.sin(np.linspace(0, 20, 200)) + 0.3*np.random.randn(200)
def acf(x, lag):
    x = x - x.mean()
    return (x[:-lag] @ x[lag:]) / (x @ x)
print([round(acf(x, k), 2) for k in (1, 5, 10)])

import numpy as np
# AR(1): x_t = phi * x_{t-1} + noise
phi, n = 0.8, 300
x = np.zeros(n)
for t in range(1, n):
    x[t] = phi*x[t-1] + np.random.randn()
# estimate phi back by regressing x_t on x_{t-1}
phi_hat = (x[:-1] @ x[1:]) / (x[:-1] @ x[:-1])
print(round(phi_hat, 2))

import numpy as np
# Decompose into trend + seasonal + residual
t = np.arange(120)
trend = 0.05*t; season = 2*np.sin(2*np.pi*t/12)
y = trend + season + 0.3*np.random.randn(120)
trend_hat = np.convolve(y, np.ones(12)/12, mode='same')   # moving average
season_hat = (y - trend_hat)
print(np.round(trend_hat[:3], 2))

import numpy as np
# Exponential smoothing forecast: weight recent points more
y = np.cumsum(np.random.randn(50)) + 10
alpha, level = 0.3, y[0]
for v in y:
    level = alpha*v + (1-alpha)*level    # update
forecast = level                          # next-step prediction
print(round(forecast, 2))

import numpy as np
# A confounder U drives both X and Y -> naive correlation is misleading
n = 5000
U = np.random.randn(n)
X = U + 0.5*np.random.randn(n)            # no real X->Y effect
Y = 2*U + 0.5*np.random.randn(n)          # Y depends on U, not X
naive = np.polyfit(X, Y, 1)[0]            # ~1, spurious
adjusted = np.polyfit(np.c_[X, U], Y, 1)  # control for U -> X coef ~0
print(round(naive, 2))

import numpy as np
# Randomization breaks confounding: treatment is independent of U
n = 10000; U = np.random.randn(n)
T = (np.random.rand(n) < 0.5).astype(float)   # randomly assigned
Y = 3*T + 2*U + np.random.randn(n)             # true effect of T is 3
ate = Y[T==1].mean() - Y[T==0].mean()          # difference in means
print(round(ate, 2))     # ~3, unbiased thanks to randomization

import numpy as np
# Backdoor adjustment: condition on confounders to block spurious paths
import pandas as pd
n=8000; Z=np.random.binomial(1,0.5,n)          # confounder
T=np.random.binomial(1,0.3+0.4*Z,n)            # Z->T
Y=1.5*T+2*Z+np.random.randn(n)                 # Z->Y, true T effect=1.5
df=pd.DataFrame({'Y':Y,'T':T,'Z':Z})
eff=sum(df[df.Z==z].groupby('T').Y.mean().diff().iloc[-1]*np.mean(Z==z)
        for z in (0,1))                         # adjust within strata of Z
print(round(eff,2))

import numpy as np
# Instrumental variable: Z affects Y only through X
n=10000; Z=np.random.randn(n); U=np.random.randn(n)
X=Z+U+0.3*np.random.randn(n)                   # Z->X (and confounder U)
Y=2*X+3*U+np.random.randn(n)                   # true X effect=2
iv = np.cov(Z,Y)[0,1] / np.cov(Z,X)[0,1]       # Cov(Z,Y)/Cov(Z,X)
print(round(iv, 2))     # ~2, unlike biased OLS

import numpy as np
# Posterior over a coin's bias on a grid (no conjugacy needed)
theta = np.linspace(0,1,201)
prior = np.ones_like(theta)                     # uniform prior
heads, tails = 8, 2
like = theta**heads * (1-theta)**tails
post = prior*like; post /= np.trapz(post, theta)
print("posterior mean:", round(np.trapz(theta*post, theta), 3))

import numpy as np
# Metropolis-Hastings samples from an unnormalized density
def logp(x): return -0.5*x**2                   # target ~ N(0,1)
x, samples = 0.0, []
for _ in range(20000):
    xp = x + np.random.normal(0, 1)             # proposal
    if np.log(np.random.rand()) < logp(xp)-logp(x):
        x = xp                                   # accept
    samples.append(x)
samples = np.array(samples[2000:])
print(round(samples.mean(),2), round(samples.std(),2))

import numpy as np
# Gibbs sampler for a correlated 2-D Gaussian (sample each coord in turn)
rho = 0.8; x = y = 0.0; S = []
for _ in range(20000):
    x = np.random.normal(rho*y, np.sqrt(1-rho**2))
    y = np.random.normal(rho*x, np.sqrt(1-rho**2))
    S.append((x, y))
S = np.array(S[2000:])
print("corr:", round(np.corrcoef(S.T)[0,1], 2))   # ~0.8

import numpy as np
# Posterior predictive: average predictions over the posterior, not a point
a_post, b_post = 10, 4                            # Beta posterior on bias
draws = np.random.beta(a_post, b_post, 10000)     # sample parameters
pred = np.random.binomial(1, draws)               # then sample data
print("P(next flip = heads) =", round(pred.mean(), 3))

import numpy as np
# Group fairness metrics from predictions and a protected attribute
y = np.array([1,0,1,1,0,1,0,0.])
pred = np.array([1,0,1,0,0,1,1,0.])
grp = np.array([0,0,0,0,1,1,1,1])           # protected group
for g in (0,1):
    m = grp==g
    print(f"group {g}: selection rate = {pred[m].mean():.2f}")  # demographic parity

import numpy as np
# Impossibility: you generally can't satisfy all fairness criteria at once
# unless base rates are equal. Show base rates differ -> tension.
base = {'A': 0.30, 'B': 0.60}     # P(y=1) per group
# Equalizing FPR and calibration simultaneously is impossible if base differs
print("base rates differ ->", base['A'] != base['B'],
      "=> calibration & equalized odds conflict")

import numpy as np
# Sampling bias: training group proportions != population -> skewed model
pop = np.r_[np.zeros(900), np.ones(100)]          # 10% group B in reality
sample = np.r_[np.zeros(500), np.ones(500)]        # 50% in (biased) sample
print("pop rate:", pop.mean(), "sample rate:", sample.mean())
# a model trained on `sample` will misjudge the population

import numpy as np
# Post-processing: pick per-group thresholds to equalize a chosen metric
scores = np.random.rand(1000); grp = np.random.randint(0,2,1000)
y = (scores + 0.2*grp + 0.2*np.random.randn(1000) > 0.6).astype(int)
thr = {g: np.quantile(scores[grp==g], 0.7) for g in (0,1)}   # group-specific
pred = np.array([scores[i] > thr[grp[i]] for i in range(1000)])
print({g: round(pred[grp==g].mean(),2) for g in (0,1)})

import numpy as np
# Kernel density estimate: sum of little Gaussians, one per data point
data = np.r_[np.random.randn(100), np.random.randn(50)+4]
def kde(x, data, h=0.4):
    return np.mean(np.exp(-0.5*((x-data)/h)**2)/(h*np.sqrt(2*np.pi)))
grid = np.linspace(-4, 8, 100)
dens = np.array([kde(x, data) for x in grid])
print("peak near:", round(grid[dens.argmax()], 1))

import numpy as np
# Robust z-score with the median absolute deviation (MAD)
x = np.r_[np.random.randn(100), [12, -9]]      # two outliers
med = np.median(x); mad = np.median(np.abs(x-med))
z = 0.6745*(x-med)/mad                          # robust score
print("flagged:", np.where(np.abs(z) > 3.5)[0])

import numpy as np
from sklearn.svm import OneClassSVM
# Learn a boundary around "normal" data; score new points as in/out
X = np.random.randn(200, 2)
oc = OneClassSVM(nu=0.05, gamma='scale').fit(X)
test = np.array([[0,0.],[6,6.]])
print(oc.predict(test))     # +1 = inlier, -1 = outlier

import numpy as np
from sklearn.ensemble import IsolationForest
# Isolation Forest: anomalies are isolated in fewer random splits
X = np.r_[np.random.randn(200,2), [[8,8]]]
iso = IsolationForest(contamination=0.01).fit(X)
print("anomaly?", iso.predict([[8,8],[0,0]]))   # -1 for the outlier

import numpy as np
# Online gradient descent: update on one example at a time
w = np.zeros(3)
for _ in range(5000):
    x = np.random.randn(3); y = (x @ [1.,-2,0.5] > 0)
    p = 1/(1+np.exp(-(w@x)))
    w -= 0.1 * (p - y) * x          # single-sample SGD step
print(np.round(w, 2))

import numpy as np
# Active learning: query the labels the model is most uncertain about
scores = np.random.rand(1000)                    # model P(y=1) on a pool
uncertainty = np.abs(scores - 0.5)               # near 0.5 = uncertain
to_label = uncertainty.argsort()[:10]            # cheapest informative 10
print("query indices:", to_label[:5])

import numpy as np
# Detect distribution drift by watching a rolling error rate
err = np.r_[np.random.rand(200) < 0.1,           # stable period
            np.random.rand(200) < 0.4]            # drift: errors jump
roll = np.convolve(err, np.ones(30)/30, 'valid')
alarm = np.argmax(roll > 0.25)                    # first crossing
print("drift detected near step", alarm)

import numpy as np
# Multi-armed bandit: epsilon-greedy, track regret vs the best arm
true = np.array([0.2, 0.5, 0.7]); Q = np.zeros(3); N = np.zeros(3)
regret = 0
for t in range(2000):
    a = np.random.randint(3) if np.random.rand()<0.1 else Q.argmax()
    r = np.random.rand() < true[a]
    N[a]+=1; Q[a]+=(r-Q[a])/N[a]
    regret += true.max() - true[a]
print("total regret:", round(regret, 1))

import numpy as np
# Sample functions from a GP prior (smooth random curves via the kernel)
X = np.linspace(0, 1, 100)
K = np.exp(-0.5*(X[:,None]-X[None])**2/0.05**2) + 1e-9*np.eye(100)
L = np.linalg.cholesky(K)
f = L @ np.random.randn(100)             # one prior sample
print(f.shape)            # draw several to see the prior over functions

import numpy as np
# GP regression: condition the prior on observed points -> mean + variance
def k(a,b,l=0.2): return np.exp(-0.5*(a[:,None]-b[None])**2/l**2)
Xtr=np.array([0.1,0.5,0.9]); ytr=np.sin(6*Xtr); Xte=np.linspace(0,1,50)
Ktt=k(Xtr,Xtr)+1e-6*np.eye(3); Kst=k(Xte,Xtr)
mu = Kst @ np.linalg.solve(Ktt, ytr)             # posterior mean
var = 1 - np.einsum('ij,jk,ik->i', Kst, np.linalg.inv(Ktt), Kst)
print(np.round(mu[:3], 2))

import numpy as np
# The kernel encodes assumptions: length-scale controls wiggliness
def rbf(X, l): return np.exp(-0.5*(X[:,None]-X[None])**2/l**2)
X = np.linspace(0,1,5)
print("short l (wiggly):\n", np.round(rbf(X,0.05),2))
print("long  l (smooth):\n", np.round(rbf(X,0.5),2))

import numpy as np
# Bayesian optimization: pick the next point by an acquisition function
# (here: upper confidence bound = mean + kappa * std from a GP posterior)
mu  = np.array([0.2, 0.5, 0.4, 0.55])     # GP posterior mean on candidates
std = np.array([0.3, 0.1, 0.05, 0.2])     # GP posterior std
ucb = mu + 2.0*std
print("evaluate candidate:", ucb.argmax())   # balances explore/exploit

import numpy as np
# Entropy in bits: expected surprise of a distribution
def entropy(p):
    p = p[p > 0]
    return -(p * np.log2(p)).sum()
print(entropy(np.array([0.5, 0.5])))      # 1.0 bit (max for 2 outcomes)
print(entropy(np.array([0.9, 0.1])))      # 0.469 bits (more predictable)

import numpy as np
# KL divergence: extra bits from using q to code data drawn from p
def kl(p, q): 
    p, q = np.asarray(p), np.asarray(q)
    return np.sum(p * np.log2(p/q))
p = np.array([0.7, 0.3]); q = np.array([0.5, 0.5])
print(round(kl(p,q),3), round(kl(q,p),3))   # asymmetric!

import numpy as np
# Mutual information from a joint distribution table
joint = np.array([[0.4, 0.1],
                  [0.1, 0.4]])              # correlated X,Y
px = joint.sum(1); py = joint.sum(0)
mi = sum(joint[i,j]*np.log2(joint[i,j]/(px[i]*py[j]))
         for i in range(2) for j in range(2) if joint[i,j]>0)
print("I(X;Y) =", round(mi, 3), "bits")

import numpy as np
# Optimal code length per symbol is about -log2 p; average >= entropy
p = np.array([0.5, 0.25, 0.125, 0.125])
code_len = -np.log2(p)                       # ideal bits per symbol
avg = (p * code_len).sum()
print("lengths:", code_len, "avg = entropy =", avg)

import numpy as np
# Random search explores more distinct values of the important knob
def objective(lr, _unused): return -(np.log10(lr)+3)**2   # peak near lr=1e-3
best = max(((10**np.random.uniform(-5,0), np.random.rand()) for _ in range(40)),
           key=lambda c: objective(*c))
print("best lr found:", f"{best[0]:.1e}")

import numpy as np
# Successive halving: many configs, little budget; keep the top half, repeat
configs = np.random.rand(16)            # "true" quality of each config
budget = 1
while len(configs) > 1:
    keep = np.argsort(configs)[::-1][:max(1,len(configs)//2)]
    configs = configs[keep]; budget *= 2
print("survivor quality:", round(configs[0], 2), "final budget x", budget)

import numpy as np
# Toy architecture search: hill-climb over (depth, width)
score = lambda d,w: -((d-5)**2 + (w-8)**2)      # unknown landscape
best = (np.random.randint(1,10), np.random.randint(1,16))
for _ in range(200):
    cand = (np.clip(best[0]+np.random.randint(-1,2),1,12),
            np.clip(best[1]+np.random.randint(-2,3),1,20))
    if score(*cand) > score(*best): best = cand
print("best architecture:", best)

import numpy as np
from sklearn.model_selection import validation_curve
from sklearn.tree import DecisionTreeClassifier
X = np.random.randn(300,5); y = (X[:,0]>0).astype(int)
tr, va = validation_curve(DecisionTreeClassifier(), X, y,
            param_name='max_depth', param_range=range(1,12), cv=5)
print("best depth:", list(range(1,12))[va.mean(1).argmax()])

import numpy as np
# Right-censored data: (time, event) where event=0 means "still alive at t"
time  = np.array([3.2, 5.1, 2.4, 6.8, 4.0])
event = np.array([1,   0,   1,   1,   0  ])   # 0 = censored
print("observed events:", event.sum(), "censored:", (event==0).sum())
# censored rows mean "true time > observed" -- never discard them

import numpy as np
# Kaplan-Meier estimator of the survival curve S(t)
time  = np.array([1,1,2,3,4,4,5]); event = np.array([1,1,1,0,1,1,1])
order = np.argsort(time); t, e = time[order], event[order]
S, surv, n = 1.0, [], len(t)
for i, (ti, ei) in enumerate(zip(t, e)):
    if ei == 1: S *= 1 - 1/(n-i)        # drop at each event
    surv.append((ti, round(S,3)))
print(surv)

import numpy as np
# Hazard -> survival: S(t) = exp(-cumulative hazard)
t = np.linspace(0, 5, 100)
hazard = 0.2 + 0.1*t                      # rising (wear-out) hazard
cum_haz = np.cumsum(hazard)*(t[1]-t[0])
S = np.exp(-cum_haz)
print("S(2.5) =", round(np.interp(2.5, t, S), 3))

import numpy as np
# Cox model: covariate scales a shared baseline hazard by exp(beta*x)
def survival(t, S0, beta, x): return S0(t) ** np.exp(beta*x)
S0 = lambda t: np.exp(-(t/3)**1.5)        # baseline survival
beta = 0.6
for x in (-1, 0, 1):
    print(f"x={x:+d}: S(2)={survival(2,S0,beta,x):.2f}  HR={np.exp(beta*x):.2f}")

import numpy as np
from sklearn.linear_model import LogisticRegression
# Self-training: label confident unlabeled points, add them, refit
Xl=np.array([[-2,0.],[2,0.]]); yl=np.array([0,1])
Xu=np.random.randn(200,2)+np.random.choice([-2,2],(200,1))
clf=LogisticRegression().fit(Xl,yl)
for _ in range(3):
    p=clf.predict_proba(Xu).max(1); conf=p>0.95
    Xl=np.vstack([Xl,Xu[conf]]); yl=np.r_[yl,clf.predict(Xu)[conf]]
    Xu=Xu[~conf]; clf.fit(Xl,yl)
print("labeled set grew to", len(yl))

import numpy as np
# Label propagation: diffuse labels through a similarity graph
X = np.r_[np.random.randn(20,2)-2, np.random.randn(20,2)+2]
f = np.zeros(40); f[0], f[20] = -1, 1            # two seed labels
W = np.exp(-0.5*((X[:,None]-X[None])**2).sum(-1))
for _ in range(50):
    f = W @ f / W.sum(1)
    f[0], f[20] = -1, 1                           # clamp seeds
print("cluster A mean label:", round(f[:20].mean(), 2))

import numpy as np
# Consistency loss: prediction should be stable under input perturbation
def model(x, w): return 1/(1+np.exp(-(x@w)))
x = np.random.randn(5, 3); w = np.random.randn(3)
p1 = model(x, w)
p2 = model(x + 0.05*np.random.randn(*x.shape), w)   # perturbed
print("consistency penalty:", round(np.mean((p1-p2)**2), 5))

import numpy as np
# Low-density separation: put the boundary in the gap between clusters
data = np.r_[np.random.randn(100)-2, np.random.randn(100)+2]
grid = np.linspace(-4, 4, 200)
dens = np.array([np.exp(-0.5*((g-data)/0.4)**2).mean() for g in grid])
boundary = grid[dens.argmin()]              # sparsest point
print("decision boundary near:", round(boundary, 2))