Advanced Probability

Bayes' theorem states P(A|B) = P(B|A) * P(A) / P(B). In probabilistic terms: posterior = likelihood * prior / evidence. The prior P(A) represents belief before observing data; the likelihood P(B|A) is the probability of data given the hypothesis; the posterior P(A|B) updates belief after observing data. In Bayesian statistics, the denominator P(B) = Σ P(B|A_i)P(A_i) normalizes the posterior. Applications: spam classification (P(spam|words)), medical diagnosis (P(disease|symptoms)), and sequential belief updating. Conjugate priors (Beta-Binomial, Gamma-Poisson) enable closed-form posterior computation without numerical integration.

# Beta-Binomial update
from scipy import stats
prior_alpha, prior_beta = 1, 1  # uniform prior
conversions, trials = 50, 200
posterior = stats.beta(prior_alpha + conversions, prior_beta + trials - conversions)
print(posterior.mean(), posterior.interval(0.95))

Law of Total Probability: for a partition {B_1,...,B_n} of the sample space, P(A) = Σ P(A|B_i) * P(B_i). It decomposes a complex probability by conditioning on mutually exclusive, exhaustive events. Law of Total Expectation (Law of Iterated Expectations): E[X] = E[E[X|Y]], which means we can compute the marginal expectation of X by first conditioning on Y and then taking expectation over Y. This is powerful for hierarchical models: compute conditional expectations within each group, then average. Extended form: Var(X) = E[Var(X|Y)] + Var(E[X|Y]) — the law of total variance, decomposing variance into within-group and between-group components.

# Law of total variance example
import numpy as np
# Var(X) = E[Var(X|Y)] + Var(E[X|Y])
mu_within = np.mean([np.var(group) for group in groups])
mu_between = np.var([np.mean(group) for group in groups])
total_var = mu_within + mu_between

The MGF of a random variable X is M_X(t) = E[e^{tX}], defined for t in a neighborhood of 0. The k-th derivative at t=0 gives the k-th raw moment: M^(k)(0) = E[X^k]. The MGF uniquely determines the distribution if it exists in a neighborhood of 0. Key property: if X and Y are independent, M_{X+Y}(t) = M_X(t) * M_Y(t), enabling derivation of sum distributions. For example, the sum of independent normals has MGF that is the product of individual MGFs, confirming the sum is normal. The characteristic function E[e^{itX}] (using imaginary unit) always exists and has the same properties without convergence issues.

from scipy import stats
# Normal MGF: M(t) = exp(mu*t + sigma^2*t^2/2)
import numpy as np
def normal_mgf(t, mu=0, sigma=1):
    return np.exp(mu*t + 0.5*sigma**2*t**2)

Markov's inequality states that for a non-negative random variable X and a > 0: P(X >= a) <= E[X] / a. It requires only that X is non-negative and E[X] exists — no other assumptions. It is tight (achieved with equality) when X takes value a with probability E[X]/a and 0 otherwise (a two-point distribution). The bound is often very loose in practice because it ignores distribution shape. Markov is the base for deriving Chebyshev's inequality (apply Markov to (X-μ)^2) and Chernoff bounds (apply Markov to the MGF). It is the weakest useful concentration inequality, requiring minimal assumptions.

# Markov's inequality: P(X >= a) <= E[X] / a
E_X = 5
a = 10
markov_bound = E_X / a  # = 0.5
print(f'P(X >= {a}) <= {markov_bound}')

Chebyshev's inequality states that for any random variable with finite mean μ and variance σ^2, P(|X - μ| >= k*σ) <= 1/k^2 for any k > 0. It applies to any distribution without assuming normality. For k=2: at most 25% of data lies more than 2 standard deviations from the mean; for k=3: at most 11%. Chebyshev is derived by applying Markov's inequality to (X-μ)^2. One-sided Cantelli's inequality: P(X - μ >= k*σ) <= 1/(1+k^2). Chebyshev is used for distribution-free confidence intervals and worst-case analysis, but is very loose compared to normal-based bounds when data is approximately normal.

def chebyshev_bound(k):
    return 1 / k**2
