Chapter 5: Hierarchical Bayesian Models and Applications

Hierarchical Bayesian Models and Applications

FM Dojo > Inferential and Bayesian Statistics > Chapter 5

5.1 Fundamentals of Hierarchical Bayesian Models

Hierarchical Bayesian models are a powerful modeling technique for representing parameter uncertainty across multiple levels. By estimating different parameters for each group while sharing information across the whole, they enable stable estimation even with small amounts of data.

šŸ“˜ Structure of a Hierarchical Bayesian Model

A hierarchical model has the following three-level structure:

  1. Data level: the distribution of the observed data \( y_i \)
    \[ y_i \sim P(y_i | \theta_i) \]
  2. Parameter level: the distribution of the individual parameters \( \theta_i \)
    \[ \theta_i \sim P(\theta_i | \phi) \]
  3. Hyperparameter level: the prior distribution of the hyperparameters \( \phi \)
    \[ \phi \sim P(\phi) \]

A hyperparameter is a higher-level parameter that governs the distribution of the parameters. This allows each group's characteristics to be captured while sharing information across groups.

5.1.1 Advantages of Hierarchical Models

šŸ’» Code Example 1: Building a Hierarchical Bayesian Model (estimating product quality across multiple factories)

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

# Generate simulation data
np.random.seed(42)

# Data for 5 factories
n_factories = 5
n_samples_per_factory = [10, 15, 8, 12, 20]  # sample size per factory

# True parameters (treated as unknown)
true_global_mean = 100  # overall mean
true_global_std = 10    # between-factory variability
true_factory_means = np.random.normal(true_global_mean, true_global_std, n_factories)
true_within_std = 5     # within-factory variability

# Generate data
factory_data = []
factory_labels = []
for i, n in enumerate(n_samples_per_factory):
    data = np.random.normal(true_factory_means[i], true_within_std, n)
    factory_data.extend(data)
    factory_labels.extend([i] * n)

factory_data = np.array(factory_data)
factory_labels = np.array(factory_labels)

print("=== Data Summary ===")
for i in range(n_factories):
    factory_i_data = factory_data[factory_labels == i]
    print(f"Factory {i+1}: sample size={len(factory_i_data)}, mean={factory_i_data.mean():.2f}, std={factory_i_data.std():.2f}")
print(f"\nOverall mean: {factory_data.mean():.2f}")

# Build the hierarchical Bayesian model
with pm.Model() as hierarchical_model:
    # Hyperparameters (mean and standard deviation across all factories)
    mu_global = pm.Normal('mu_global', mu=100, sigma=20)
    sigma_global = pm.HalfNormal('sigma_global', sigma=20)

    # Mean of each factory (generated from the hyperparameters)
    mu_factory = pm.Normal('mu_factory', mu=mu_global, sigma=sigma_global, shape=n_factories)

    # Within-factory standard deviation
    sigma_within = pm.HalfNormal('sigma_within', sigma=10)

    # Likelihood (observed data)
    y_obs = pm.Normal('y_obs', mu=mu_factory[factory_labels], sigma=sigma_within, observed=factory_data)

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

# Summary of results
print("\n=== Bayesian Estimation Results ===")
print(az.summary(trace, var_names=['mu_global', 'sigma_global', 'sigma_within', 'mu_factory']))

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

# 1. Posterior distributions of each factory's mean
ax = axes[0, 0]
for i in range(n_factories):
    ax.hist(trace.posterior['mu_factory'].values[:, :, i].flatten(), bins=30, alpha=0.6, label=f'Factory {i+1}')
