Chapter 4: Fundamentals of Bayesian Inference and MCMC | Inferential and Bayesian Statistics

4.1 Basic Concepts of Bayesian Inference

Bayesian inference is a statistical framework that combines prior knowledge with data to quantify uncertainty. Unlike frequentist inferential statistics, it treats parameters as random variables.

📘 Bayes' Theorem

For a parameter $ \theta $ and data $ D $:

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

Meaning of each term:

Posterior distribution $ P(\theta | D) $: the probability distribution of the parameter after observing the data
Likelihood $ P(D | \theta) $: the probability of obtaining data $ D $ given the parameter $ \theta $
Prior distribution $ P(\theta) $: beliefs about the parameter before observing the data
Marginal likelihood $ P(D) $: the probability of the data (normalizing constant)

In practice:

$$ P(\theta | D) \propto P(D | \theta) P(\theta) $$

The posterior distribution is proportional to the product of the likelihood and the prior distribution.

💻 Code Example 1: Implementing Bayes' Theorem (Coin-Toss Problem)

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

# Coin-toss problem: estimate the probability p that the coin lands heads
# Data: tossed 10 times, 7 heads
n_trials = 10
n_success = 7

# Prior distribution: Beta(2, 2) (weakly informative prior)
alpha_prior = 2
beta_prior = 2

# Likelihood: Binomial(n_success | n_trials, p)
# By the conjugate prior property, the posterior is also a Beta distribution
# Beta(alpha_prior + n_success, beta_prior + n_trials - n_success)
alpha_post = alpha_prior + n_success
beta_post = beta_prior + (n_trials - n_success)

print("=== Bayesian inference: coin-toss problem ===")
print(f"Data: {n_success} heads out of {n_trials} tosses")
print(f"Prior distribution: Beta({alpha_prior}, {beta_prior})")
print(f"Posterior distribution: Beta({alpha_post}, {beta_post})")
print(f"\nPosterior mean: {alpha_post/(alpha_post + beta_post):.4f}")
print(f"Posterior mode: {(alpha_post-1)/(alpha_post+beta_post-2):.4f}")

# 95% credible interval
ci_lower = stats.beta.ppf(0.025, alpha_post, beta_post)
ci_upper = stats.beta.ppf(0.975, alpha_post, beta_post)
print(f"95% credible interval: [{ci_lower:.4f}, {ci_upper:.4f}]")

# Comparison with the maximum likelihood estimate
mle = n_success / n_trials
print(f"\nFrequentist MLE: {mle:.4f}")

# Visualization
p_values = np.linspace(0, 1, 200)
prior_pdf = stats.beta.pdf(p_values, alpha_prior, beta_prior)
likelihood = stats.binom.pmf(n_success, n_trials, p_values)
posterior_pdf = stats.beta.pdf(p_values, alpha_post, beta_post)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Prior distribution
axes[0, 0].plot(p_values, prior_pdf, 'b-', linewidth=2)
axes[0, 0].fill_between(p_values, 0, prior_pdf, alpha=0.3, color='blue')
axes[0, 0].set_xlabel('p (probability of heads)')
axes[0, 0].set_ylabel('Probability density')
axes[0, 0].set_title(f'Prior distribution Beta({alpha_prior}, {beta_prior})')
axes[0, 0].grid(True, alpha=0.3)

# Likelihood function
axes[0, 1].plot(p_values, likelihood, 'g-', linewidth=2)
axes[0, 1].fill_between(p_values, 0, likelihood, alpha=0.3, color='green')
axes[0, 1].axvline(mle, color='red', linestyle='--', linewidth=2,
                   label=f'MLE: {mle:.2f}')
axes[0, 1].set_xlabel('p')
axes[0, 1].set_ylabel('Likelihood')
axes[0, 1].set_title(f'Likelihood function Bin({n_success}|{n_trials}, p)')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Posterior distribution
axes[1, 0].plot(p_values, posterior_pdf, 'r-', linewidth=2)
axes[1, 0].fill_between(p_values, 0, posterior_pdf, alpha=0.3, color='red')
axes[1, 0].axvline(alpha_post/(alpha_post+beta_post), color='blue',
                   linestyle='--', linewidth=2, label='Posterior mean')
axes[1, 0].axvspan(ci_lower, ci_upper, alpha=0.2, color='yellow',
                   label='95% credible interval')
axes[1, 0].set_xlabel('p')
axes[1, 0].set_ylabel('Probability density')
axes[1, 0].set_title(f'Posterior distribution Beta({alpha_post}, {beta_post})')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Overlay the three distributions
axes[1, 1].plot(p_values, prior_pdf / prior_pdf.max(), 'b-',
                linewidth=2, label='Prior (normalized)')
axes[1, 1].plot(p_values, likelihood / likelihood.max(), 'g--',
                linewidth=2, label='Likelihood (normalized)')