print(f'P(|X-μ| >= 2σ) <= {chebyshev_bound(2):.2f}')  # 0.25

Hoeffding's inequality provides tight exponential tail bounds for sums of bounded independent random variables. If X_1,...,X_n are independent with a_i <= X_i <= b_i, then P(|S_n - E[S_n]| >= t) <= 2*exp(-2t^2 / Σ(b_i-a_i)^2). It is much tighter than Chebyshev for bounded variables because it exploits the bounded range rather than just the variance. Applications: PAC learning theory (sample complexity bounds), bandit algorithm regret bounds, and confidence intervals for empirical means. For Bernoulli variables with range [0,1]: P(|X̄ - μ| >= ε) <= 2*exp(-2nε^2). It forms the basis of many learning theory proofs.

import numpy as np
def hoeffding_bound(n, epsilon, a=0, b=1):
    return 2 * np.exp(-2 * n * epsilon**2 / (b-a)**2)
print(f'P(|mean - mu| >= 0.1) <= {hoeffding_bound(1000, 0.1):.4f}')

The Central Limit Theorem (CLT) states that the standardized sum of n i.i.d. random variables with finite mean μ and variance σ^2 converges in distribution to a standard normal: √n(X̄_n - μ)/σ → N(0,1). Conditions: independence, identical distribution (relaxable via Lindeberg-Feller CLT for non-identical but small contributions), and finite variance. The rate of convergence depends on the third moment (Berry-Esseen theorem: error ≤ C*ρ/(σ^3 * √n), where ρ = E[|X|^3]). The CLT justifies normal approximations in A/B testing and confidence intervals for large samples. Failure: heavy-tailed distributions with infinite variance (Lévy-stable distributions) do not converge to normal.

import numpy as np
from scipy import stats
# Empirical CLT verification
sample_means = [np.mean(np.random.exponential(2, 100)) for _ in range(10000)]
stats.normaltest(sample_means)  # should be p > 0.05

The Beta distribution Beta(α, β) is defined on [0,1] with PDF proportional to x^{α-1}(1-x)^{β-1}. Mean = α/(α+β), Variance = αβ/((α+β)^2(α+β+1)). Special cases: Beta(1,1) = Uniform[0,1]. It is the conjugate prior for the Bernoulli and Binomial likelihood, making it the natural prior for proportions and conversion rates. α and β can be interpreted as pseudo-counts of successes and failures. As α, β increase, the distribution concentrates near its mean. Used in Bayesian A/B testing, Thompson Sampling for bandits, Bayesian neural networks (uncertainty quantification), and modeling proportional data in regression via beta regression.

from scipy import stats
import numpy as np
# Posterior after 30 successes in 100 trials (prior: Beta(1,1))
posterior = stats.beta(31, 71)
print(f'Mean: {posterior.mean():.3f}, 95% CI: {posterior.interval(0.95)}')

The Dirichlet distribution Dir(α_1,...,α_K) is the multivariate generalization of the Beta distribution, defined over probability simplices (vectors summing to 1). Its concentration parameter α determines how peaked or diffuse the distribution is: small α produces sparse distributions (near the corners), large α produces uniform-like distributions. It is the conjugate prior for the multinomial distribution. In NLP, it appears as the prior in Latent Dirichlet Allocation (LDA) topic models for both topic-per-document and word-per-topic distributions. Symmetric Dirichlet Dir(α,...,α) with small α (e.g., 0.1) encourages sparse topic mixtures, capturing that documents discuss few topics.

from scipy import stats
import numpy as np
# Sample from symmetric Dirichlet
alpha = np.full(5, 0.1)  # sparse topics
samples = stats.dirichlet(alpha).rvs(10)
print(samples.sum(axis=1))  # each sums to 1

The Poisson distribution models the number of events in a fixed interval when events occur at a constant average rate λ and independently. PMF: P(X=k) = λ^k * e^{-λ} / k!. Mean = Variance = λ. A Poisson process is a continuous-time stochastic process where events occur at constant rate λ: inter-arrival times are i.i.d. Exponential(λ), the number of events in [0,t] is Poisson(λt), and disjoint intervals are independent. The Poisson distribution approximates Binomial(n,p) when n is large and p is small (np = λ). Used for: server request rates, rare disease occurrences, customer arrivals, and A/B tests with count metrics.

