Bayesian methods

Bayes' theorem states P(H|E) = P(E|H)·P(H) / P(E). The four components are: the posterior P(H|E) — probability of hypothesis given evidence; the likelihood P(E|H) — probability of observing evidence if the hypothesis is true; the prior P(H) — initial belief in the hypothesis before data; and the evidence P(E) — total probability of observing the data. It allows us to rationally update beliefs as new evidence arrives.

The prior distribution encodes our belief about a parameter before seeing any data; the posterior distribution reflects our updated belief after incorporating the observed data via the likelihood. The posterior combines prior and likelihood through Bayes' theorem: P(θ|data) ∝ P(data|θ)·P(θ). As we gather more data, the posterior becomes less influenced by the prior and more driven by the evidence, eventually converging toward the true parameter value.

The likelihood function L(θ|x) measures how probable the observed data x is for different values of the parameter θ; it is a function of θ, not of x. Unlike a probability distribution, likelihood values do not sum or integrate to 1 over θ. In Bayes' theorem, the likelihood P(data|θ) is the bridge that links prior belief about θ to the posterior belief after observing data.

A conjugate prior is a prior distribution that, when combined with a specific likelihood, produces a posterior in the same distributional family. For the Binomial likelihood (counting k successes in n trials), if the prior on the success probability θ is Beta(α, β), the posterior is Beta(α+k, β+n−k) — still a Beta distribution. This mathematical convenience makes computations tractable: α and β act as pseudo-counts of prior successes and failures.

Frequentists interpret probability as the long-run frequency of an event over infinitely repeated experiments; parameters are fixed constants and data is random. Bayesians treat probability as a degree of belief, so parameters themselves can have probability distributions. This means Bayesian methods can incorporate prior knowledge and produce direct probability statements about parameters (e.g., 'there is a 90% probability θ > 0.5'), which frequentist methods cannot legitimately make.

After observing k=7 heads and n−k=3 tails with prior Beta(5,5), the posterior is Beta(5+7, 5+3) = Beta(12,8). The posterior mean is 12/(12+8) = 0.60, which lies between the MLE of 0.70 and the prior mean of 0.50. The prior acts as regularization, pulling the estimate toward 0.5. With more data the influence of the prior diminishes, and the posterior mean approaches the true bias.

MAP estimation finds the parameter value that maximizes the posterior P(θ|data) ∝ P(data|θ)·P(θ), effectively incorporating the prior. MLE ignores the prior and only maximizes the likelihood P(data|θ). MAP can be seen as regularized MLE: a Gaussian prior on parameters corresponds to L2 (Ridge) regularization. When the prior is uniform over all values, MAP and MLE give identical results.

The naive assumption is that all features are conditionally independent given the class label: P(x₁,…,xₙ|C) = ∏ P(xᵢ|C). In practice, features are often correlated (e.g., word co-occurrences in text), so this assumption is frequently wrong. Despite this, Naive Bayes performs well in many real-world tasks like spam detection because the classifier only needs to rank classes correctly, not estimate probabilities precisely. It is also extremely fast and works well with limited training data.

In sequential Bayesian updating the posterior from one observation becomes the prior for the next. After the first positive: unnormalized scores are 0.9×0.5=0.45 (works) and 0.2×0.5=0.10 (fails), giving P(works|+₁) ≈ 0.818. This becomes the new prior and the process repeats for the second and third positives, rapidly increasing confidence. The final posterior is identical whether you process all three results at once or one at a time, because multiplication is commutative — Bayesian updating is order-independent.

The posterior predictive distribution P(x̃|x) = ∫ P(x̃|θ)P(θ|x) dθ averages predictions over all possible parameter values weighted by their posterior probability. This accounts for parameter uncertainty: instead of predicting as if θ is known exactly, we integrate out our uncertainty. A point estimate such as MAP or MLE ignores this uncertainty and can be overconfident, producing intervals that are too narrow. The posterior predictive naturally gives wider and more honest uncertainty bounds on future observations.