axes[1, 1].plot(p_values, posterior_pdf / posterior_pdf.max(), 'r-',
                linewidth=2, label='Posterior (normalized)')
axes[1, 1].axvline(mle, color='orange', linestyle=':', linewidth=2,
                   label=f'MLE: {mle:.2f}')
axes[1, 1].set_xlabel('p')
axes[1, 1].set_ylabel('Normalized probability density')
axes[1, 1].set_title('Visualization of Bayesian updating')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Demonstration of sequential Bayesian updating
print("\n=== Sequential Bayesian updating ===")
data_sequence = [1, 1, 0, 1, 1, 1, 0, 1, 1, 0]  # 1=heads, 0=tails
alpha_seq = alpha_prior
beta_seq = beta_prior

print(f"Initial prior distribution: Beta({alpha_seq}, {beta_seq})")
for i, outcome in enumerate(data_sequence, 1):
    if outcome == 1:
        alpha_seq += 1
    else:
        beta_seq += 1
    mean = alpha_seq / (alpha_seq + beta_seq)
    print(f"  After {i} data points: Beta({alpha_seq}, {beta_seq}), mean={mean:.4f}")

📌 Bayesian inference vs. frequentist inference

Bayesian: can be interpreted as "there is a 95% probability that p lies in the range 0.6-0.8"
Frequentist: "if such interval estimation is performed 100 times, 95 of the intervals will contain the true p"

Bayesian inference has the advantages that its probabilistic statements are intuitive and that prior knowledge can be incorporated formally.

4.2 Conjugate Prior Distributions

📘 Properties of Conjugate Priors

When the prior and posterior distributions belong to the same distribution family, that prior is called a conjugate prior.

Main conjugate pairs:

Binomial likelihood + Beta prior → Beta posterior
Poisson likelihood + Gamma prior → Gamma posterior
Normal likelihood (known variance) + Normal prior → Normal posterior
Normal likelihood (unknown variance) + Normal-Gamma prior → Normal-Gamma posterior

Using a conjugate prior lets us compute the posterior analytically, with the advantage that MCMC becomes unnecessary.

💻 Code Example 2: Conjugate Prior (Beta-Binomial)

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

# Compare Bayesian updating with various priors
n_trials = 20
n_success = 15

# Three types of priors
priors = [
    {"name": "Uninformative Uniform", "alpha": 1, "beta": 1},
    {"name": "Weakly informative Beta(2,2)", "alpha": 2, "beta": 2},
    {"name": "Strongly informative Beta(8,2)", "alpha": 8, "beta": 2}  # belief that heads is more likely
]

p_values = np.linspace(0, 1, 200)

fig, axes = plt.subplots(len(priors), 3, figsize=(15, 10))

for i, prior in enumerate(priors):
    alpha_pr = prior["alpha"]
    beta_pr = prior["beta"]
    alpha_po = alpha_pr + n_success
    beta_po = beta_pr + (n_trials - n_success)

    # Prior distribution
    prior_pdf = stats.beta.pdf(p_values, alpha_pr, beta_pr)
    axes[i, 0].plot(p_values, prior_pdf, 'b-', linewidth=2)
    axes[i, 0].fill_between(p_values, 0, prior_pdf, alpha=0.3, color='blue')
    axes[i, 0].set_title(f'{prior["name"]}\nBeta({alpha_pr}, {beta_pr})')
    axes[i, 0].set_ylabel('Probability density')
    axes[i, 0].grid(True, alpha=0.3)

    # Likelihood
    likelihood = stats.binom.pmf(n_success, n_trials, p_values)
    axes[i, 1].plot(p_values, likelihood, 'g-', linewidth=2)
    axes[i, 1].fill_between(p_values, 0, likelihood, alpha=0.3, color='green')
    axes[i, 1].set_title(f'Likelihood\nBin({n_success}|{n_trials}, p)')
    axes[i, 1].grid(True, alpha=0.3)

    # Posterior distribution
    posterior_pdf = stats.beta.pdf(p_values, alpha_po, beta_po)
    axes[i, 2].plot(p_values, posterior_pdf, 'r-', linewidth=2)
    axes[i, 2].fill_between(p_values, 0, posterior_pdf, alpha=0.3, color='red')
    post_mean = alpha_po / (alpha_po + beta_po)
    axes[i, 2].axvline(post_mean, color='blue', linestyle='--',
                       linewidth=2, label=f'Mean: {post_mean:.3f}')
    axes[i, 2].set_title(f'Posterior distribution\nBeta({alpha_po}, {beta_po})')
    axes[i, 2].legend()
    axes[i, 2].grid(True, alpha=0.3)

    if i == len(priors) - 1:
        axes[i, 0].set_xlabel('p')
        axes[i, 1].set_xlabel('p')
        axes[i, 2].set_xlabel('p')

