Bayesian Inference & Modeling

Bayesian statistics provides a coherent mathematical framework for updating the probability of a hypothesis as more evidence or information becomes available. Unlike the frequentist approach, which treats parameters as fixed and data as random, the Bayesian paradigm treats parameters themselves as random variables, allowing for a more natural integration of prior knowledge and uncertainty.

1. Frequentist vs. Bayesian Philosophies

The fundamental divide in statistical inference rests on the interpretation of probability:

• Frequentist Probability: Defined as the limit of an event’s relative frequency in a large number of trials. If we say a coin has $P(\text{Heads})=0.5$, we mean that as the number of flips $n\to \infty$, the proportion of heads converges to $0.5$. Parameters $\theta$ are unknown but fixed constants.
• Bayesian Probability: Defined as a measure of “degree of belief” or “plausibility.” It is a quantifyable state of knowledge. This allows us to assign probabilities to one-off events (e.g., “the probability that it rained on Mars yesterday”) where frequentist repetition is impossible. Parameters $\theta$ are random variables characterized by probability distributions.

2. The Bayesian Framework

The engine of Bayesian inference is Bayes’ Theorem. Let $y$ represent the observed data and $\theta$ denote the parameters of the model.

The Posterior Distribution

The goal is to compute the Posterior Distribution $P(\theta |y)$, which represents our updated belief about the parameters after observing the data:

$P(\theta |y)=\frac{P(y|\theta )P(\theta )}{P(y)}$

Where:

• Prior $P(\theta )$: The distribution representing information about $\theta$ before the data $y$ is collected.
• Likelihood $P(y|\theta )$: The probability of the data $y$ occurring given a specific parameter value $\theta$.
• Evidence (Marginal Likelihood) $P(y)$: The total probability of the data, integrated over all possible parameters: $P(y)={\int }_{\Theta }P(y|\theta )P(\theta )d\theta$. This serves as a normalizing constant.

In most applications, we focus on the kernel of the distribution: $P(\theta |y)\propto P(y|\theta )\times P(\theta )$ This captures the intuition that the posterior is a compromise between the evidence provided by the data (likelihood) and our pre-existing knowledge (prior).

3. Conjugate Priors

A significant portion of analytic Bayesian statistics relies on Conjugacy. A prior $P(\theta )$ is conjugate to a likelihood $P(y|\theta )$ if the resulting posterior $P(\theta |y)$ belongs to the same family of distributions as the prior.

Beta-Binomial Model

Consider $y$ successes in $n$ Bernoulli trials. The likelihood is Binomial: $P(y|\theta )=(\genfrac{}{}{0ex}{}{n}{y}){\theta }^y(1-\theta {)}^{n-y}$ If we choose a Beta distribution as our prior, $\theta \sim \text{Beta}(\alpha ,\beta )$: $P(\theta )=\frac{1}{B(\alpha ,\beta )}{\theta }^{\alpha -1}(1-\theta {)}^{\beta -1}$ The posterior will be: $P(\theta |y)\propto {\theta }^{y+\alpha -1}(1-\theta {)}^{n-y+\beta -1}$ Which is recognized as $\text{Beta}(\alpha +y,\beta +n-y)$. This provides a simple rule: $\alpha$ and $\beta$ can be interpreted as “prior successes” and “prior failures.”

Normal-Normal Model

If $y\sim \mathcal{N}(\theta ,{\sigma }^2)$ and the prior $\theta \sim \mathcal{N}({\mu }_0,{\tau }_0^2)$, the posterior is also Normal with: ${\mu }_{post}=\left(\frac{{\mu }_0}{{\tau }_0^2}+\frac{\sum y_i}{{\sigma }^2}\right)/\left(\frac{1}{{\tau }_0^2}+\frac{n}{{\sigma }^2}\right)$ This demonstrates that the posterior mean is a precision-weighted average of the prior mean and the sample mean.

4. Uninformative and Objective Priors

When we have no prior information, we want a prior that adds minimal information.

• Principle of Indifference: Using a flat (Uniform) prior. However, a uniform prior on $\sigma$ is not uniform on ${\sigma }^2$ (variance).
• Jeffreys’ Prior: An objective prior that is invariant under reparameterization. It is proportional to the square root of the determinant of the Fisher Information $I(\theta )$: $P(\theta )\propto \sqrt{|I(\theta )|}$ For a Binomial distribution, the Jeffreys prior is $\text{Beta}(0.5,0.5)$.

5. Point Estimation: Mean, Median, and MAP

Unlike frequentist statistics which gives a point estimate (MLE), Bayes gives a distribution. If a single number is required, we use decision theory:

1. Posterior Mean: $E[\theta |y]$. Minimizes the expected squared error loss.
2. Posterior Median: The value $m$ such that $P(\theta \le m|y)=0.5$. Minimizes the expected absolute error loss.
3. Maximum A Posteriori (MAP): The mode of the posterior. ${\hat{\theta }}_{MAP}=\arg {\max }_{\theta }[P(y|\theta )P(\theta )]$ MAP is the Bayesian analog to Maximum Likelihood Estimation.

6. Credible Intervals: HPD Intervals

A Bayesian $100(1-\alpha )\%$ Credible Interval is an interval in which the true parameter lies with probability $1-\alpha$.

• Equal-Tailed Interval: The interval between the $\alpha /2$ and $1-\alpha /2$ quantiles.
• Highest Posterior Density (HPD) Interval: The shortest possible interval such that the probability $1-\alpha$ is contained within it. In an HPD interval, every point inside has a higher posterior density than any point outside. For asymmetric distributions, HPD intervals are superior to equal-tailed ones.

7. Model Choice: Bayes Factors and BIC

Bayes Factors

To compare model $M_1$ and $M_2$, we use the ratio of the marginal likelihoods (the evidence): $K=\frac{P(y|M_1)}{P(y|M_2)}=\frac{\int P(y|{\theta }_1,M_1)P({\theta }_1|M_1)d{\theta }_1}{\int P(y|{\theta }_2,M_2)P({\theta }_2|M_2)d{\theta }_2}$ A Bayes Factor $K>10$ is generally considered “strong” evidence for $M_1$.

BIC

The Bayesian Information Criterion is an asymptotic approximation to the Bayes Factor: $\text{BIC}=k\ln (n)-2\ln (\hat{L})$ where $k$ is the number of parameters and $\hat{L}$ is the maximum likelihood. Lower BIC indicates a better model.

8. Computational Bayesian Statistics: MCMC

For complex models, the integral in the denominator of Bayes’ Theorem ($P(y)$) is intractable. We use Markov Chain Monte Carlo (MCMC) to sample from the posterior without knowing the normalizing constant.

The Metropolis-Hastings Algorithm

Metropolis-Hastings generates a sequence of samples such that the distribution of these samples converges to the posterior.

1. Start at ${\theta }^{(0)}$.
2. Propose ${\theta }^{\ast }$ from a proposal distribution $q({\theta }^{\ast }|{\theta }^{(t)})$.
3. Calculate the acceptance ratio: $a=\frac{P(y|{\theta }^{\ast })P({\theta }^{\ast })q({\theta }^{(t)}|{\theta }^{\ast })}{P(y|{\theta }^{(t)})P({\theta }^{(t)})q({\theta }^{\ast }|{\theta }^{(t)})}$
4. Accept ${\theta }^{\ast }$ with probability $\min (1,a)$.

Python: Metropolis-Hastings From Scratch

The following script estimates the mean $\mu$ of a Normal distribution with a known variance ${\sigma }^2=1$, using a wide Normal prior.

import numpy as np

def metropolis_hastings(data, p_mu, p_sd, n_iter, prop_sd):
    """
    Estimates the mean of a Normal distribution.
    data: Observed data points
    p_mu, p_sd: Prior mean and standard deviation
    n_iter: Number of MCMC iterations
    prop_sd: Standard deviation of the proposal distribution (tuning parameter)
    """
    theta_curr = 0.0 # Initial guess
    chain = []
    
    def log_post_unnorm(theta):
        # Log-Likelihood: sum of log-PDFs of Normal(theta, 1)
        log_lik = -0.5 * np.sum((data - theta)**2)
        # Log-Prior: log-PDF of Normal(p_mu, p_sd)
        log_pri = -0.5 * ((theta - p_mu) / p_sd)**2
        return log_lik + log_pri

    for _ in range(n_iter):
        # 1. Propose new theta
        theta_prop = np.random.normal(theta_curr, prop_sd)
        
        # 2. Acceptance ratio
        log_acc = log_post_unnorm(theta_prop) - log_post_unnorm(theta_curr)
        
        # 3. Decision
        if np.log(np.random.rand()) < log_acc:
            theta_curr = theta_prop
        
        chain.append(theta_curr)
        
    return np.array(chain)

# Setup
np.random.seed(42)
true_mean = 3.5
data = np.random.normal(true_mean, 1, size=100)

# Execution
trace = metropolis_hastings(data, p_mu=0, p_sd=10, n_iter=5000, prop_sd=0.2)

# Results
burn_in = 1000
final_estimate = np.mean(trace[burn_in:])
print(f'Empirical Mean: {np.mean(data):.3f}')
print(f'Bayesian Posterior Mean: {final_estimate:.3f}')