from scipy import stats
lambda_rate = 10
dist = stats.poisson(lambda_rate)
print(f'P(X=8) = {dist.pmf(8):.4f}')
print(f'P(X<=12) = {dist.cdf(12):.4f}')

The Gamma distribution Gamma(α, β) has PDF ∝ x^{α-1}*e^{-x/β}. Mean = αβ, Variance = αβ^2. Special cases: Gamma(1, β) = Exponential(1/β); Gamma(n/2, 2) = Chi-squared(n). The sum of n i.i.d. Exponential(λ) variables follows Gamma(n, 1/λ). The Gamma is the conjugate prior for the Poisson rate and the precision (1/variance) of a normal. In Bayesian inference, the Normal-Gamma conjugate family jointly models the mean and precision. Gamma distributions model waiting times (time for k events in a Poisson process), service time distributions, and rainfall amounts.

from scipy import stats
# Time to 5th event in Poisson process with rate 2
k, lam = 5, 2
time_to_kth = stats.gamma(a=k, scale=1/lam)
print(f'Mean: {time_to_kth.mean():.2f}')

The Negative Binomial (NB) distribution models the number of failures before r successes in Bernoulli trials, or equivalently the number of trials needed for r successes. PMF: P(X=k) = C(k+r-1,r-1)*p^r*(1-p)^k. Mean = r(1-p)/p, Variance = r(1-p)/p^2. Since Variance > Mean when p < 1, NB is overdispersed relative to Poisson (which has Var = Mean). NB is widely used as a count distribution for overdispersed data: RNA-seq read counts, website traffic, customer purchase counts, and word frequencies. NB regression (glm.nb in R, statsmodels' NegativeBinomial) models overdispersed count outcomes where Poisson regression underestimates standard errors.

from scipy import stats
# NB as overdispersed Poisson alternative
nb = stats.nbinom(r=10, p=0.5)
print(f'Mean: {nb.mean()}, Var: {nb.var()}')  # Var > Mean

A martingale is a stochastic process {X_t} where E[X_{t+1} | X_1,...,X_t] = X_t — today's value is the best prediction of tomorrow's value given the full history. Submartingale: E[X_{t+1}|...] >= X_t (tends to increase). Supermartingale: <= X_t. The Optional Stopping Theorem (OST) states that if τ is a stopping time with E[τ] < ∞ and certain regularity conditions hold, then E[X_τ] = E[X_0]. This means expected wealth is preserved in a fair game regardless of when you stop — no gambling strategy can beat a fair game. Used in sequential testing, the SPRT likelihood ratio is a martingale, and e-values provide always-valid confidence sequences.

A Markov chain is a stochastic process where the next state depends only on the current state (Markov property). Transition matrix P has P_{ij} = P(X_{t+1}=j | X_t=i). A stationary distribution π satisfies π = πP: it is a left eigenvector of P with eigenvalue 1. The chain converges to π from any initial state if it is ergodic (irreducible — all states reachable from all others — and aperiodic). The convergence rate is determined by the second-largest eigenvalue. Stationary distributions are computed by solving πP = π, π*1 = 1. Applications: PageRank (web link graph Markov chain), MCMC sampling, queueing models, and natural language models.

import numpy as np
P = np.array([[0.9, 0.1], [0.2, 0.8]])
# Stationary: solve pi @ P = pi
from scipy.linalg import eig
vals, vecs = eig(P.T)
stat = vecs[:, np.argmax(vals.real)].real
stat /= stat.sum()
print(stat)

Detailed balance (reversibility) is a sufficient condition for a distribution π to be stationary for a Markov chain: π(x) * P(x→y) = π(y) * P(y→x) for all x, y. It means the probability flux from x to y equals that from y to x, ensuring equilibrium. In MCMC, the Metropolis-Hastings algorithm is designed to satisfy detailed balance with the target distribution π by accepting or rejecting proposed moves according to the MH acceptance ratio α = min(1, π(y)*q(x|y)/(π(x)*q(y|x))). Gibbs sampling satisfies detailed balance by sampling each variable from its full conditional. Detailed balance is sufficient but not necessary for stationarity — some efficient samplers violate it intentionally.

import numpy as np
def metropolis_step(x, log_target, proposal_std):
    y = x + np.random.normal(0, proposal_std)
    log_alpha = log_target(y) - log_target(x)
    if np.log(np.random.uniform()) < log_alpha:
        return y
    return x

The k-th order statistic X_{(k)} of a sample X_1,...,X_n is the k-th smallest value. The PDF of X_{(k)} from a continuous distribution with PDF f and CDF F is: f_{(k)}(x) = n!/(k-1)!(n-k)! * F(x)^{k-1} * (1-F(x))^{n-k} * f(x). Minimum = X_{(1)}, Maximum = X_{(n)}, Median = X_{((n+1)/2)}. Order statistics arise in: reliability theory (time to first failure = min), extreme value theory, distribution-free tests (Wilcoxon), bootstrap inference, and record-breaking probability. The expected maximum of n i.i.d. standard normals grows like √(2 log n) — important for multiple testing.

from scipy import stats
# Expected maximum of n uniform samples
n = 100
max_dist = stats.beta(n, 1)  # Maximum of n Uniform(0,1)
print(f'E[max of {n} uniforms]: {max_dist.mean():.4f}')

The Weak Law of Large Numbers (WLLN) states that the sample mean X̄_n converges in probability to the population mean μ as n → ∞: P(|X̄_n - μ| > ε) → 0 for any ε > 0. The Strong LLN (SLLN) states that X̄_n → μ almost surely (with probability 1). Both require finite mean; SLLN requires slightly stronger conditions. The LLN tells us the sample mean converges to the true mean but says nothing about the distribution of the error. The CLT goes further: it characterizes the distribution of the error √n(X̄_n - μ)/σ → N(0,1). LLN = convergence (first-order); CLT = distributional behavior of the error (second-order).

Convergence in probability: X_n →p X if P(|X_n - X| > ε) → 0 for all ε > 0. Convergence in distribution (weak convergence): X_n →d X if F_{X_n}(x) → F_X(x) at all continuity points of F_X. Convergence in probability is stronger: X_n →p X implies X_n →d X. The CLT proves convergence in distribution (to N(0,1)) but not in probability (the limiting random variable is not a fixed number). Slutsky's theorem: if X_n →d X and Y_n →p c (constant), then X_n + Y_n →d X + c and X_n * Y_n →d cX. This is used extensively in asymptotic statistics for constructing test statistics from estimated quantities.

The Chernoff bound is an exponentially tighter tail bound obtained by applying Markov's inequality to e^{tX} for the optimal t. For a sum S of n independent Bernoulli(p) variables: P(S >= (1+δ)np) <= (e^δ/(1+δ)^{1+δ})^{np}. This decays exponentially in the deviation, while Chebyshev decays only polynomially (1/k^2). Chernoff outperforms Chebyshev whenever variables are bounded and the desired tail probability is small. The bound requires: independence, bounded support, and the MGF must exist. Chernoff bounds are fundamental in: algorithm analysis (hash tables, load balancing), machine learning theory (PAC learning sample complexity), and randomized algorithms.

import numpy as np
def chernoff_bound_binomial(n, p, delta):
    mu = n * p
    return (np.exp(delta) / (1+delta)**(1+delta))**mu

The characteristic function (CF) of X is φ_X(t) = E[e^{itX}] (using imaginary unit i), the Fourier transform of the PDF. Unlike the MGF, the CF always exists for any distribution because |e^{itX}| = 1 — no convergence issues. The CF uniquely determines the distribution, enabling proof of the CLT via Lévy's continuity theorem. CF of the sum of independent variables is the product of individual CFs. The k-th cumulant is obtained from the cumulant generating function K(t) = log(MGF). CFs are essential for heavy-tailed distributions (Cauchy, Pareto) where the MGF does not exist due to heavy tails. The inversion theorem recovers the PDF from the CF.

import numpy as np
# Empirical characteristic function at t
def ecf(x, t):
    return np.mean(np.exp(1j * t * x))
print(ecf(np.random.normal(0,1,1000), t=1.0))

KL divergence measures how much probability distribution Q differs from P: KL(P||Q) = Σ P(x) log(P(x)/Q(x)). Properties: non-negative (KL >= 0 by Jensen's inequality, equals 0 iff P=Q), asymmetric (KL(P||Q) ≠ KL(Q||P)), not a metric (no triangle inequality). KL(P||Q) is the expected log-likelihood ratio under P — it appears in the cross-entropy: H(P,Q) = H(P) + KL(P||Q). Forward KL (KL(P||Q)) minimization leads to moment-matching; reverse KL (KL(Q||P)) minimization leads to mode-seeking. Used in variational inference, information-theoretic bounds, ELBO in VAEs, and policy gradient methods.

import numpy as np
def kl_divergence(p, q):
    p, q = np.asarray(p), np.asarray(q)
    return np.sum(p * np.log(p / q + 1e-15))

Shannon entropy H(X) = -Σ P(x) log P(x) measures the average uncertainty or information content of a random variable. It is maximized by the uniform distribution and minimized (= 0) by a deterministic variable. For a binary variable: H = -p*log2(p) - (1-p)*log2(1-p), peaking at 1 bit for p=0.5. Differential entropy h(X) = -∫f(x)log f(x)dx extends to continuous distributions but is not always positive. Mutual Information I(X;Y) = H(X) + H(Y) - H(X,Y) measures shared information between X and Y: zero for independent variables, equals H(X) when Y=X. Used in feature selection, decision trees (information gain), compression theory, and channel capacity.

import numpy as np
def entropy(p):
    p = np.asarray(p)
    p = p[p > 0]
    return -np.sum(p * np.log2(p))

def mutual_information(pxy):
    px = pxy.sum(axis=1, keepdims=True)
    py = pxy.sum(axis=0, keepdims=True)
    return np.sum(pxy * np.log2(pxy / (px * py) + 1e-15))

A copula is a joint distribution function on [0,1]^d that couples univariate marginals to form a multivariate distribution. Sklar's theorem: any joint CDF F(x,y) can be written as C(F_X(x), F_Y(y)) for a copula C. This separates the marginal distributions from the dependence structure. Common copulas: Gaussian (linear correlation), t-copula (tail dependence), Clayton (lower tail dependence), and Gumbel (upper tail dependence). Copulas enable modeling variables with different marginal distributions (e.g., normal and exponential) while capturing complex non-linear dependence. Used in financial risk modeling (CDO pricing), insurance, and extreme event co-occurrence analysis.

from scipy import stats
from scipy.stats import norm
# Gaussian copula sample
def gaussian_copula_sample(n, rho):
    cov = [[1, rho], [rho, 1]]
    z = np.random.multivariate_normal([0,0], cov, n)
    return norm.cdf(z)  # uniform marginals with Gaussian dependence

Inverse transform sampling generates samples from any continuous distribution F by: (1) sample U ~ Uniform(0,1), (2) return X = F^{-1}(U). This works because if U is uniform, F^{-1}(U) has CDF F. It requires a closed-form or numerically computable quantile function. For distributions without simple inverses (normal, gamma), alternatives include rejection sampling, Box-Muller transform (normal), or the Ziggurat algorithm. Inverse transform sampling is used for simulating non-standard distributions, generating correlated samples via copulas (transform through marginal CDFs), and implementing random variable generation in probabilistic programming systems.

import numpy as np
# Exponential samples via inverse transform
def sample_exponential(lam, n):
    U = np.random.uniform(0, 1, n)
    return -np.log(1 - U) / lam

samples = sample_exponential(2.0, 1000)

Rejection sampling generates samples from a target distribution p(x) by: (1) sampling a proposal x ~ q(x) and u ~ Uniform(0,1), (2) accepting x if u <= p(x)/(M*q(x)), where M = max p(x)/q(x) is the normalization constant. Accepted samples are distributed exactly according to p. Acceptance rate = 1/M, so choose q close to p for efficiency. Rejection sampling works for any distribution with a known unnormalized density and a bounding distribution q. Limitations: inefficiency in high dimensions (acceptance rate → 0 exponentially); requires M to be finite. For high-dimensional distributions, MCMC methods (Metropolis-Hastings, HMC) are preferred. Rejection sampling is used in probabilistic programming for small-dimensional posteriors.

import numpy as np
def rejection_sample(target_pdf, proposal_sampler, M, n):
    samples = []
    while len(samples) < n:
        x = proposal_sampler()
        u = np.random.uniform()
        if u <= target_pdf(x) / (M * proposal_pdf(x)):
            samples.append(x)
    return np.array(samples)

The Cauchy distribution has PDF f(x) = 1/(π*γ*(1+(x-x0)^2/γ^2)) — a ratio of two standard normals. It is heavy-tailed: tails decay as 1/x^2, just too slowly for the mean integral to converge. E[|X|] = ∞, so the mean is undefined. By the LLN, the sample mean of Cauchy observations does not converge to any value — it follows the same Cauchy distribution regardless of sample size. This makes standard statistical inference (t-tests, confidence intervals for means) completely invalid for Cauchy data. Robust statistics (median, IQR) are required. The Cauchy is the t-distribution with 1 degree of freedom and appears in physics as the Lorentzian line shape.

from scipy import stats
cauchy = stats.cauchy(loc=0, scale=1)
# Sample mean doesn't converge
for n in [100, 1000, 10000]:
    print(n, np.mean(cauchy.rvs(n)))  # stays Cauchy distributed

Importance sampling (IS) estimates E_p[f(X)] by sampling from a proposal q and weighting: E_p[f(X)] ≈ Σ w_i*f(x_i) / Σ w_i, where w_i = p(x_i)/q(x_i). IS works even when we cannot sample from p directly. Effective Sample Size (ESS) = (Σ w_i)^2 / Σ w_i^2 measures the equivalent number of i.i.d. samples from p. Low ESS (ESS << n) indicates the proposal q is a poor approximation of p, with a few high-weight samples dominating. Self-normalized IS avoids needing to know the normalizing constant. IS is used in Bayesian computation, sequential Monte Carlo (particle filters), and policy evaluation in reinforcement learning.

import numpy as np
def importance_sampling(f, log_target, log_proposal, samples):
    log_w = log_target(samples) - log_proposal(samples)
    log_w -= log_w.max()
    w = np.exp(log_w)
    w /= w.sum()
    ESS = 1 / np.sum(w**2)
    return np.dot(w, f(samples)), ESS

The Poisson limit theorem states that the Binomial(n, p) distribution converges to Poisson(λ) as n → ∞ and p → 0 with np = λ fixed. This approximation is accurate when n > 20 and p < 0.05. It justifies using the Poisson distribution for rare events: modeling defects in manufacturing, network packet arrivals, number of typos per page, or rare disease incidence. The bound on the approximation error is min(p, λ/n) * e * λ (via Chen-Stein method). The Poisson approximation also extends to sums of weakly dependent Bernoullis (local dependence), which is key for approximations in combinatorics, birthday problems, and random graph theory.

from scipy import stats
import numpy as np
n, p = 1000, 0.005  # lambda = 5
binom_pmf = stats.binom(n, p).pmf
poisson_pmf = stats.poisson(n*p).pmf
for k in range(10):
    print(f'k={k}: Binom={binom_pmf(k):.4f}, Pois={poisson_pmf(k):.4f}')

The birthday problem asks: how many people must be in a room so there is a > 50% probability that at least two share a birthday? The answer is 23. The probability of no collision among k people from n equally likely birthdays is n!/(n-k)! / n^k ≈ exp(-k(k-1)/(2n)) using the Poisson approximation. For 50% probability: k ≈ √(2n*ln 2) ≈ 1.177*√n. For n=365: k≈23. The birthday paradox reveals that collision probability grows as O(√n), not O(n). This has applications in hash function security (birthday attacks need O(√n) queries), cryptography, and approximate counting algorithms. It quantifies why collisions happen far sooner than intuition suggests.

import numpy as np
def birthday_prob(k, n=365):
    log_prob = sum(np.log((n-i)/n) for i in range(k))
    return 1 - np.exp(log_prob)
print(f'P(collision, k=23): {birthday_prob(23):.3f}')

A GMM models data as a mixture of K Gaussian distributions: p(x) = Σ_k π_k * N(x; μ_k, Σ_k), where π_k are mixing weights summing to 1. Introducing a latent variable Z ∈ {1,...,K} for the component: P(Z=k) = π_k and X|Z=k ~ N(μ_k, Σ_k). The EM algorithm alternates between E-step (compute posterior P(Z=k|X)) and M-step (update μ_k, Σ_k, π_k). GMMs are used for density estimation, soft clustering, outlier detection, and as priors in Bayesian models. The number of components K must be chosen via BIC, AIC, or cross-validation. GMMs generalize to mixture of experts and variational autoencoders.

from sklearn.mixture import GaussianMixture
gmm = GaussianMixture(n_components=3, covariance_type='full', random_state=42)
gmm.fit(X)
log_lik = gmm.score(X)
posterior = gmm.predict_proba(X)

The Bernstein-von Mises (BvM) theorem states that under regularity conditions, the posterior distribution converges to a normal distribution centered at the MLE with covariance equal to the inverse Fisher information matrix, as sample size n → ∞. Formally: √n(θ - θ_MLE) | data → N(0, I(θ_0)^{-1}). This provides the theoretical justification for asymptotic equivalence of Bayesian and frequentist inference for large samples — posterior credible intervals approximate frequentist confidence intervals. BvM fails in: infinite-dimensional models, misspecified models (true distribution not in the family), and certain nonparametric settings. It justifies the Laplace approximation for posterior computation.

The Cramér-Rao lower bound (CRLB) states that the variance of any unbiased estimator θ̂ of a parameter θ is bounded below by the inverse of the Fisher information: Var(θ̂) >= 1/I(θ), where I(θ) = E[(∂ log L/∂θ)^2]. An estimator achieving this bound is called efficient. The MLE is asymptotically efficient — it achieves the CRLB in large samples. The bound illustrates fundamental limits on estimation precision from data. It is used to evaluate estimators, guide experimental design (maximizing I(θ) via D-optimal design), and quantify the minimum variance achievable for any unbiased procedure.

from scipy import stats
import numpy as np
# Fisher info for Normal mean (sigma known)
def fisher_info_normal_mean(n, sigma):
    return n / sigma**2
CRLB = 1 / fisher_info_normal_mean(100, 1)
print(f'CRLB for mean: {CRLB:.4f}')

A statistic T(X) is sufficient for θ if the conditional distribution of X given T(X) does not depend on θ — T captures all information in the data about θ. The Fisher-Neyman factorization theorem: T(X) is sufficient for θ iff the likelihood factors as L(θ;x) = g(T(x),θ) * h(x), where h does not depend on θ. For exponential family distributions, the natural sufficient statistics are the ones that appear in the exponent. Sufficiency implies data reduction without loss of information: any inference should depend on data only through T. Minimal sufficient statistics capture all relevant information with minimum dimensionality.

The delta method approximates the distribution of a function g(X̄) of a sample mean X̄. If √n(X̄ - μ) →d N(0, σ^2), then by Taylor expansion: √n(g(X̄) - g(μ)) →d N(0, σ^2 * g'(μ)^2). This enables computing asymptotic standard errors for nonlinear transformations without simulation. Examples: variance of log(X̄) ≈ σ^2/(n*μ^2); variance of √X̄ ≈ σ^2/(4n*μ); variance of AUC estimates. The multivariate delta method uses the gradient vector: Var(g(X̄)) ≈ ∇g(μ)^T * Σ/n * ∇g(μ). Used for ratio metrics in A/B tests (delta method for CTR = clicks/impressions), survival analysis, and confidence intervals for nonlinear estimators.

import numpy as np
def delta_method_variance(n, mu, sigma, g_prime):
    """Asymptotic variance of g(Xbar)"""
    return sigma**2 * g_prime(mu)**2 / n

# Variance of log(Xbar)
mu, sigma_sq, n = 5, 4, 100
g_prime_log = lambda x: 1/x
print(delta_method_variance(n, mu, sigma_sq, g_prime_log))

The Gumbel distribution is one of three extreme value distributions (alongside Fréchet and Weibull) in the Fisher-Tippett-Gnedenko theorem: the maximum of a large number of i.i.d. variables converges in distribution to one of these three types, depending on the tail behavior of the original distribution. The Gumbel arises for thin-tailed distributions (normal, exponential). CDF: exp(-exp(-(x-μ)/β)). Mean = μ + β*γ (γ = Euler-Mascheroni constant ≈ 0.5772). The Gumbel distribution is used in: reliability engineering (max load analysis), records analysis, and the Gumbel-softmax (concrete distribution) in machine learning for differentiable discrete variable sampling.

from scipy import stats
gumbel = stats.gumbel_r(loc=0, scale=1)
print(f'Mean: {gumbel.mean():.4f}, Var: {gumbel.var():.4f}')

Cumulants are alternative descriptors of a distribution defined via the cumulant generating function (CGF): K(t) = log M(t) = log E[e^{tX}]. The first cumulant κ_1 = E[X] (mean), κ_2 = Var(X), κ_3 = E[(X-μ)^3] (third central moment = skewness * σ^3), κ_4 = excess kurtosis * σ^4. For the normal distribution, all cumulants of order >= 3 are zero — making cumulants a measure of non-normality. For sums of independent variables, cumulants add: κ_n(X+Y) = κ_n(X) + κ_n(Y). This additive property makes cumulants more convenient than moments for analyzing sums. Used in the Gram-Charlier and Edgeworth expansions for non-normal approximations.

from scipy import stats
import numpy as np
data = np.random.exponential(2, 1000)
print(f'Skewness: {stats.skew(data):.3f}')
print(f'Excess kurtosis: {stats.kurtosis(data):.3f}')

The Jeffreys prior is a non-informative prior proportional to the square root of the Fisher information determinant: π(θ) ∝ √det(I(θ)). It is invariant under reparameterization: if you transform the parameter, the Jeffreys prior on the new parameter is also the Jeffreys prior. This makes it a principled default for objective Bayesian inference. For a Binomial proportion: Jeffreys prior is Beta(0.5, 0.5) — concentrating at 0 and 1. For a Poisson rate: Jeffreys prior is Gamma(0.5, 0) — improper but proper posteriors result. The Jeffreys prior can be improper (does not integrate to 1) but still yields proper posteriors in regular models. It is the default choice when prior information is absent.

# Jeffreys prior for Binomial proportion: Beta(0.5, 0.5)
from scipy import stats
jeffreys = stats.beta(0.5, 0.5)
# Posterior after 30 successes in 100 trials
posterior = stats.beta(0.5+30, 0.5+70)
print(posterior.interval(0.95))

A Gaussian process is a distribution over functions: any finite collection of function values f(x_1),...,f(x_n) follows a multivariate normal distribution, fully specified by a mean function m(x) and covariance (kernel) function k(x,x'). The kernel encodes smoothness assumptions: RBF/squared exponential kernel gives infinitely differentiable functions; Matérn kernel gives functions with specified differentiability. GP regression (Kriging) provides a posterior over functions given observations, yielding predictions with uncertainty estimates. GPs are used in Bayesian optimization (GPyOpt, BoTorch), spatial statistics, emulation of expensive simulations, and as priors in non-parametric Bayesian models.

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel
kernel = ConstantKernel(1.0) * RBF(length_scale=1.0)
gpr = GaussianProcessRegressor(kernel=kernel, alpha=1e-6)
gpr.fit(X_train, y_train)
y_pred, y_std = gpr.predict(X_test, return_std=True)