plt.suptitle(f'How the posterior changes with different priors (Data: {n_success} successes out of {n_trials})',
             fontsize=14, y=1.00)
plt.tight_layout()
plt.show()

# Numerical summary
print("=== Influence of the prior distribution ===")
print(f"Data: {n_success} successes out of {n_trials} (MLE={n_success/n_trials:.3f})\n")
for prior in priors:
    alpha_pr = prior["alpha"]
    beta_pr = prior["beta"]
    alpha_po = alpha_pr + n_success
    beta_po = beta_pr + (n_trials - n_success)

    prior_mean = alpha_pr / (alpha_pr + beta_pr)
    post_mean = alpha_po / (alpha_po + beta_po)

    print(f"{prior['name']}:")
    print(f"  Prior mean: {prior_mean:.4f}")
    print(f"  Posterior mean: {post_mean:.4f}")
    print(f"  Change: {post_mean - prior_mean:+.4f}\n")

💻 Code Example 3: Conjugate Prior for the Normal Distribution

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

# Estimating the mean μ of normally distributed data (variance σ² is known)
sigma_known = 10  # known standard deviation
n = 15
np.random.seed(42)
true_mu = 100
data = np.random.normal(true_mu, sigma_known, n)

# Prior distribution: N(μ0, τ0²)
mu_0 = 95  # prior mean
tau_0 = 8  # prior standard deviation

# Sample statistics
data_mean = np.mean(data)
data_se = sigma_known / np.sqrt(n)

# Posterior distribution (conjugacy of the normal distribution)
# Working with precision (the inverse of variance) is convenient
precision_prior = 1 / tau_0**2
precision_likelihood = n / sigma_known**2
precision_post = precision_prior + precision_likelihood

mu_post = (precision_prior * mu_0 + precision_likelihood * data_mean) / precision_post
tau_post = 1 / np.sqrt(precision_post)

print("=== Conjugate Bayesian inference for the normal distribution ===")
print(f"Data: n={n}, sample mean={data_mean:.2f}, known σ={sigma_known}")
print(f"\nPrior distribution: N({mu_0}, {tau_0}²)")
print(f"Posterior distribution: N({mu_post:.2f}, {tau_post:.2f}²)")
print(f"\nFrequentist estimate:")
print(f"  Sample mean: {data_mean:.2f}")
print(f"  Standard error: {data_se:.2f}")

# Visualization
x = np.linspace(70, 120, 300)
prior_pdf = stats.norm.pdf(x, mu_0, tau_0)
likelihood_pdf = stats.norm.pdf(x, data_mean, data_se)
posterior_pdf = stats.norm.pdf(x, mu_post, tau_post)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# Overlay of the distributions
axes[0].plot(x, prior_pdf, 'b-', linewidth=2, label=f'Prior N({mu_0}, {tau_0}²)')
axes[0].plot(x, likelihood_pdf, 'g--', linewidth=2,
             label=f'Likelihood N({data_mean:.1f}, {data_se:.1f}²)')
axes[0].plot(x, posterior_pdf, 'r-', linewidth=2,
             label=f'Posterior N({mu_post:.1f}, {tau_post:.1f}²)')
axes[0].axvline(true_mu, color='black', linestyle=':', linewidth=2,
                label=f'True value: {true_mu}')
axes[0].set_xlabel('μ')
axes[0].set_ylabel('Probability density')
axes[0].set_title('Conjugate Bayesian inference for the normal distribution')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Effect of data size
sample_sizes = [5, 10, 20, 50, 100]
post_means = []
post_stds = []

np.random.seed(42)
for ns in sample_sizes:
    data_temp = np.random.normal(true_mu, sigma_known, ns)
    dm = np.mean(data_temp)
    prec_lik = ns / sigma_known**2
    prec_po = precision_prior + prec_lik
    mu_po = (precision_prior * mu_0 + prec_lik * dm) / prec_po
    tau_po = 1 / np.sqrt(prec_po)
    post_means.append(mu_po)
    post_stds.append(tau_po)

axes[1].errorbar(sample_sizes, post_means,
                 yerr=[2*s for s in post_stds],
                 fmt='o-', markersize=8, capsize=8, linewidth=2,
                 label='Posterior mean ±2SD')
axes[1].axhline(mu_0, color='blue', linestyle='--', linewidth=2,
                label=f'Prior mean: {mu_0}')
axes[1].axhline(true_mu, color='red', linestyle='--', linewidth=2,
                label=f'True value: {true_mu}')