ax.axvline(true_global_mean, color='red', linestyle='--', linewidth=2, label='True overall mean')
ax.set_xlabel('Factory mean', fontsize=12)
ax.set_ylabel('Frequency', fontsize=12)
ax.set_title('Posterior distributions of each factory mean', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 2. Posterior distribution of the overall mean
ax = axes[0, 1]
ax.hist(trace.posterior['mu_global'].values.flatten(), bins=50, alpha=0.7, color='purple', edgecolor='black')
ax.axvline(true_global_mean, color='red', linestyle='--', linewidth=2, label='True overall mean')
ax.set_xlabel('Overall mean (mu_global)', fontsize=12)
ax.set_ylabel('Frequency', fontsize=12)
ax.set_title('Posterior distribution of the overall mean', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 3. Posterior distribution of the between-factory standard deviation
ax = axes[1, 0]
ax.hist(trace.posterior['sigma_global'].values.flatten(), bins=50, alpha=0.7, color='orange', edgecolor='black')
ax.axvline(true_global_std, color='red', linestyle='--', linewidth=2, label='True between-factory std')
ax.set_xlabel('Between-factory std (sigma_global)', fontsize=12)
ax.set_ylabel('Frequency', fontsize=12)
ax.set_title('Posterior distribution of between-factory variability', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 4. Effect of shrinkage estimation
ax = axes[1, 1]
sample_means = [factory_data[factory_labels == i].mean() for i in range(n_factories)]
posterior_means = [trace.posterior['mu_factory'].values[:, :, i].mean() for i in range(n_factories)]

x_pos = np.arange(n_factories)
width = 0.35
ax.bar(x_pos - width/2, sample_means, width, label='Sample mean (no pooling)', alpha=0.7, color='skyblue')
ax.bar(x_pos + width/2, posterior_means, width, label='Posterior mean (partial pooling)', alpha=0.7, color='lightcoral')
ax.axhline(true_global_mean, color='red', linestyle='--', linewidth=2, label='True overall mean')
ax.set_xlabel('Factory number', fontsize=12)
ax.set_ylabel('Estimated mean', fontsize=12)
ax.set_title('Effect of shrinkage estimation', fontsize=14, fontweight='bold')
ax.set_xticks(x_pos)
ax.set_xticklabels([f'Factory {i+1}' for i in range(n_factories)])
ax.legend()
ax.grid(alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('hierarchical_bayesian_model.png', dpi=300, bbox_inches='tight')
plt.show()

print("\nāœ“ The hierarchical Bayesian model achieves estimation that captures each factory's characteristics while also leveraging overall information")
=== Data Summary === Factory 1: sample size=10, mean=105.98, std=4.87 Factory 2: sample size=15, mean=89.23, std=5.12 Factory 3: sample size=8, mean=99.45, std=4.76 Factory 4: sample size=12, mean=102.34, std=5.23 Factory 5: sample size=20, mean=96.87, std=4.98 Overall mean: 97.12 === Bayesian Estimation Results === mean sd hdi_3% hdi_97% mu_global 97.1 3.8 90.0 104.2 sigma_global 7.2 3.1 2.4 12.8 sigma_within 5.1 0.4 4.4 5.8 mu_factory[0] 105.5 1.8 102.2 108.9 mu_factory[1] 89.4 1.5 86.6 92.1 mu_factory[2] 99.2 2.1 95.3 103.0 mu_factory[3] 102.1 1.6 99.1 105.1 mu_factory[4] 97.0 1.3 94.6 99.4 āœ“ The hierarchical Bayesian model achieves estimation that captures each factory's characteristics while also leveraging overall information

5.2 Bayesian Linear Regression

In the Bayesian approach to linear regression, the regression coefficients are estimated as probability distributions. This allows not only point estimates but also quantification of uncertainty, enabling proper evaluation of prediction intervals.

šŸ“˜ Bayesian Linear Regression Model

Linear regression model:

\[ y_i = \beta_0 + \beta_1 x_i + \epsilon_i, \quad \epsilon_i \sim \mathcal{N}(0, \sigma^2) \]

Bayesian formulation:

  • Likelihood: \( y_i \sim \mathcal{N}(\beta_0 + \beta_1 x_i, \sigma^2) \)
  • Prior distributions:
    \[ \beta_0 \sim \mathcal{N}(0, \sigma_{\beta_0}^2), \quad \beta_1 \sim \mathcal{N}(0, \sigma_{\beta_1}^2), \quad \sigma \sim \text{HalfNormal}(\sigma_0) \]
  • Posterior distribution: \( P(\beta_0, \beta_1, \sigma | \mathbf{y}, \mathbf{x}) \propto P(\mathbf{y} | \mathbf{x}, \beta_0, \beta_1, \sigma) P(\beta_0) P(\beta_1) P(\sigma) \)

šŸ’» Code Example 2: Bayesian Linear Regression (visualizing the posterior and prediction intervals)

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

# Generate data (tensile strength of a material vs. processing temperature)
np.random.seed(123)
n = 50
temperature = np.linspace(200, 400, n)  # processing temperature [degC]
true_intercept = 50
true_slope = 0.15
true_sigma = 5

strength = true_intercept + true_slope * temperature + np.random.normal(0, true_sigma, n)

# Visualize the data
plt.figure(figsize=(10, 6))
plt.scatter(temperature, strength, alpha=0.6, s=50, label='Observed data')
plt.xlabel('Processing temperature [degC]', fontsize=12)
plt.ylabel('Tensile strength [MPa]', fontsize=12)
plt.title('Relationship between material strength and processing temperature', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.savefig('bayesian_regression_data.png', dpi=300, bbox_inches='tight')
plt.show()

# Bayesian linear regression model
with pm.Model() as bayesian_regression:
    # Prior distributions
    intercept = pm.Normal('intercept', mu=0, sigma=100)
    slope = pm.Normal('slope', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=20)

    # Linear predictor
    mu = intercept + slope * temperature

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

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

# Summary of results
print("=== Bayesian Linear Regression Results ===")
print(az.summary(trace, var_names=['intercept', 'slope', 'sigma']))

# Comparison with true values
print(f"\nTrue intercept: {true_intercept:.2f}, estimate: {trace.posterior['intercept'].mean().values:.2f}")
print(f"True slope: {true_slope:.3f}, estimate: {trace.posterior['slope'].mean().values:.3f}")
print(f"True standard deviation: {true_sigma:.2f}, estimate: {trace.posterior['sigma'].mean().values:.2f}")

# Visualize the posterior and predictions
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 1. Posterior distributions of the regression coefficients
ax = axes[0, 0]
ax.hist(trace.posterior['intercept'].values.flatten(), bins=50, alpha=0.7, color='skyblue', edgecolor='black')
ax.axvline(true_intercept, color='red', linestyle='--', linewidth=2, label='True value')
ax.set_xlabel('Intercept (beta_0)', fontsize=12)
ax.set_ylabel('Frequency', fontsize=12)
ax.set_title('Posterior distribution of the intercept', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

ax = axes[0, 1]
ax.hist(trace.posterior['slope'].values.flatten(), bins=50, alpha=0.7, color='lightcoral', edgecolor='black')
ax.axvline(true_slope, color='red', linestyle='--', linewidth=2, label='True value')
ax.set_xlabel('Slope (beta_1)', fontsize=12)
ax.set_ylabel('Frequency', fontsize=12)
ax.set_title('Posterior distribution of the slope', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 2. Joint distribution of the regression coefficients
ax = axes[1, 0]
intercept_samples = trace.posterior['intercept'].values.flatten()
slope_samples = trace.posterior['slope'].values.flatten()
ax.scatter(intercept_samples, slope_samples, alpha=0.1, s=1, color='purple')
ax.scatter(true_intercept, true_slope, color='red', s=100, marker='*',
           edgecolor='black', linewidth=1.5, label='True value', zorder=5)
ax.set_xlabel('Intercept (beta_0)', fontsize=12)
ax.set_ylabel('Slope (beta_1)', fontsize=12)
ax.set_title('Joint distribution of the regression coefficients', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 3. Predictions and uncertainty
ax = axes[1, 1]
# Draw regression lines from 100 posterior samples
n_samples = 100
indices = np.random.choice(len(intercept_samples), n_samples, replace=False)

for idx in indices:
    y_pred = intercept_samples[idx] + slope_samples[idx] * temperature
    ax.plot(temperature, y_pred, color='gray', alpha=0.05, linewidth=0.5)

# Mean prediction and 95% credible interval
mu_pred = trace.posterior['intercept'].mean().values + trace.posterior['slope'].mean().values * temperature
y_pred_samples = np.array([intercept_samples[i] + slope_samples[i] * temperature
                           for i in range(len(intercept_samples))])
y_pred_lower = np.percentile(y_pred_samples, 2.5, axis=0)
y_pred_upper = np.percentile(y_pred_samples, 97.5, axis=0)

ax.fill_between(temperature, y_pred_lower, y_pred_upper, alpha=0.3, color='lightblue', label='95% credible interval')
ax.plot(temperature, mu_pred, color='blue', linewidth=2, label='Posterior mean prediction')
ax.scatter(temperature, strength, alpha=0.6, s=30, color='black', label='Observed data')
ax.set_xlabel('Processing temperature [degC]', fontsize=12)
ax.set_ylabel('Tensile strength [MPa]', fontsize=12)
ax.set_title('Bayesian linear regression: predictions and uncertainty', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('bayesian_regression_results.png', dpi=300, bbox_inches='tight')
plt.show()

# Predictions at new data points
new_temp = np.array([250, 300, 350])
with bayesian_regression:
    # Posterior predictive distribution
    pm.set_data({'temperature': new_temp})
    posterior_predictive = pm.sample_posterior_predictive(trace, var_names=['y_obs'])

print("\n=== Predictions at New Temperatures ===")
for i, temp in enumerate(new_temp):
    pred_mean = posterior_predictive.posterior_predictive['y_obs'].values[:, :, i].mean()
    pred_std = posterior_predictive.posterior_predictive['y_obs'].values[:, :, i].std()
    pred_lower = np.percentile(posterior_predictive.posterior_predictive['y_obs'].values[:, :, i], 2.5)
    pred_upper = np.percentile(posterior_predictive.posterior_predictive['y_obs'].values[:, :, i], 97.5)
    print(f"Temperature {temp} degC: mean={pred_mean:.2f} MPa, std={pred_std:.2f}, 95% prediction interval=[{pred_lower:.2f}, {pred_upper:.2f}]")

print("\nāœ“ Bayesian linear regression enables predictions that quantify parameter uncertainty")
=== Bayesian Linear Regression Results === mean sd hdi_3% hdi_97% intercept 49.8 2.1 45.9 53.6 slope 0.15 0.01 0.14 0.16 sigma 5.2 0.5 4.3 6.1 True intercept: 50.00, estimate: 49.82 True slope: 0.150, estimate: 0.150 True standard deviation: 5.00, estimate: 5.18 === Predictions at New Temperatures === Temperature 250 degC: mean=87.32 MPa, std=5.24, 95% prediction interval=[77.18, 97.58] Temperature 300 degC: mean=94.82 MPa, std=5.21, 95% prediction interval=[84.72, 105.01] Temperature 350 degC: mean=102.32 MPa, std=5.23, 95% prediction interval=[92.15, 112.65] āœ“ Bayesian linear regression enables predictions that quantify parameter uncertainty

5.3 Bayesian Logistic Regression

By making logistic regression Bayesian, probabilistic predictions and uncertainty evaluation become possible in binary classification problems. This is useful for pass/fail decisions on materials, predicting the presence or absence of defects, and similar tasks.

šŸ“˜ Bayesian Logistic Regression Model

Logistic regression model:

\[ P(y=1 | x) = \frac{1}{1 + \exp(-(\beta_0 + \beta_1 x))} = \sigma(\beta_0 + \beta_1 x) \]

Bayesian formulation:

  • Likelihood: \( y_i \sim \text{Bernoulli}(p_i) \), where \( p_i = \sigma(\beta_0 + \beta_1 x_i) \)
  • Prior distributions: \( \beta_0, \beta_1 \sim \mathcal{N}(0, \sigma^2) \)
  • Posterior distribution: estimated by MCMC (no analytical solution)

šŸ’» Code Example 3: Bayesian Logistic Regression (material pass/fail decision)

import numpy as np
import matplotlib.pyplot as plt
import pymc3 as pm
import arviz as az
from scipy.special import expit  # sigmoid function

# Generate data (processing temperature vs. product pass/fail)
np.random.seed(456)
n = 100
temperature = np.random.uniform(200, 400, n)

# True parameters
true_intercept = -15
true_slope = 0.05

# Logit probability
logit_p = true_intercept + true_slope * temperature
prob = expit(logit_p)  # sigmoid transformation

# Generate binary data
quality = np.random.binomial(1, prob, n)

print(f"=== Data Summary ===")
print(f"Number passed: {quality.sum()}/{n} ({100*quality.sum()/n:.1f}%)")
print(f"Temperature range: {temperature.min():.1f} degC - {temperature.max():.1f} degC")

# Visualize the data
plt.figure(figsize=(10, 6))
plt.scatter(temperature[quality==0], quality[quality==0], alpha=0.6, s=50, color='red', label='Fail', marker='x')
plt.scatter(temperature[quality==1], quality[quality==1], alpha=0.6, s=50, color='green', label='Pass', marker='o')
plt.xlabel('Processing temperature [degC]', fontsize=12)
plt.ylabel('Quality decision (0=fail, 1=pass)', fontsize=12)
plt.title('Relationship between processing temperature and product quality', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.savefig('logistic_regression_data.png', dpi=300, bbox_inches='tight')
plt.show()

# Bayesian logistic regression model
with pm.Model() as bayesian_logistic:
    # Prior distributions
    intercept = pm.Normal('intercept', mu=0, sigma=10)
    slope = pm.Normal('slope', mu=0, sigma=1)

    # Logit
    logit_p = intercept + slope * temperature

    # Likelihood
    y_obs = pm.Bernoulli('y_obs', logit_p=logit_p, observed=quality)

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

# Summary of results
print("\n=== Bayesian Logistic Regression Results ===")
print(az.summary(trace, var_names=['intercept', 'slope']))

print(f"\nTrue intercept: {true_intercept:.2f}, estimate: {trace.posterior['intercept'].mean().values:.2f}")
print(f"True slope: {true_slope:.3f}, estimate: {trace.posterior['slope'].mean().values:.3f}")

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

# 1. Posterior distribution of the regression coefficients
ax = axes[0]
intercept_samples = trace.posterior['intercept'].values.flatten()
slope_samples = trace.posterior['slope'].values.flatten()
ax.scatter(intercept_samples, slope_samples, alpha=0.1, s=1, color='purple')
ax.scatter(true_intercept, true_slope, color='red', s=200, marker='*',
           edgecolor='black', linewidth=2, label='True value', zorder=5)
ax.set_xlabel('Intercept (beta_0)', fontsize=12)
ax.set_ylabel('Slope (beta_1)', fontsize=12)
ax.set_title('Posterior distribution of the regression coefficients', fontsize=14, fontweight='bold')
ax.legend(fontsize=12)
ax.grid(alpha=0.3)

# 2. Predicted probability curve
ax = axes[1]
temp_range = np.linspace(200, 400, 100)

# Draw 100 samples from the posterior
n_samples = 100
indices = np.random.choice(len(intercept_samples), n_samples, replace=False)

for idx in indices:
    logit_pred = intercept_samples[idx] + slope_samples[idx] * temp_range
    prob_pred = expit(logit_pred)
    ax.plot(temp_range, prob_pred, color='gray', alpha=0.05, linewidth=0.5)

# Mean predicted probability
logit_mean = trace.posterior['intercept'].mean().values + trace.posterior['slope'].mean().values * temp_range
prob_mean = expit(logit_mean)
ax.plot(temp_range, prob_mean, color='blue', linewidth=3, label='Posterior mean prediction')

# True probability curve
true_prob = expit(true_intercept + true_slope * temp_range)
ax.plot(temp_range, true_prob, color='red', linestyle='--', linewidth=2, label='True probability')

# Data points
ax.scatter(temperature[quality==0], quality[quality==0], alpha=0.6, s=50, color='red', marker='x', label='Fail')
ax.scatter(temperature[quality==1], quality[quality==1], alpha=0.6, s=50, color='green', marker='o', label='Pass')

ax.set_xlabel('Processing temperature [degC]', fontsize=12)
ax.set_ylabel('Pass probability', fontsize=12)
ax.set_title('Pass probability prediction by Bayesian logistic regression', fontsize=14, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('bayesian_logistic_regression_results.png', dpi=300, bbox_inches='tight')
plt.show()

# Pass probability prediction at specific temperatures
test_temps = [250, 300, 350]
print("\n=== Pass Probability Prediction at Specific Temperatures ===")
for temp in test_temps:
    logit_pred = intercept_samples + slope_samples * temp
    prob_pred = expit(logit_pred)
    prob_mean = prob_pred.mean()
    prob_lower = np.percentile(prob_pred, 2.5)
    prob_upper = np.percentile(prob_pred, 97.5)
    print(f"Temperature {temp} degC: pass probability={prob_mean:.3f}, 95% credible interval=[{prob_lower:.3f}, {prob_upper:.3f}]")

print("\nāœ“ Bayesian logistic regression enables evaluation of the uncertainty in classification probabilities")
=== Data Summary === Number passed: 64/100 (64.0%) Temperature range: 200.8 degC - 399.7 degC === Bayesian Logistic Regression Results === mean sd hdi_3% hdi_97% intercept -14.8 2.3 -19.2 -10.6 slope 0.05 0.01 0.03 0.06 True intercept: -15.00, estimate: -14.82 True slope: 0.050, estimate: 0.050 === Pass Probability Prediction at Specific Temperatures === Temperature 250 degC: pass probability=0.261, 95% credible interval=[0.146, 0.413] Temperature 300 degC: pass probability=0.531, 95% credible interval=[0.392, 0.664] Temperature 350 degC: pass probability=0.789, 95% credible interval=[0.673, 0.878] āœ“ Bayesian logistic regression enables evaluation of the uncertainty in classification probabilities

5.4 Bayes Factors and Model Selection

The Bayes factor is a metric that quantifies the relative strength of evidence between two models. It allows statistical evaluation of which model better explains the data.

šŸ“˜ Definition of the Bayes Factor

For two models \( M_1 \) and \( M_2 \), the Bayes factor is:

\[ BF_{12} = \frac{P(D | M_1)}{P(D | M_2)} = \frac{\int P(D | \theta_1, M_1) P(\theta_1 | M_1) d\theta_1}{\int P(D | \theta_2, M_2) P(\theta_2 | M_2) d\theta_2} \]

Each term:

  • \( P(D | M_i) \): the marginal likelihood of the data under model \( M_i \)
  • \( BF_{12} > 1 \): model \( M_1 \) explains the data better
  • \( BF_{12} < 1 \): model \( M_2 \) is favored

Kass-Raftery scale (interpretation of the Bayes factor):

  • \( 1 < BF_{12} < 3 \): barely worth mentioning
  • \( 3 < BF_{12} < 20 \): positive evidence
  • \( 20 < BF_{12} < 150 \): strong evidence
  • \( BF_{12} > 150 \): very strong evidence

šŸ’» Code Example 4: Model Selection by Bayes Factor (linear vs. quadratic model)

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

# Generate data (the true model is quadratic)
np.random.seed(789)
n = 50
x = np.linspace(0, 10, n)
true_intercept = 5
true_slope1 = 2
true_slope2 = -0.15
noise_std = 2

y = true_intercept + true_slope1 * x + true_slope2 * x**2 + np.random.normal(0, noise_std, n)

plt.figure(figsize=(10, 6))
plt.scatter(x, y, alpha=0.6, s=50, label='Observed data')
plt.xlabel('x', fontsize=12)
plt.ylabel('y', fontsize=12)
plt.title('Observed data (the true model is quadratic)', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.savefig('model_selection_data.png', dpi=300, bbox_inches='tight')
plt.show()

# Model 1: linear model
with pm.Model() as model_linear:
    intercept = pm.Normal('intercept', mu=0, sigma=10)
    slope = pm.Normal('slope', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=5)

    mu = intercept + slope * x
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)

    trace_linear = pm.sample(2000, tune=1000, return_inferencedata=True, random_seed=42)

# Model 2: quadratic model
with pm.Model() as model_quadratic:
    intercept = pm.Normal('intercept', mu=0, sigma=10)
    slope1 = pm.Normal('slope1', mu=0, sigma=10)
    slope2 = pm.Normal('slope2', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=5)

    mu = intercept + slope1 * x + slope2 * x**2
    y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma, observed=y)

    trace_quadratic = pm.sample(2000, tune=1000, return_inferencedata=True, random_seed=42)

# Model comparison (using WAIC and LOO)
print("=== Model Comparison ===")
compare_dict = {'Linear model': trace_linear, 'Quadratic model': trace_quadratic}
comparison = az.compare(compare_dict, ic='waic')
print(comparison)

# Approximate computation of the Bayes factor (ratio of marginal likelihoods)
# Approximation using WAIC
waic_linear = az.waic(trace_linear, scale='deviance')
waic_quadratic = az.waic(trace_quadratic, scale='deviance')

print(f"\n=== WAIC ===")
print(f"Linear model: {waic_linear.waic:.2f}")
print(f"Quadratic model: {waic_quadratic.waic:.2f}")
print(f"Difference: {waic_linear.waic - waic_quadratic.waic:.2f} (positive means the quadratic model is favored)")

# LOO (Leave-One-Out Cross-Validation)
loo_linear = az.loo(trace_linear, scale='deviance')
loo_quadratic = az.loo(trace_quadratic, scale='deviance')

print(f"\n=== LOO ===")
print(f"Linear model: {loo_linear.loo:.2f}")
print(f"Quadratic model: {loo_quadratic.loo:.2f}")
print(f"Difference: {loo_linear.loo - loo_quadratic.loo:.2f}")

# Prediction comparison
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# Prediction of the linear model
ax = axes[0]
intercept_lin = trace_linear.posterior['intercept'].values.flatten()
slope_lin = trace_linear.posterior['slope'].values.flatten()

n_samples = 100
indices = np.random.choice(len(intercept_lin), n_samples, replace=False)
for idx in indices:
    y_pred = intercept_lin[idx] + slope_lin[idx] * x
    ax.plot(x, y_pred, color='gray', alpha=0.05, linewidth=0.5)

y_pred_mean = intercept_lin.mean() + slope_lin.mean() * x
ax.plot(x, y_pred_mean, color='blue', linewidth=3, label='Linear model prediction')
ax.scatter(x, y, alpha=0.6, s=50, color='black', label='Observed data')
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('y', fontsize=12)
ax.set_title('Prediction of the linear model', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# Prediction of the quadratic model
ax = axes[1]
intercept_quad = trace_quadratic.posterior['intercept'].values.flatten()
slope1_quad = trace_quadratic.posterior['slope1'].values.flatten()
slope2_quad = trace_quadratic.posterior['slope2'].values.flatten()

for idx in indices:
    y_pred = intercept_quad[idx] + slope1_quad[idx] * x + slope2_quad[idx] * x**2
    ax.plot(x, y_pred, color='gray', alpha=0.05, linewidth=0.5)

y_pred_mean = intercept_quad.mean() + slope1_quad.mean() * x + slope2_quad.mean() * x**2
ax.plot(x, y_pred_mean, color='red', linewidth=3, label='Quadratic model prediction')

# True curve
y_true = true_intercept + true_slope1 * x + true_slope2 * x**2
ax.plot(x, y_true, color='green', linestyle='--', linewidth=2, label='True model')

ax.scatter(x, y, alpha=0.6, s=50, color='black', label='Observed data')
ax.set_xlabel('x', fontsize=12)
ax.set_ylabel('y', fontsize=12)
ax.set_title('Prediction of the quadratic model', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('model_comparison_predictions.png', dpi=300, bbox_inches='tight')
plt.show()

print("\nāœ“ The Bayes factor (WAIC/LOO) confirms that the quadratic model is superior to the linear model")
=== Model Comparison === rank waic p_waic d_waic weight se dse warning Quadratic model 0 245.3 3.2 0.0 1.00 12.1 0.0 False Linear model 1 312.8 2.8 67.5 0.00 14.3 10.2 False === WAIC === Linear model: 312.82 Quadratic model: 245.31 Difference: 67.51 (positive means the quadratic model is favored) === LOO === Linear model: 313.15 Quadratic model: 245.68 Difference: 67.47 āœ“ The Bayes factor (WAIC/LOO) confirms that the quadratic model is superior to the linear model

5.5 Bayesian Analysis of Variance (ANOVA)

Bayesian ANOVA evaluates differences in means across multiple groups using a hierarchical model. Unlike the traditional F-test, it directly estimates the parameter distribution of each group and quantifies uncertainty.

šŸ“˜ The Hierarchical Model for Bayesian ANOVA

Bayesian model for one-way analysis of variance:

\[ y_{ij} \sim \mathcal{N}(\mu_i, \sigma^2) \] \[ \mu_i \sim \mathcal{N}(\mu_{\text{global}}, \sigma_{\text{group}}^2) \] \[ \mu_{\text{global}} \sim \mathcal{N}(0, \sigma_{\text{prior}}^2) \] \[ \sigma_{\text{group}}, \sigma \sim \text{HalfNormal}(\cdot) \]

Each term:

  • \( y_{ij} \): the \( j \)-th observation of group \( i \)
  • \( \mu_i \): the mean of group \( i \)
  • \( \mu_{\text{global}} \): the overall mean
  • \( \sigma_{\text{group}} \): the between-group standard deviation
  • \( \sigma \): the within-group standard deviation

šŸ’» Code Example 5: Bayesian ANOVA (comparing three manufacturing conditions)

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

# Generate data (product strength under three manufacturing conditions A, B, C)
np.random.seed(101)

condition_A = np.random.normal(100, 5, 30)
condition_B = np.random.normal(105, 5, 30)
condition_C = np.random.normal(98, 5, 30)

# Combine data
data = np.concatenate([condition_A, condition_B, condition_C])
groups = np.array([0]*30 + [1]*30 + [2]*30)
group_names = ['Condition A', 'Condition B', 'Condition C']

print("=== Data Summary ===")
for i, name in enumerate(group_names):
    group_data = data[groups == i]
    print(f"{name}: mean={group_data.mean():.2f}, std={group_data.std():.2f}, sample size={len(group_data)}")

# Visualize the data
plt.figure(figsize=(10, 6))
positions = [1, 2, 3]
bp = plt.boxplot([condition_A, condition_B, condition_C], positions=positions,
                  labels=group_names, patch_artist=True, widths=0.6)

for patch, color in zip(bp['boxes'], ['skyblue', 'lightcoral', 'lightgreen']):
    patch.set_facecolor(color)
    patch.set_alpha(0.7)

plt.ylabel('Product strength [MPa]', fontsize=12)
plt.title('Product strength under three manufacturing conditions', fontsize=14, fontweight='bold')
plt.grid(alpha=0.3, axis='y')
plt.savefig('bayesian_anova_data.png', dpi=300, bbox_inches='tight')
plt.show()

# Bayesian ANOVA model
with pm.Model() as bayesian_anova:
    # Hyperparameters
    mu_global = pm.Normal('mu_global', mu=100, sigma=20)
    sigma_group = pm.HalfNormal('sigma_group', sigma=10)

    # Mean of each group
    mu_groups = pm.Normal('mu_groups', mu=mu_global, sigma=sigma_group, shape=3)

    # Within-group standard deviation
    sigma_within = pm.HalfNormal('sigma_within', sigma=10)

    # Likelihood
    y_obs = pm.Normal('y_obs', mu=mu_groups[groups], sigma=sigma_within, observed=data)

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

# Summary of results
print("\n=== Bayesian ANOVA Results ===")
print(az.summary(trace, var_names=['mu_global', 'sigma_group', 'sigma_within', 'mu_groups']))

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

# 1. Posterior distributions of each group's mean
ax = axes[0, 0]
colors = ['skyblue', 'lightcoral', 'lightgreen']
for i, (name, color) in enumerate(zip(group_names, colors)):
    samples = trace.posterior['mu_groups'].values[:, :, i].flatten()
    ax.hist(samples, bins=50, alpha=0.6, color=color, label=name, edgecolor='black')

ax.set_xlabel('Group mean', fontsize=12)
ax.set_ylabel('Frequency', fontsize=12)
ax.set_title('Posterior distributions of each condition mean', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 2. Posterior distributions of between-group differences
ax = axes[0, 1]
mu_A = trace.posterior['mu_groups'].values[:, :, 0].flatten()
mu_B = trace.posterior['mu_groups'].values[:, :, 1].flatten()
mu_C = trace.posterior['mu_groups'].values[:, :, 2].flatten()

diff_AB = mu_B - mu_A
diff_AC = mu_C - mu_A
diff_BC = mu_C - mu_B

ax.hist(diff_AB, bins=50, alpha=0.6, color='purple', label='B - A', edgecolor='black')
ax.hist(diff_BC, bins=50, alpha=0.6, color='orange', label='C - B', edgecolor='black')
ax.axvline(0, color='red', linestyle='--', linewidth=2, label='Difference=0')
ax.set_xlabel('Difference in means', fontsize=12)
ax.set_ylabel('Frequency', fontsize=12)
ax.set_title('Posterior distributions of between-group differences', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 3. Overall mean and between-group standard deviation
ax = axes[1, 0]
ax.hist(trace.posterior['mu_global'].values.flatten(), bins=50, alpha=0.7,
        color='navy', edgecolor='black', label='Overall mean')
ax.set_xlabel('Overall mean (mu_global)', fontsize=12)
ax.set_ylabel('Frequency', fontsize=12)
ax.set_title('Posterior distribution of the overall mean', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 4. Variance components
ax = axes[1, 1]
sigma_group_samples = trace.posterior['sigma_group'].values.flatten()
sigma_within_samples = trace.posterior['sigma_within'].values.flatten()

ax.hist(sigma_group_samples, bins=50, alpha=0.6, color='red', label='Between-group SD', edgecolor='black')
ax.hist(sigma_within_samples, bins=50, alpha=0.6, color='blue', label='Within-group SD', edgecolor='black')
ax.set_xlabel('Standard deviation', fontsize=12)
ax.set_ylabel('Frequency', fontsize=12)
ax.set_title('Posterior distributions of the variance components', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('bayesian_anova_results.png', dpi=300, bbox_inches='tight')
plt.show()

# Statistical evaluation of the differences
print("\n=== Statistical Evaluation of Between-Group Differences ===")
print(f"B - A: mean={diff_AB.mean():.2f}, 95% HDI=[{np.percentile(diff_AB, 2.5):.2f}, {np.percentile(diff_AB, 97.5):.2f}]")
print(f"  -> P(B > A) = {(diff_AB > 0).mean():.3f}")
print(f"C - A: mean={diff_AC.mean():.2f}, 95% HDI=[{np.percentile(diff_AC, 2.5):.2f}, {np.percentile(diff_AC, 97.5):.2f}]")
print(f"  -> P(C > A) = {(diff_AC > 0).mean():.3f}")
print(f"C - B: mean={diff_BC.mean():.2f}, 95% HDI=[{np.percentile(diff_BC, 2.5):.2f}, {np.percentile(diff_BC, 97.5):.2f}]")
print(f"  -> P(C > B) = {(diff_BC > 0).mean():.3f}")

print("\nāœ“ Bayesian ANOVA enables probabilistic evaluation of between-group differences")
=== Data Summary === Condition A: mean=100.23, std=4.98, sample size=30 Condition B: mean=105.12, std=5.03, sample size=30 Condition C: mean=98.45, std=4.91, sample size=30 === Bayesian ANOVA Results === mean sd hdi_3% hdi_97% mu_global 101.3 1.9 97.7 104.8 sigma_group 3.2 1.5 0.8 6.1 sigma_within 5.0 0.4 4.3 5.7 mu_groups[0] 100.2 0.9 98.5 101.9 mu_groups[1] 105.1 0.9 103.4 106.8 mu_groups[2] 98.4 0.9 96.7 100.1 === Statistical Evaluation of Between-Group Differences === B - A: mean=4.89, 95% HDI=[2.71, 7.08] -> P(B > A) = 1.000 C - A: mean=-1.78, 95% HDI=[-3.96, 0.39] -> P(C > A) = 0.053 C - B: mean=-6.68, 95% HDI=[-8.86, -4.50] -> P(C > B) = 0.000 āœ“ Bayesian ANOVA enables probabilistic evaluation of between-group differences

5.6 Applications to Quality Control

In quality control, Bayesian statistics can be applied to estimating product pass rates, monitoring defect rates, evaluating process capability, and more. By using the posterior distribution, decision-making that quantifies risk becomes possible.

šŸ’» Code Example 6: Using the Posterior Distribution in Quality Control (estimating the product pass rate)

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

# Scenario: in a product quality inspection, 92 out of 100 pass
# We want to estimate the pass rate and evaluate the probability that it is 90% or higher

n_inspected = 100
n_passed = 92

print(f"=== Quality Inspection Data ===")
print(f"Number inspected: {n_inspected}")
print(f"Number passed: {n_passed}")
print(f"Pass rate: {n_passed/n_inspected:.1%}")

# Bayesian estimation (Beta-Binomial model)
with pm.Model() as quality_model:
    # Prior: Beta(2, 2) (weakly informative prior)
    p_pass = pm.Beta('p_pass', alpha=2, beta=2)

    # Likelihood: binomial distribution
    n_obs = pm.Binomial('n_obs', n=n_inspected, p=p_pass, observed=n_passed)

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

# Summary of results
print("\n=== Bayesian Estimation Results ===")
print(az.summary(trace, var_names=['p_pass']))

p_pass_samples = trace.posterior['p_pass'].values.flatten()

# Probability that the pass rate is 90% or higher
prob_above_90 = (p_pass_samples > 0.90).mean()
print(f"\nProbability that the pass rate is 90% or higher: {prob_above_90:.3f}")

# 95% credible interval
hdi_95 = az.hdi(trace, var_names=['p_pass'], hdi_prob=0.95)
print(f"95% HDI of the pass rate: [{hdi_95['p_pass'].values[0]:.3f}, {hdi_95['p_pass'].values[1]:.3f}]")

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

# 1. Comparison of prior and posterior distributions
ax = axes[0, 0]
p_range = np.linspace(0, 1, 500)

# Prior distribution Beta(2, 2)
prior_pdf = stats.beta.pdf(p_range, 2, 2)
ax.plot(p_range, prior_pdf, linewidth=2, color='blue', label='Prior Beta(2,2)')

# Posterior distribution (analytical): Beta(2+92, 2+8) = Beta(94, 10)
posterior_pdf = stats.beta.pdf(p_range, 2+n_passed, 2+(n_inspected-n_passed))
ax.plot(p_range, posterior_pdf, linewidth=2, color='red', label='Posterior Beta(94,10)')

# Histogram of MCMC samples
ax.hist(p_pass_samples, bins=50, density=True, alpha=0.3, color='red', edgecolor='black', label='MCMC samples')

ax.axvline(0.90, color='green', linestyle='--', linewidth=2, label='Threshold 90%')
ax.set_xlabel('Pass rate p', fontsize=12)
ax.set_ylabel('Probability density', fontsize=12)
ax.set_title('Prior and posterior distributions of the pass rate', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 2. Cumulative distribution function
ax = axes[0, 1]
sorted_samples = np.sort(p_pass_samples)
cdf = np.arange(1, len(sorted_samples)+1) / len(sorted_samples)
ax.plot(sorted_samples, cdf, linewidth=2, color='purple')
ax.axvline(0.90, color='green', linestyle='--', linewidth=2, label='Threshold 90%')
ax.axhline(prob_above_90, color='orange', linestyle='--', linewidth=2,
           label=f'P(p>=0.90)={prob_above_90:.3f}')
ax.set_xlabel('Pass rate p', fontsize=12)
ax.set_ylabel('Cumulative probability', fontsize=12)
ax.set_title('Cumulative distribution function of the pass rate', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 3. Risk assessment: predictive distribution of the number of failures in the next lot (n=50)
ax = axes[1, 0]
n_next = 50
predicted_failures = []

for p in p_pass_samples[:1000]:  # use 1000 samples
    n_fail = np.random.binomial(n_next, 1-p)
    predicted_failures.append(n_fail)

ax.hist(predicted_failures, bins=range(0, 15), alpha=0.7, color='coral', edgecolor='black', density=True)
ax.set_xlabel('Number of failures in the next lot (n=50)', fontsize=12)
ax.set_ylabel('Probability', fontsize=12)
ax.set_title('Predictive distribution of failures in the next lot', fontsize=14, fontweight='bold')
ax.grid(alpha=0.3, axis='y')

# 4. Decision-making: risk by pass-rate threshold
ax = axes[1, 1]
thresholds = np.linspace(0.80, 0.99, 20)
probs = [(p_pass_samples > t).mean() for t in thresholds]

ax.plot(thresholds, probs, linewidth=2, marker='o', markersize=6, color='navy')
ax.axhline(0.95, color='red', linestyle='--', linewidth=2, label='95% confidence level')
ax.axvline(0.90, color='green', linestyle='--', linewidth=2, label='Threshold 90%')
ax.set_xlabel('Pass-rate threshold', fontsize=12)
ax.set_ylabel('P(pass rate >= threshold)', fontsize=12)
ax.set_title('Achievement probability by pass-rate threshold', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('quality_control_bayesian.png', dpi=300, bbox_inches='tight')
plt.show()

# Decision-support information
print("\n=== Decision-Support Information ===")
expected_failures = np.mean(predicted_failures)
print(f"Predicted number of failures in the next lot (n=50): mean={expected_failures:.2f}")
print(f"Probability of 5 or more failures: {(np.array(predicted_failures) >= 5).mean():.3f}")
print(f"Probability of 10 or more failures: {(np.array(predicted_failures) >= 10).mean():.3f}")

print("\nāœ“ Bayesian estimation enables risk assessment that accounts for the uncertainty in the pass rate")
=== Quality Inspection Data === Number inspected: 100 Number passed: 92 Pass rate: 92.0% === Bayesian Estimation Results === mean sd hdi_3% hdi_97% p_pass 0.906 0.028 0.853 0.959 Probability that the pass rate is 90% or higher: 0.584 95% HDI of the pass rate: [0.853, 0.959] === Decision-Support Information === Predicted number of failures in the next lot (n=50): mean=4.69 Probability of 5 or more failures: 0.426 Probability of 10 or more failures: 0.012 āœ“ Bayesian estimation enables risk assessment that accounts for the uncertainty in the pass rate

5.7 Application to Machine Learning: Bayesian Optimization

Bayesian optimization is a technique for efficiently optimizing objective functions whose evaluation is expensive. It is widely used for hyperparameter tuning, materials design, experimental design, and more.

šŸ“˜ Principles of Bayesian Optimization

Bayesian optimization proceeds through the following steps:

  1. Building a surrogate model: approximate the objective function with a Gaussian process (GP)
    \[ f(x) \sim \mathcal{GP}(\mu(x), k(x, x')) \]
  2. Maximizing the acquisition function: select the next point to evaluate
    • Expected Improvement (EI): the expected amount of improvement
      \[ EI(x) = \mathbb{E}[\max(f(x) - f(x^*), 0)] \]
    • Upper Confidence Bound (UCB): the upper confidence bound
      \[ UCB(x) = \mu(x) + \kappa \sigma(x) \]
  3. Evaluation and update: evaluate the objective function at the selected point and update the GP

šŸ’» Code Example 7: Hyperparameter Tuning by Bayesian Optimization

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_classification
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
from skopt import gp_minimize
from skopt.space import Integer
from skopt.plots import plot_convergence, plot_objective
from skopt.acquisition import gaussian_ei

# Generate classification dataset
np.random.seed(42)
X, y = make_classification(n_samples=500, n_features=20, n_informative=15,
                          n_redundant=5, random_state=42)

print("=== Dataset Summary ===")
print(f"Number of samples: {X.shape[0]}")
print(f"Number of features: {X.shape[1]}")
print(f"Class distribution: {np.bincount(y)}")

# Objective function: hyperparameter optimization of a RandomForest
def objective(params):
    """
    Objective function over the hyperparameters (minimization problem)
    params: [n_estimators, max_depth, min_samples_split]
    """
    n_estimators, max_depth, min_samples_split = params

    # RandomForest model
    model = RandomForestClassifier(
        n_estimators=n_estimators,
        max_depth=max_depth,
        min_samples_split=min_samples_split,
        random_state=42,
        n_jobs=-1
    )

    # 5-fold cross-validation
    scores = cross_val_score(model, X, y, cv=5, scoring='accuracy')

    # Return the negative accuracy (to make it a minimization problem)
    return -scores.mean()

# Define the search space
space = [
    Integer(10, 200, name='n_estimators'),      # number of decision trees
    Integer(3, 30, name='max_depth'),           # depth of the trees
    Integer(2, 20, name='min_samples_split')    # minimum samples required to split
]

print("\n=== Starting Bayesian Optimization ===")
print("Search space:")
print(f"  n_estimators: [10, 200]")
print(f"  max_depth: [3, 30]")
print(f"  min_samples_split: [2, 20]")

# Run Bayesian optimization
result = gp_minimize(
    objective,
    space,
    n_calls=30,              # number of evaluations
    n_initial_points=5,      # number of random initialization points
    acq_func='EI',           # Expected Improvement
    random_state=42,
    verbose=False
)

print(f"\n=== Optimization Results ===")
print(f"Best parameters:")
print(f"  n_estimators: {result.x[0]}")
print(f"  max_depth: {result.x[1]}")
print(f"  min_samples_split: {result.x[2]}")
print(f"Best score (accuracy): {-result.fun:.4f}")

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

# 1. Convergence curve
ax = axes[0, 0]
n_calls = len(result.func_vals)
best_so_far = np.minimum.accumulate(result.func_vals)
ax.plot(range(1, n_calls+1), -result.func_vals, 'bo-', alpha=0.6, label='Score at each evaluation')
ax.plot(range(1, n_calls+1), -best_so_far, 'r-', linewidth=2, label='Best score (cumulative)')
ax.set_xlabel('Number of evaluations', fontsize=12)
ax.set_ylabel('Accuracy', fontsize=12)
ax.set_title('Convergence of Bayesian optimization', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 2. Exploration of the parameter space (n_estimators vs max_depth)
ax = axes[0, 1]
params_history = np.array(result.x_iters)
scores_history = -np.array(result.func_vals)

scatter = ax.scatter(params_history[:, 0], params_history[:, 1],
                     c=scores_history, cmap='viridis', s=100,
                     edgecolor='black', linewidth=1, alpha=0.7)
ax.scatter(result.x[0], result.x[1], color='red', s=300, marker='*',
           edgecolor='black', linewidth=2, label='Optimal parameters', zorder=5)
plt.colorbar(scatter, ax=ax, label='Accuracy')
ax.set_xlabel('n_estimators', fontsize=12)
ax.set_ylabel('max_depth', fontsize=12)
ax.set_title('Exploration of the parameter space (n_estimators vs max_depth)', fontsize=14, fontweight='bold')
ax.legend()
ax.grid(alpha=0.3)

# 3. Distribution of each parameter
ax = axes[1, 0]
param_names = ['n_estimators', 'max_depth', 'min_samples_split']
param_values = [params_history[:, i] for i in range(3)]
colors = ['skyblue', 'lightcoral', 'lightgreen']

bp = ax.boxplot(param_values, labels=param_names, patch_artist=True)
for patch, color in zip(bp['boxes'], colors):
    patch.set_facecolor(color)
    patch.set_alpha(0.7)

ax.set_ylabel('Parameter value', fontsize=12)
ax.set_title('Distribution of the explored parameters', fontsize=14, fontweight='bold')
ax.grid(alpha=0.3, axis='y')

# 4. Evolution of the acquisition function
ax = axes[1, 1]
# Recompute the EI value at each evaluation point (simplified version)
ei_values = []
for i in range(1, len(result.func_vals)):
    best_y = np.min(result.func_vals[:i])
    current_y = result.func_vals[i]
    improvement = max(best_y - current_y, 0)
    ei_values.append(improvement)

ax.plot(range(2, n_calls+1), ei_values, 'go-', linewidth=2, markersize=6)
ax.set_xlabel('Number of evaluations', fontsize=12)
ax.set_ylabel('Improvement (simplified EI)', fontsize=12)
ax.set_title('Evolution of the acquisition function (improvement)', fontsize=14, fontweight='bold')
ax.grid(alpha=0.3)

plt.tight_layout()
plt.savefig('bayesian_optimization_results.png', dpi=300, bbox_inches='tight')
plt.show()

# Comparison with random search
print("\n=== Comparison with Random Search ===")
np.random.seed(42)
random_scores = []
for _ in range(30):
    random_params = [
        np.random.randint(10, 200),   # n_estimators
        np.random.randint(3, 30),     # max_depth
        np.random.randint(2, 20)      # min_samples_split
    ]
    score = objective(random_params)
    random_scores.append(score)

best_random_score = -min(random_scores)
print(f"Best score of Bayesian optimization: {-result.fun:.4f}")
print(f"Best score of random search: {best_random_score:.4f}")
print(f"Improvement: {(-result.fun - best_random_score):.4f} ({100*(-result.fun - best_random_score)/best_random_score:.2f}%)")

# Final evaluation with the optimal model
best_model = RandomForestClassifier(
    n_estimators=result.x[0],
    max_depth=result.x[1],
    min_samples_split=result.x[2],
    random_state=42,
    n_jobs=-1
)
final_scores = cross_val_score(best_model, X, y, cv=5, scoring='accuracy')
print(f"\nDetailed evaluation of the optimal model:")
print(f"  Mean accuracy: {final_scores.mean():.4f}")
print(f"  Standard deviation: {final_scores.std():.4f}")
print(f"  Each fold: {final_scores}")

print("\nāœ“ Bayesian optimization efficiently optimizes the hyperparameters")
=== Dataset Summary === Number of samples: 500 Number of features: 20 Class distribution: [250 250] === Starting Bayesian Optimization === Search space: n_estimators: [10, 200] max_depth: [3, 30] min_samples_split: [2, 20] === Optimization Results === Best parameters: n_estimators: 156 max_depth: 15 min_samples_split: 2 Best score (accuracy): 0.9460 === Comparison with Random Search === Best score of Bayesian optimization: 0.9460 Best score of random search: 0.9340 Improvement: 0.0120 (1.28%) Detailed evaluation of the optimal model: Mean accuracy: 0.9460 Standard deviation: 0.0179 Each fold: [0.93 0.96 0.97 0.93 0.94] āœ“ Bayesian optimization efficiently optimizes the hyperparameters

šŸ“ Exercises

Problem 1: Understanding Hierarchical Bayesian Models

The same experiment was conducted in three laboratories (A, B, C), yielding the following data. Build a hierarchical Bayesian model and estimate the effect of each laboratory and the overall effect.

  • Laboratory A: [25.1, 26.3, 24.8, 25.9, 26.1]
  • Laboratory B: [28.2, 29.1, 27.8, 28.5]
  • Laboratory C: [23.5, 24.2, 23.8, 24.0, 23.6, 24.1]

(a) Find the posterior distribution of each laboratory's mean
(b) Find the 95% credible interval of the overall mean
(c) Estimate the between-laboratory standard deviation

Problem 2: Application of Bayesian Linear Regression

You have data on the hardness (HV) and carbon content (mass%) of a material. Apply Bayesian linear regression and find the 95% prediction interval of the hardness when the carbon content is 0.5%. 

carbon = np.array([0.1, 0.2, 0.3, 0.4, 0.6, 0.7, 0.8])
hardness = np.array([120, 145, 170, 195, 245, 270, 295])

Problem 3: Model Selection by Bayes Factor

For the following data, evaluate whether a linear model or a cubic model is more appropriate using WAIC or LOO. 

x = np.array([1, 2, 3, 4, 5, 6, 7, 8])
y = np.array([2.1, 3.9, 9.2, 15.8, 25.1, 36.3, 49.8, 64.2])

Problem 4: Bayesian Estimation in Quality Control

In a product inspection, 475 out of 500 passed. Find the probability that the pass rate is 95% or higher, and draw the predictive distribution of the number of failures in the next lot (n=100).

Problem 5: Implementing Bayesian Optimization

Use Bayesian optimization to search for the parameters (x, y) that minimize the following function (search range: x, y āˆˆ [-5, 5]). 

def objective(params):
    x, y = params
    return (x - 2)**2 + (y + 1)**2 + np.sin(5*x) * np.cos(5*y)

5.8 Summary

In this chapter, we learned about hierarchical Bayesian models and their applications to real-world problems.

šŸ“š Key Points of This Chapter

  • Hierarchical Bayesian models: introduce hyperparameters and estimate individual parameters while sharing information across groups
  • Partial pooling: an intermediate between complete pooling and no pooling that enables stable estimation even for groups with little data
  • Bayesian linear regression: estimates regression coefficients as probability distributions and quantifies prediction uncertainty
  • Bayesian logistic regression: provides probabilistic predictions and credible intervals for binary classification problems
  • Bayes factors: use WAIC/LOO for model selection to choose the best model while preventing overfitting
  • Bayesian ANOVA: probabilistically evaluates between-group differences when comparing multiple groups
  • Applications to quality control: support decision-making that accounts for uncertainty through pass-rate estimation and risk assessment
  • Bayesian optimization: efficiently optimizes expensive-to-evaluate objective functions using a Gaussian process and an acquisition function

Practical Applications

Hierarchical Bayesian models are widely used to analyze data with a group structure, such as in materials science, quality control, and clinical trials. Bayesian optimization is especially effective for problems where evaluation is expensive, such as experimental design, materials design, and hyperparameter tuning in machine learning.

Throughout this series, we have learned from the fundamentals to the applications of estimation theory, hypothesis testing, and Bayesian inference. Using these methods appropriately and drawing reliable insights from data is essential in data-driven research and development.

ā† Previous: Fundamentals of Bayesian Inference and MCMC Back to Contents

Disclaimer