In most real problems the posterior P(θ|data) is analytically intractable because the normalizing constant (marginal likelihood P(data)) requires integrating over all possible parameter values — a high-dimensional integral that is infeasible. MCMC constructs a Markov chain whose stationary distribution equals the target posterior, generating samples without ever computing the normalizing constant. These samples allow estimation of posterior means, credible intervals, and other quantities numerically. Common MCMC algorithms include Metropolis-Hastings and Gibbs sampling.

The acceptance probability is α = min(1, [P(θ*|data)·q(θ|θ*)] / [P(θ|data)·q(θ*|θ)]), where θ* is the proposed state and q is the proposal distribution. If the proposed state has higher posterior density it is always accepted; if it has lower density, it is accepted with probability equal to the posterior density ratio. This 'detailed balance' condition guarantees the Markov chain's stationary distribution equals the target posterior. Crucially, the normalizing constant cancels in the ratio so it never needs to be computed.

Gibbs sampling updates one parameter at a time by sampling from its full conditional distribution P(θᵢ|θ_{-i}, data) — the distribution of θᵢ given all other parameters and data. Unlike Metropolis-Hastings, Gibbs proposals are always accepted (acceptance rate = 1) because samples are drawn exactly from the conditional. It is preferred when these conditional distributions are known and easy to sample from, as in models with conjugate priors. However, Gibbs can be slow when parameters are highly correlated, since it takes many small steps to explore the joint posterior.

A 95% Bayesian credible interval [a, b] means: given the observed data, there is a 95% probability that the true parameter θ lies in [a, b]. A frequentist 95% confidence interval means: if we repeated the experiment infinitely many times, 95% of the constructed intervals would contain the true θ. The Bayesian interval makes a direct probability statement about the parameter, which is what most people intuitively want. The frequentist interval treats θ as fixed and can never attach a probability to it being in a specific interval.

The Bayes factor BF₁₂ = P(data|M₁) / P(data|M₂) compares the marginal likelihood of two models — how well each model predicts the observed data averaged over its prior on parameters. A BF₁₂ > 10 is generally considered strong evidence for M₁. Unlike p-values, Bayes factors can also provide evidence for a simpler (null) model when BF₁₂ < 1. They automatically penalize model complexity because overly complex models spread their prior over many parameter values, reducing the average likelihood.