axes[1].set_xlabel('Data size n')
axes[1].set_ylabel('Estimate of μ')
axes[1].set_title('Effect of data size (convergence from prior to true value)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
axes[1].set_xscale('log')

plt.tight_layout()
plt.show()

4.3 Markov Chain Monte Carlo (MCMC)

Complex posterior distributions often cannot be computed analytically, so sampling methods are required. MCMC is a technique that uses Markov chains to generate samples from the posterior distribution.

📘 Metropolis-Hastings Algorithm

A general-purpose algorithm for sampling from an arbitrary probability distribution $ p(\theta) $:

From the current state $ \theta^{(t)} $, generate a candidate $ \theta^* $ using the proposal distribution $ q(\theta^* | \theta^{(t)}) $
Compute the acceptance probability:
$$ \alpha = \min\left(1, \frac{p(\theta^*) q(\theta^{(t)} | \theta^*)}{p(\theta^{(t)}) q(\theta^* | \theta^{(t)})}\right) $$
Accept the candidate with probability $ \alpha $ ($ \theta^{(t+1)} = \theta^* $); otherwise reject it ($ \theta^{(t+1)} = \theta^{(t)} $)

For a symmetric proposal distribution ($ q(\theta^* | \theta) = q(\theta | \theta^*) $), the acceptance probability is:

$$ \alpha = \min\left(1, \frac{p(\theta^*)}{p(\theta^{(t)})}\right) $$

💻 Code Example 4: Implementing the Metropolis-Hastings Algorithm

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

# Target distribution: a mixture of two normal distributions (an example hard to handle analytically)
def target_distribution(theta):
    """Mixture of normals 0.3*N(0,1) + 0.7*N(5,1.5)"""
    component1 = 0.3 * stats.norm.pdf(theta, 0, 1)
    component2 = 0.7 * stats.norm.pdf(theta, 5, 1.5)
    return component1 + component2

# Metropolis-Hastings algorithm
def metropolis_hastings(target_func, n_samples, proposal_std=1.0, initial=0):
    samples = np.zeros(n_samples)
    samples[0] = initial
    n_accept = 0

    for t in range(1, n_samples):
        current = samples[t-1]

        # Proposal (symmetric normal distribution)
        proposal = current + np.random.normal(0, proposal_std)

        # Acceptance probability
        p_current = target_func(current)
        p_proposal = target_func(proposal)
        alpha = min(1, p_proposal / p_current) if p_current > 0 else 1

        # Accept/reject decision
        if np.random.rand() < alpha:
            samples[t] = proposal
            n_accept += 1
        else:
            samples[t] = current

    acceptance_rate = n_accept / (n_samples - 1)
    return samples, acceptance_rate

# Run MCMC
np.random.seed(42)
n_samples = 10000
samples, acc_rate = metropolis_hastings(target_distribution, n_samples,
                                        proposal_std=2.0, initial=0)

print("=== Metropolis-Hastings algorithm ===")
print(f"Number of samples: {n_samples}")
print(f"Acceptance rate: {acc_rate:.3f}")
print(f"Mean after burn-in: {np.mean(samples[1000:]):.4f}")
print(f"Standard deviation after burn-in: {np.std(samples[1000:]):.4f}")

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Trace plot
axes[0, 0].plot(samples[:500], 'b-', linewidth=0.5, alpha=0.7)
axes[0, 0].axhline(0, color='red', linestyle='--', linewidth=1, alpha=0.5)
axes[0, 0].axhline(5, color='green', linestyle='--', linewidth=1, alpha=0.5)
axes[0, 0].set_xlabel('Iteration')
axes[0, 0].set_ylabel('θ')
axes[0, 0].set_title(f'Trace plot (first 500 samples)\nAcceptance rate: {acc_rate:.2%}')
axes[0, 0].grid(True, alpha=0.3)

# Histogram and true distribution
burnin = 1000
x = np.linspace(-5, 10, 300)
true_pdf = target_distribution(x)

axes[0, 1].hist(samples[burnin:], bins=60, density=True, alpha=0.7,
                color='skyblue', edgecolor='black', label='MCMC samples')
axes[0, 1].plot(x, true_pdf, 'r-', linewidth=2, label='True distribution')
axes[0, 1].set_xlabel('θ')
axes[0, 1].set_ylabel('Probability density')
axes[0, 1].set_title(f'Estimation of the posterior (burn-in: {burnin} samples)')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Autocorrelation
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(samples[burnin:], lags=50, ax=axes[1, 0], alpha=0.05)
axes[1, 0].set_title('Autocorrelation function (checking convergence and independence)')
axes[1, 0].grid(True, alpha=0.3)

# Effect of the proposal distribution's standard deviation
proposal_stds = [0.1, 0.5, 1.0, 2.0, 5.0]
acceptance_rates = []

for pstd in proposal_stds:
    _, ar = metropolis_hastings(target_distribution, 5000,
                                proposal_std=pstd, initial=0)
    acceptance_rates.append(ar)

axes[1, 1].plot(proposal_stds, acceptance_rates, 'bo-',
                markersize=8, linewidth=2)
axes[1, 1].axhline(0.234, color='red', linestyle='--', linewidth=2,
                   label='Ideal acceptance rate (≈23.4%)')
axes[1, 1].set_xlabel('Standard deviation of the proposal distribution')
axes[1, 1].set_ylabel('Acceptance rate')
axes[1, 1].set_title('Tuning the proposal distribution')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].set_xscale('log')

plt.tight_layout()
plt.show()

📌 Practical points for MCMC

Burn-in: discard the first several thousand samples to remove the influence of the initial value
Acceptance rate: 20-40% is a good target (too low means slow convergence, too high means insufficient exploration)
Autocorrelation: check the independence between samples (thin them out if it is high)
Convergence diagnostics: multiple chains, the Gelman-Rubin statistic, trace plots

💻 Code Example 5: Bayesian Inference via Gibbs Sampling

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

# Sampling from a joint distribution of two variables
# Example: jointly estimating the mean μ and precision τ (inverse of variance) of a normal distribution

# Data generation
np.random.seed(42)
true_mu = 50
true_sigma = 10
n = 30
data = np.random.normal(true_mu, true_sigma, n)

# Prior distribution parameters
# μ ~ N(μ0, (λ0*τ)^{-1})
mu_0 = 45
lambda_0 = 0.1

# τ ~ Gamma(α0, β0)
alpha_0 = 2
beta_0 = 20

# Gibbs Sampling
def gibbs_sampling_normal(data, n_iter=5000, burnin=1000):
    n = len(data)
    data_mean = np.mean(data)
    data_sum_sq = np.sum((data - data_mean)**2)

    # Initial values
    mu = data_mean
    tau = 1 / np.var(data)

    # Store samples
    mu_samples = np.zeros(n_iter)
    tau_samples = np.zeros(n_iter)

    for i in range(n_iter):
        # Sample μ from its conditional posterior
        lambda_n = lambda_0 + n * tau
        mu_n = (lambda_0 * mu_0 + n * tau * data_mean) / lambda_n
        mu = np.random.normal(mu_n, 1/np.sqrt(lambda_n))

        # Sample τ from its conditional posterior
        alpha_n = alpha_0 + n/2
        beta_n = beta_0 + 0.5 * (np.sum((data - mu)**2) + lambda_0 * (mu - mu_0)**2)
        tau = np.random.gamma(alpha_n, 1/beta_n)

        mu_samples[i] = mu
        tau_samples[i] = tau

    return mu_samples[burnin:], tau_samples[burnin:]

# Run Gibbs Sampling
mu_samples, tau_samples = gibbs_sampling_normal(data, n_iter=10000, burnin=2000)
sigma_samples = 1 / np.sqrt(tau_samples)

print("=== Gibbs Sampling: parameter estimation for a normal distribution ===")
print(f"Data: n={n}, sample mean={np.mean(data):.2f}, sample SD={np.std(data):.2f}")
print(f"True values: μ={true_mu}, σ={true_sigma}")
print(f"\nPosterior estimates (mean):")
print(f"  μ: {np.mean(mu_samples):.2f} (95%CI: [{np.percentile(mu_samples, 2.5):.2f}, {np.percentile(mu_samples, 97.5):.2f}])")
print(f"  σ: {np.mean(sigma_samples):.2f} (95%CI: [{np.percentile(sigma_samples, 2.5):.2f}, {np.percentile(sigma_samples, 97.5):.2f}])")

# Visualization
fig, axes = plt.subplots(2, 3, figsize=(16, 10))

# Trace plot for μ
axes[0, 0].plot(mu_samples[:1000], 'b-', linewidth=0.5, alpha=0.7)
axes[0, 0].axhline(true_mu, color='red', linestyle='--', linewidth=2)
axes[0, 0].set_xlabel('Iteration')
axes[0, 0].set_ylabel('μ')
axes[0, 0].set_title('Trace plot for μ')
axes[0, 0].grid(True, alpha=0.3)

# Posterior distribution of μ
axes[0, 1].hist(mu_samples, bins=50, density=True, alpha=0.7,
                color='skyblue', edgecolor='black')
axes[0, 1].axvline(true_mu, color='red', linestyle='--',
                   linewidth=2, label=f'True value: {true_mu}')
axes[0, 1].axvline(np.mean(mu_samples), color='blue', linestyle='-',
                   linewidth=2, label=f'Posterior mean: {np.mean(mu_samples):.1f}')
axes[0, 1].set_xlabel('μ')
axes[0, 1].set_ylabel('Probability density')
axes[0, 1].set_title('Posterior distribution of μ')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