The classifier computes P(spam|words) ∝ P(spam)·∏ P(wordᵢ|spam) using word frequency counts from labelled training emails as the likelihoods. Words like 'FREE' or 'WINNER' have high P(word|spam) and push the posterior toward spam. If a word never appeared in training spam emails, P(word|spam) = 0, which collapses the entire product to zero regardless of all other words — this is the zero-probability problem, which is fixed by Laplace smoothing (adding a pseudocount of 1 to every word's count).

Laplace smoothing adds a pseudocount of 1 to every word's frequency before computing probabilities: P(wordᵢ|class) = (count(wordᵢ,class)+1)/(total words in class + vocabulary size). Without smoothing, any word absent from training gets P=0, zeroing out the entire product (the zero-probability problem). Smoothing ensures every word gets a small nonzero probability. Mathematically it is equivalent to placing a Dirichlet prior with all pseudo-counts equal to 1 on the word probability vector.

An informative prior encodes genuine domain knowledge about a parameter before seeing data — for example, Beta(10,10) encodes strong belief that a coin is fair. A non-informative or flat prior assigns roughly equal probability to all values (e.g., Beta(1,1) = Uniform[0,1]) and lets the data dominate. Flat priors appear objective but can still carry implicit assumptions under reparametrisation, a problem addressed by Jeffreys priors. Using weakly informative priors is often a good practical compromise, providing mild regularisation without overwhelming the data.

OLS produces a single best-fit weight vector w by minimising squared error, yielding a point estimate. Bayesian linear regression instead places a prior distribution over w (often Gaussian), combines it with the Gaussian likelihood of the data via Bayes' theorem, and computes a full posterior distribution over w. Every prediction comes with a credible interval reflecting both data noise and parameter uncertainty, which is especially valuable when data is scarce or when honest uncertainty quantification matters for decision-making.

A hierarchical (multilevel) model places priors on hyperparameters that govern the priors for individual-level parameters, so groups borrow statistical strength from each other. For example, to estimate click-through rates for 50 websites, instead of treating each independently (overfitting small groups) or pooling all data (ignoring group differences), a hierarchical model uses a shared Beta hyperprior so data-sparse sites are pulled toward the global average. This shrinkage avoids both overfitting and underfitting and is ideal when group sizes vary greatly.

In Bayesian A/B testing each variant's conversion rate is modelled with a Beta prior, updated to a posterior from observed conversions. We can directly compute P(variant B rate > variant A rate) by comparing posterior samples — a concrete probability statement. Frequentist A/B testing computes a p-value under the null hypothesis of no difference, but cannot state 'there is an 85% chance B is better'. Bayesian tests also allow principled early stopping: once P(B > A) crosses a threshold the test can be halted without inflating error rates the same way optional stopping does in frequentist testing.

Let D = disease, T+ = positive test. Given P(D)=0.01, P(T+|D)=0.99, and P(T+|no D)=0.05 (1-specificity). By total probability: P(T+) = 0.99×0.01 + 0.05×0.99 = 0.0099+0.0495 = 0.0594. Bayes' theorem gives P(D|T+) = 0.0099/0.0594 ≈ 0.167, or only about 17%. Despite the highly sensitive test, the low disease prevalence means most positives are false positives — this base-rate fallacy shows why Bayesian reasoning is critical in medical diagnosis.

The mean of Beta(α,β) is α/(α+β) and the variance is αβ/[(α+β)²(α+β+1)]. When α=β the distribution is symmetric around 0.5; when α>β it skews toward 1, and when α<β it skews toward 0. Larger α+β (more total pseudo-counts) makes the distribution narrower, reflecting higher certainty. For example, Beta(6, 4) has mean 6/(6+4) = 0.60 and is skewed toward 1 since α > β. Special cases: Beta(1,1) is the uniform distribution on [0,1], and Beta(α,β) with both α,β < 1 is U-shaped (bimodal at 0 and 1).

A Bayesian network is a directed acyclic graph (DAG) where nodes represent random variables and edges represent direct probabilistic dependencies. The joint distribution factorises as P(X₁,…,Xₙ) = ∏ P(Xᵢ | parents(Xᵢ)). A node is conditionally independent of all non-descendants given its parents, which drastically reduces the number of parameters required. For example, a network 'Rain → Wet Grass ← Sprinkler' encodes that rain and the sprinkler each independently cause wet grass.

The Jeffreys prior is P(θ) ∝ √det(I(θ)), where I(θ) is the Fisher information matrix. Its key property is invariance under reparametrisation: transforming θ→φ(θ) and applying Jeffreys to φ gives the same result as transforming the Jeffreys prior on θ. A flat uniform prior is NOT invariant — flat in one parametrisation is not flat in another (e.g., flat on p is not flat on log-odds of p). This invariance makes the Jeffreys prior a principled 'objective' choice for representing complete prior ignorance, though it can be improper (not integrating to 1) in some models.

The Beta distribution is the conjugate prior for Binomial data, so after observing k=7 heads in n=10 flips the posterior is Beta(α+k, β+n−k) = Beta(9,5). The prior mean is 0.50 and the posterior mean is 9/(9+5) ≈ 0.643, demonstrating how the evidence shifts the belief toward the observed frequency. Plotting both curves makes the Bayesian update visually intuitive.

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

alpha_prior, beta_prior = 2, 2
k, n = 7, 10          # 7 heads in 10 flips
alpha_post = alpha_prior + k
beta_post  = beta_prior + (n - k)

x = np.linspace(0, 1, 300)
plt.plot(x, beta.pdf(x, alpha_prior, beta_prior),
         label=f'Prior Beta({alpha_prior},{beta_prior})')
plt.plot(x, beta.pdf(x, alpha_post, beta_post),
         label=f'Posterior Beta({alpha_post},{beta_post})')
plt.xlabel('θ (coin bias)')
plt.ylabel('Density')
plt.legend()
plt.title('Bayesian Update: Coin Flip')
plt.show()
print(f"Posterior mean: {alpha_post / (alpha_post + beta_post):.3f}")

The Naive Bayes classifier computes log P(class) + Σ log P(xᵢ|class) for each class and picks the highest log-posterior. For binary (0/1) features, each feature contributes log P(xᵢ=1|class) when xᵢ=1 and log P(xᵢ=0|class) when xᵢ=0 — both must be included so that absent features are not silently ignored. Using logarithms avoids numerical underflow when multiplying many small probabilities. Laplace smoothing (adding 1 to counts, 2 to the denominator for binary features) prevents zero probabilities for unseen feature–class combinations.

import numpy as np

class NaiveBayes:
    def fit(self, X, y):
        self.classes = np.unique(y)
        self.log_prior = {}
        self.log_likelihood = {}  # log_likelihood[c] shape: (2, n_features)
        for c in self.classes:
            Xc = X[y == c]
            self.log_prior[c] = np.log(len(Xc) / len(X))
            # Laplace smoothing: add 1 to counts, +2 to denominator (binary features)
            p1 = (Xc.sum(axis=0) + 1) / (len(Xc) + 2)
            # Row 0: log P(x=0|c), Row 1: log P(x=1|c)
            self.log_likelihood[c] = np.log(np.vstack([1 - p1, p1]))

    def predict(self, X):
        results = []
        for x in X:
            scores = {
                c: self.log_prior[c] + np.sum(
                    self.log_likelihood[c][0] * (1 - x) +
                    self.log_likelihood[c][1] * x
                )
                for c in self.classes
            }
            results.append(max(scores, key=scores.get))
        return np.array(results)

The Beta distribution is the conjugate prior for a binomial (two-category) probability; the Dirichlet distribution Dir(α₁,…,αₖ) generalises this to k categories, representing a distribution over probability vectors (p₁,…,pₖ) that sum to 1. The mean of component i is αᵢ/Σαⱼ, and larger αᵢ values express stronger prior belief in category i. The Dirichlet is the conjugate prior for the Multinomial distribution and is central to topic modelling algorithms such as Latent Dirichlet Allocation (LDA).

The marginal likelihood P(data) = ∫ P(data|θ)P(θ) dθ is the probability of the observed data averaged over all parameter values with respect to the prior; it serves as the normalising constant that makes the posterior integrate to 1. For conjugate models this integral has a closed-form solution, but for complex models with many parameters it requires numerical methods such as MCMC or variational inference. The marginal likelihood is also fundamental to Bayesian model comparison via Bayes factors, since it measures how well a model predicts data under its entire prior.

Prior sensitivity analysis tests how much the posterior and downstream conclusions change when different reasonable priors are used. If the posterior is robust across a range of priors, confidence in the conclusions is higher. If conclusions reverse with a different prior, it signals that the data alone are insufficient to determine the answer, and either more data or a better-justified prior is needed. In practice one runs the analysis with several priors (flat, weakly informative, strongly informative) and compares resulting posteriors and credible intervals.

The log-odds form is: log[P(H|E)/P(¬H|E)] = log[P(E|H)/P(E|¬H)] + log[P(H)/P(¬H)], i.e., log posterior odds = log likelihood ratio + log prior odds. Each new piece of evidence simply adds a log likelihood ratio term, turning multiplicative Bayesian updating into additive updates in log space. A positive log likelihood ratio supports H; a negative one supports ¬H. This additive form underlies logistic regression, where the model outputs a log-odds that is a linear function of features.

With limited data MLE can overfit drastically — if you flip a coin three times and get 3 heads, MLE gives θ̂=1.0, claiming the coin always lands heads. A Bayesian estimate with a reasonable prior such as Beta(3,3) yields a posterior Beta(3+3, 3+0) = Beta(6,3) with mean 6/9 ≈ 0.67 — far more sensible than 1.0. The prior acts as regularisation, preventing extreme estimates from small samples. As data grows both methods converge to the true parameter value.

Variational inference (VI) approximates the intractable posterior P(θ|data) by finding the distribution q(θ) from a simpler family (e.g., Gaussian) that is closest to the true posterior measured by KL divergence. Instead of sampling as in MCMC, VI turns inference into an optimisation problem: minimise KL[q(θ) || P(θ|data)], equivalently maximise the Evidence Lower Bound (ELBO). VI is typically much faster and more scalable than MCMC but produces an approximation rather than exact samples; MCMC is asymptotically exact but computationally expensive. Modern deep learning uses VI extensively in variational autoencoders (VAEs).

In fully Bayesian inference hyperparameters (parameters of the prior) are themselves given priors and inferred jointly with all other unknowns. Empirical Bayes instead estimates hyperparameters from the data by maximising the marginal likelihood P(data|hyperparams), then plugs these estimates into the prior before computing the posterior. It is a pragmatic compromise: more principled than pure MLE, and easier to compute than fully Bayesian, but it 'double-dips' the data — using it both to set the prior and to update it — which can underestimate posterior uncertainty. Empirical Bayes is widely used in genomics and hierarchical models.

The Metropolis algorithm proposes a new state from a symmetric proposal (e.g., Gaussian step), computes the acceptance ratio as the ratio of posterior densities, and accepts or rejects. Because the normalizing constant cancels, only the unnormalised posterior (prior × likelihood) is needed. The collected samples after a burn-in period approximate the posterior, and their mean estimates the posterior mean. The snippet below targets a Beta(3,7) posterior (proportional to θ²(1−θ)⁶) and verifies the sample mean ≈ 3/10 = 0.30.

import numpy as np

def log_posterior(theta):
    """Unnormalised: log of Beta(3,7) ∝ theta^2 * (1-theta)^6"""
    if theta <= 0 or theta >= 1:
        return -np.inf
    return 2 * np.log(theta) + 6 * np.log(1 - theta)

def metropolis(n_samples=10_000, step=0.1, burn_in=1_000):
    theta = 0.5  # starting value
    samples = []
    for _ in range(n_samples + burn_in):
        proposal = theta + np.random.normal(0, step)
        log_ratio = log_posterior(proposal) - log_posterior(theta)
        if np.log(np.random.rand()) < log_ratio:
            theta = proposal
        samples.append(theta)
    return np.array(samples[burn_in:])   # discard burn-in

samples = metropolis()
print(f"Sample mean: {samples.mean():.3f}  (true mean = 0.300)")

Explaining away (also called Berkson's paradox or inter-causal reasoning) occurs in a Bayesian network when two independent causes of the same effect become dependent once we observe the effect. For example, in the network 'Rain → Wet Grass ← Sprinkler', rain and the sprinkler are a priori independent. But if we know the grass is wet, observing that the sprinkler was ON makes rain less probable — the sprinkler 'explains away' the wetness, reducing our posterior belief in rain. This counter-intuitive dependency arises only through the shared effect node (a v-structure or collider).

Placing an independent Gaussian prior N(0, σ²) on each model weight w and then performing MAP estimation is equivalent to minimising the negative log-likelihood plus the term λ||w||² (L2 regularisation), where λ = 1/(2σ²). Smaller σ (tighter prior around zero) corresponds to stronger regularisation. This shows that regularisation is not just a computational trick but a Bayesian statement: we believe weights are probably small. Similarly, a Laplace (double-exponential) prior corresponds to L1 (Lasso) regularisation, which encourages sparsity.