# Trace plot for σ
axes[0, 2].plot(sigma_samples[:1000], 'g-', linewidth=0.5, alpha=0.7)
axes[0, 2].axhline(true_sigma, color='red', linestyle='--', linewidth=2)
axes[0, 2].set_xlabel('Iteration')
axes[0, 2].set_ylabel('σ')
axes[0, 2].set_title('Trace plot for σ')
axes[0, 2].grid(True, alpha=0.3)

# Posterior distribution of σ
axes[1, 0].hist(sigma_samples, bins=50, density=True, alpha=0.7,
                color='lightgreen', edgecolor='black')
axes[1, 0].axvline(true_sigma, color='red', linestyle='--',
                   linewidth=2, label=f'True value: {true_sigma}')
axes[1, 0].axvline(np.mean(sigma_samples), color='green', linestyle='-',
                   linewidth=2, label=f'Posterior mean: {np.mean(sigma_samples):.1f}')
axes[1, 0].set_xlabel('σ')
axes[1, 0].set_ylabel('Probability density')
axes[1, 0].set_title('Posterior distribution of σ')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Joint posterior distribution
axes[1, 1].hexbin(mu_samples, sigma_samples, gridsize=50, cmap='Blues')
axes[1, 1].plot(true_mu, true_sigma, 'r*', markersize=15, label='True value')
axes[1, 1].set_xlabel('μ')
axes[1, 1].set_ylabel('σ')
axes[1, 1].set_title('Joint posterior distribution of μ and σ')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

# Posterior predictive distribution
n_pred = 1000
pred_samples = np.random.normal(mu_samples[:n_pred], sigma_samples[:n_pred])

axes[1, 2].hist(data, bins=15, density=True, alpha=0.5,
                color='orange', edgecolor='black', label='Observed data')
axes[1, 2].hist(pred_samples, bins=50, density=True, alpha=0.5,
                color='skyblue', edgecolor='black', label='Posterior predictive')
axes[1, 2].set_xlabel('Value')
axes[1, 2].set_ylabel('Density')
axes[1, 2].set_title('Posterior predictive distribution')
axes[1, 2].legend()
axes[1, 2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

4.4 Bayesian Inference with PyMC3

💻 Code Example 6: Bayesian Inference Using PyMC3

import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import arviz as az

# Data generation
np.random.seed(42)
true_alpha = 2.5
true_beta = 1.5
n = 50
x = np.linspace(0, 10, n)
y_true = true_alpha + true_beta * x
y = y_true + np.random.normal(0, 2, n)

print("=== Bayesian linear regression with PyMC3 ===")
print(f"Number of data points: {n}")
print(f"True parameters: α={true_alpha}, β={true_beta}")

# Define the PyMC3 model
with pm.Model() as model:
    # Prior distributions
    alpha = pm.Normal('alpha', mu=0, sd=10)
    beta = pm.Normal('beta', mu=0, sd=10)
    sigma = pm.HalfNormal('sigma', sd=5)

    # Likelihood
    mu = alpha + beta * x
    y_obs = pm.Normal('y_obs', mu=mu, sd=sigma, observed=y)

    # Sampling
    trace = pm.sample(2000, tune=1000, return_inferencedata=True,
                      random_seed=42, progressbar=False)

# Summary of results
print("\nSummary of the posterior distribution:")
print(az.summary(trace, var_names=['alpha', 'beta', 'sigma']))

# Visualization
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Data and regression line
axes[0, 0].scatter(x, y, alpha=0.5, color='blue', s=50, label='Data')
axes[0, 0].plot(x, y_true, 'r-', linewidth=2, label='True line')

# Regression lines sampled from the posterior
alpha_samples = trace.posterior['alpha'].values.flatten()
beta_samples = trace.posterior['beta'].values.flatten()

for i in np.random.choice(len(alpha_samples), 100):
    y_pred = alpha_samples[i] + beta_samples[i] * x
    axes[0, 0].plot(x, y_pred, 'gray', alpha=0.05)

# Regression line of the posterior mean
y_mean = np.mean(alpha_samples) + np.mean(beta_samples) * x
axes[0, 0].plot(x, y_mean, 'g--', linewidth=2,
                label='Posterior mean line')
axes[0, 0].set_xlabel('x')
axes[0, 0].set_ylabel('y')
axes[0, 0].set_title('Bayesian linear regression')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Posterior distributions (α, β, σ)
az.plot_posterior(trace, var_names=['alpha', 'beta', 'sigma'],
                  ref_val=[true_alpha, true_beta, 2], ax=axes[0, 1])
axes[0, 1].set_title('Posterior distributions of the parameters')

# Trace plot
az.plot_trace(trace, var_names=['alpha', 'beta', 'sigma'],
              compact=False)
plt.suptitle('Trace plot (convergence diagnostics)', y=1.02)

plt.tight_layout()
plt.show()

# Posterior predictive check
with model:
    ppc = pm.sample_posterior_predictive(trace, random_seed=42)

fig, ax = plt.subplots(figsize=(10, 6))
az.plot_ppc(az.from_pymc3(posterior_predictive=ppc, model=model), ax=ax)
ax.set_title('Posterior predictive check (assessing model fit)')
plt.tight_layout()
plt.show()

print("\nGelman-Rubin statistic (convergence diagnostic, closer to 1.0 is better):")
print(az.rhat(trace))

4.5 Bayesian Estimation of Material Properties

💻 Code Example 7: Bayesian Estimation and Prediction of Material Strength Data

import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import arviz as az
from scipy import stats

# Material strength data (small sample)
np.random.seed(42)
n_samples = 15
true_mean = 480
true_std = 25
strength_data = np.random.normal(true_mean, true_std, n_samples)

print("=== Bayesian estimation of material strength ===")
print(f"Number of data points: {n_samples} (small sample)")
print(f"Sample mean: {np.mean(strength_data):.2f} MPa")
print(f"Sample SD: {np.std(strength_data, ddof=1):.2f} MPa")

# PyMC3 model
with pm.Model() as strength_model:
    # Prior distributions (from past experience)
    mu = pm.Normal('mu', mu=500, sd=50)  # around 500±50 MPa from past data
    sigma = pm.HalfNormal('sigma', sd=30)  # variability up to about 30 MPa

    # Likelihood
    y_obs = pm.Normal('y_obs', mu=mu, sd=sigma, observed=strength_data)

    # Sampling
    trace = pm.sample(3000, tune=1500, return_inferencedata=True,
                      random_seed=42, progressbar=False)

# Results
mu_samples = trace.posterior['mu'].values.flatten()
sigma_samples = trace.posterior['sigma'].values.flatten()

mu_mean = np.mean(mu_samples)
mu_ci = np.percentile(mu_samples, [2.5, 97.5])
sigma_mean = np.mean(sigma_samples)
sigma_ci = np.percentile(sigma_samples, [2.5, 97.5])

print(f"\nBayesian estimation results:")
print(f"  Mean strength μ: {mu_mean:.2f} MPa (95%CI: [{mu_ci[0]:.2f}, {mu_ci[1]:.2f}])")
print(f"  Standard deviation σ: {sigma_mean:.2f} MPa (95%CI: [{sigma_ci[0]:.2f}, {sigma_ci[1]:.2f}])")

# Comparison with the frequentist estimate
freq_mean = np.mean(strength_data)
freq_std = np.std(strength_data, ddof=1)
freq_ci = stats.t.interval(0.95, n_samples-1,
                            loc=freq_mean,
                            scale=freq_std/np.sqrt(n_samples))

print(f"\nFrequentist estimate (for reference):")
print(f"  Mean strength: {freq_mean:.2f} MPa (95%CI: [{freq_ci[0]:.2f}, {freq_ci[1]:.2f}])")
print(f"  Standard deviation: {freq_std:.2f} MPa")

# Posterior predictive distribution (predicting the strength of the next sample)
n_pred = 5000
idx_random = np.random.choice(len(mu_samples), n_pred)
predicted_strength = np.random.normal(mu_samples[idx_random],
                                      sigma_samples[idx_random])

pred_mean = np.mean(predicted_strength)
pred_ci = np.percentile(predicted_strength, [2.5, 97.5])

print(f"\nStrength prediction for the next sample (posterior predictive distribution):")
print(f"  Predicted mean: {pred_mean:.2f} MPa")
print(f"  95% prediction interval: [{pred_ci[0]:.2f}, {pred_ci[1]:.2f}] MPa")

# Estimation of the failure probability
design_strength = 450  # design reference strength
prob_failure = np.mean(predicted_strength < design_strength)
print(f"\nProbability of being at or below the design reference strength of {design_strength} MPa: {prob_failure:.4f} ({prob_failure*100:.2f}%)")

# Visualization
fig = plt.figure(figsize=(16, 10))
gs = fig.add_gridspec(3, 2, hspace=0.3, wspace=0.3)

# Posterior distributions (μ, σ)
ax1 = fig.add_subplot(gs[0, 0])
ax1.hist(mu_samples, bins=50, density=True, alpha=0.7,
         color='skyblue', edgecolor='black')
ax1.axvline(true_mean, color='red', linestyle='--',
            linewidth=2, label=f'True value: {true_mean}')
ax1.axvline(mu_mean, color='blue', linestyle='-',
            linewidth=2, label=f'Posterior mean: {mu_mean:.1f}')
ax1.axvspan(mu_ci[0], mu_ci[1], alpha=0.2, color='blue',
            label='95% credible interval')
ax1.set_xlabel('Mean strength μ [MPa]')
ax1.set_ylabel('Probability density')
ax1.set_title('Posterior distribution of the mean strength')
ax1.legend()
ax1.grid(True, alpha=0.3)

ax2 = fig.add_subplot(gs[0, 1])
ax2.hist(sigma_samples, bins=50, density=True, alpha=0.7,
         color='lightgreen', edgecolor='black')
ax2.axvline(true_std, color='red', linestyle='--',
            linewidth=2, label=f'True value: {true_std}')
ax2.axvline(sigma_mean, color='green', linestyle='-',
            linewidth=2, label=f'Posterior mean: {sigma_mean:.1f}')
ax2.set_xlabel('Standard deviation σ [MPa]')
ax2.set_ylabel('Probability density')
ax2.set_title('Posterior distribution of the standard deviation')
ax2.legend()
ax2.grid(True, alpha=0.3)

# Joint posterior distribution
ax3 = fig.add_subplot(gs[1, :])
ax3.hexbin(mu_samples, sigma_samples, gridsize=50, cmap='Blues')
ax3.plot(true_mean, true_std, 'r*', markersize=20, label='True value')
ax3.set_xlabel('μ [MPa]')
ax3.set_ylabel('σ [MPa]')
ax3.set_title('Joint posterior distribution of μ and σ')
ax3.legend()
ax3.grid(True, alpha=0.3)

# Posterior predictive distribution and observed data
ax4 = fig.add_subplot(gs[2, 0])
ax4.hist(predicted_strength, bins=60, density=True, alpha=0.6,
         color='lightblue', edgecolor='black', label='Posterior predictive')
ax4.hist(strength_data, bins=10, density=True, alpha=0.6,
         color='orange', edgecolor='black', label='Observed data')
ax4.axvline(design_strength, color='red', linestyle='--',
            linewidth=2, label=f'Design reference: {design_strength} MPa')
ax4.axvspan(pred_ci[0], pred_ci[1], alpha=0.2, color='green',
            label='95% prediction interval')
ax4.set_xlabel('Strength [MPa]')
ax4.set_ylabel('Probability density')
ax4.set_title('Posterior predictive distribution')
ax4.legend()
ax4.grid(True, alpha=0.3)

# Visualization of the failure probability
ax5 = fig.add_subplot(gs[2, 1])
threshold_range = np.linspace(400, 550, 100)
failure_probs = []

for threshold in threshold_range:
    prob = np.mean(predicted_strength < threshold)
    failure_probs.append(prob)

ax5.plot(threshold_range, failure_probs, 'b-', linewidth=2)
ax5.axvline(design_strength, color='red', linestyle='--',
            linewidth=2, label=f'Design reference: {design_strength} MPa')
ax5.axhline(prob_failure, color='orange', linestyle=':',
            linewidth=2, label=f'Failure probability: {prob_failure:.3f}')
ax5.fill_between(threshold_range, 0, failure_probs,
                  where=(np.array(threshold_range) <= design_strength),
                  alpha=0.3, color='red')
ax5.set_xlabel('Strength threshold [MPa]')
ax5.set_ylabel('Failure probability')
ax5.set_title('Relationship between strength threshold and failure probability')
ax5.legend()
ax5.grid(True, alpha=0.3)

plt.suptitle(f'Bayesian estimation and prediction of material strength (n={n_samples})', fontsize=14)
plt.tight_layout()
plt.show()

# Practical summary
print("\n=== Practical interpretation ===")
print(f"✓ Even with a small sample (n={n_samples}), prior knowledge yields stable estimates")
print(f"✓ The mean strength lies in the range {mu_ci[0]:.1f}-{mu_ci[1]:.1f} MPa with 95% probability")
print(f"✓ The next sample lies in the range {pred_ci[0]:.1f}-{pred_ci[1]:.1f} MPa with 95% probability")
print(f"✓ The probability of failure at or below the design reference of {design_strength} MPa is {prob_failure*100:.2f}%")

📝 Exercises

Explain the difference between Bayesian inference and frequentist inference from the perspective of the interpretation of probability.
Discuss what the posterior distribution becomes when the prior is uninformative (a uniform distribution).
Explain the problems that arise in the Metropolis-Hastings algorithm when the acceptance rate is extremely high (>90%) and when it is extremely low (<10%).
Write code that uses PyMC3 to perform Bayesian estimation of the λ parameter of a Poisson distribution.

Summary

Bayesian inference integrates prior knowledge with data to quantify uncertainty
The posterior distribution is proportional to the product of the likelihood and the prior distribution
Using a conjugate prior allows the posterior distribution to be computed analytically
MCMC is a powerful technique for sampling from complex posterior distributions
The Metropolis-Hastings algorithm and Gibbs Sampling are representative MCMC algorithms
PyMC3 lets you concisely specify and perform inference for complex Bayesian models
In materials science, reliable estimation is possible even with small samples by leveraging prior knowledge
The posterior predictive distribution is useful for predicting new data and assessing uncertainty