Chapter 2: Interval Estimation and Confidence Intervals | Inferential and Bayesian Statistics

2.1 The Concept of a Confidence Interval

Point estimation gives a single estimate of a parameter, but it carries no information about how reliable that estimate is. Interval estimation presents, in probabilistic terms, a range that is likely to contain the parameter.

📘 Definition of a Confidence Interval

A confidence interval at confidence level $ 1-\alpha $ for a parameter $ \theta $ is, using the statistics $ L(X_1, \ldots, X_n) $ and $ U(X_1, \ldots, X_n) $:

$$ P(L \leq \theta \leq U) = 1 - \alpha $$

the interval $ [L, U] $ satisfying this. Typically $ \alpha = 0.05 $ (a 95% confidence interval) is used.

📌 An Important Interpretation
A "95% confidence interval" does not mean "the probability that the parameter lies in this interval is 95%." Rather, it means "if we perform this kind of interval estimation 100 times, about 95 of the resulting intervals will contain the true parameter" (the frequentist interpretation).

2.2 Confidence Interval for the Mean of a Normal Population

2.2.1 Case of Known Population Variance (z-interval)

📘 Confidence Interval When the Population Variance Is Known

When the population is $ N(\mu, \sigma^2) $ and $ \sigma^2 $ is known:

$$ \bar{X} \pm z_{\alpha/2} \frac{\sigma}{\sqrt{n}} $$

where $ z_{\alpha/2} $ is the upper $ \alpha/2 $ point of the standard normal distribution (e.g., for a 95% confidence interval, $ z_{0.025} = 1.96 $).

💻 Code Example 1: Confidence Interval for the Mean of a Normal Population (Known Population Variance)

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

# True parameters
mu_true = 100
sigma_true = 15  # known

# Data generation
np.random.seed(42)
n = 25
data = np.random.normal(mu_true, sigma_true, n)

# Sample mean
x_bar = np.mean(data)

# 95% confidence interval (z-interval)
alpha = 0.05
z_critical = stats.norm.ppf(1 - alpha/2)
margin_of_error = z_critical * sigma_true / np.sqrt(n)
ci_lower = x_bar - margin_of_error
ci_upper = x_bar + margin_of_error

print("=== Confidence Interval for the Population Mean (Known Variance) ===")
print(f"Sample size: {n}")
print(f"Sample mean: {x_bar:.2f}")
print(f"Population standard deviation (known): {sigma_true}")
print(f"z_{{{alpha/2}}} = {z_critical:.4f}")
print(f"Margin of error: ±{margin_of_error:.2f}")
print(f"95% confidence interval: [{ci_lower:.2f}, {ci_upper:.2f}]")
print(f"True population mean {mu_true} is inside the interval: {ci_lower <= mu_true <= ci_upper}")

# Verify coverage probability via simulation
n_simulations = 1000
coverage_count = 0
ci_data = []

np.random.seed(123)
for _ in range(n_simulations):
    sample = np.random.normal(mu_true, sigma_true, n)
    sample_mean = np.mean(sample)
    ci_l = sample_mean - margin_of_error
    ci_u = sample_mean + margin_of_error
    ci_data.append((ci_l, ci_u))
    if ci_l <= mu_true <= ci_upper:
        coverage_count += 1

coverage_rate = coverage_count / n_simulations
print(f"\nSimulation ({n_simulations} runs):")
print(f"Coverage probability: {coverage_rate:.3f} (theoretical value: 0.95)")

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

# Distribution and confidence interval
axes[0].hist(data, bins=15, density=True, alpha=0.7,
             color='skyblue', edgecolor='black', label='Data')
x_range = np.linspace(mu_true - 4*sigma_true, mu_true + 4*sigma_true, 200)
axes[0].plot(x_range, stats.norm.pdf(x_range, mu_true, sigma_true),
             'r-', linewidth=2, label=f'True population N({mu_true}, {sigma_true}²)')
axes[0].axvline(x_bar, color='blue', linestyle='--', linewidth=2,
                label=f'Sample mean: {x_bar:.1f}')
axes[0].axvspan(ci_lower, ci_upper, alpha=0.3, color='green',
                label=f'95% confidence interval')
axes[0].axvline(mu_true, color='red', linestyle=':', linewidth=2,
                label=f'True population mean: {mu_true}')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Probability density')
axes[0].set_title('95% Confidence Interval for the Population Mean')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Simulation of confidence intervals (showing the first 50 runs)
for i in range(min(50, len(ci_data))):
    ci_l, ci_u = ci_data[i]
    color = 'green' if ci_l <= mu_true <= ci_u else 'red'
    axes[1].plot([ci_l, ci_u], [i, i], color=color, linewidth=1, alpha=0.6)
    axes[1].plot((ci_l + ci_u)/2, i, 'o', color=color, markersize=2)

axes[1].axvline(mu_true, color='blue', linestyle='--', linewidth=2,
                label=f'True population mean: {mu_true}')
axes[1].set_xlabel('Interval')
axes[1].set_ylabel('Simulation number')
axes[1].set_title(f'Simulation of Confidence Intervals (coverage: {coverage_rate:.1%})')
axes[1].legend()
axes[1].grid(True, alpha=0.3, axis='x')

plt.tight_layout()
plt.show()

2.2.2 Case of Unknown Population Variance (t-interval)

📘 Confidence Interval When the Population Variance Is Unknown (t-distribution)

When the population is $ N(\mu, \sigma^2) $ and $ \sigma^2 $ is unknown:

$$ \bar{X} \pm t_{\alpha/2, n-1} \frac{S}{\sqrt{n}} $$

where $ S = \sqrt{\frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X})^2} $ is the sample standard deviation, and

$ t_{\alpha/2, n-1} $ is the upper $ \alpha/2 $ point of the t-distribution with $ n-1 $ degrees of freedom.

💻 Code Example 2: Confidence Interval Using the t-distribution (Small Sample)

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

# True parameters
mu_true = 75
sigma_true = 12

# Small-sample case
np.random.seed(42)
n = 10  # small sample
data = np.random.normal(mu_true, sigma_true, n)

# Sample statistics
x_bar = np.mean(data)
s = np.std(data, ddof=1)  # unbiased standard deviation

# 95% confidence interval (t-interval)
alpha = 0.05
t_critical = stats.t.ppf(1 - alpha/2, df=n-1)
margin_of_error_t = t_critical * s / np.sqrt(n)
ci_lower_t = x_bar - margin_of_error_t
ci_upper_t = x_bar + margin_of_error_t

# For comparison, also compute the z-interval (the incorrect method)
z_critical = stats.norm.ppf(1 - alpha/2)
margin_of_error_z = z_critical * s / np.sqrt(n)
ci_lower_z = x_bar - margin_of_error_z
ci_upper_z = x_bar + margin_of_error_z

print("=== Confidence Interval for the Population Mean (Unknown Variance, Small Sample) ===")
print(f"Sample size: {n}")
print(f"Sample mean: {x_bar:.2f}")
print(f"Sample standard deviation: {s:.2f}")
print(f"\nt-distribution (correct method):")
print(f"  t_{{{alpha/2}, {n-1}}} = {t_critical:.4f}")
print(f"  95% confidence interval: [{ci_lower_t:.2f}, {ci_upper_t:.2f}]")
print(f"  Interval width: {ci_upper_t - ci_lower_t:.2f}")
print(f"\nz-distribution (incorrect method, for reference):")
print(f"  z_{{{alpha/2}}} = {z_critical:.4f}")
print(f"  95% confidence interval: [{ci_lower_z:.2f}, {ci_upper_z:.2f}]")
print(f"  Interval width: {ci_upper_z - ci_lower_z:.2f}")
print(f"\nRatio of t-interval to z-interval width: {(ci_upper_t - ci_lower_t)/(ci_upper_z - ci_lower_z):.3f}")

# Relationship between sample size and confidence-interval width
sample_sizes = np.array([5, 10, 20, 30, 50, 100, 200])
ci_widths_t = []
ci_widths_z = []

np.random.seed(42)
for ns in sample_sizes:
    sample = np.random.normal(mu_true, sigma_true, ns)
    xb = np.mean(sample)
    sb = np.std(sample, ddof=1)

    t_crit = stats.t.ppf(1 - alpha/2, df=ns-1)
    z_crit = stats.norm.ppf(1 - alpha/2)

    width_t = 2 * t_crit * sb / np.sqrt(ns)
    width_z = 2 * z_crit * sb / np.sqrt(ns)

    ci_widths_t.append(width_t)
    ci_widths_z.append(width_z)

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

# Comparison of t-distribution and z-distribution
x = np.linspace(-4, 4, 200)
axes[0].plot(x, stats.norm.pdf(x), 'b-', linewidth=2, label='Standard normal distribution (z)')
for df in [2, 5, 10, 30]:
    axes[0].plot(x, stats.t.pdf(x, df), linewidth=2, label=f't-distribution (df={df})')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Probability density')
axes[0].set_title('Comparison of the t-distribution and the Standard Normal Distribution')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Sample size and confidence-interval width
axes[1].plot(sample_sizes, ci_widths_t, 'ro-', markersize=8,
             linewidth=2, label='t-interval (correct)')
axes[1].plot(sample_sizes, ci_widths_z, 'b^--', markersize=8,
             linewidth=2, label='z-interval (case of known variance)')
axes[1].set_xlabel('Sample size')
axes[1].set_ylabel('Confidence-interval width')
axes[1].set_title('Sample Size and 95% Confidence-Interval Width')
axes[1].set_xscale('log')
axes[1].legend()
axes[1].grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.show()

📌 Key Points
For small samples (n < 30 or so), the t-distribution must be used. The smaller the degrees of freedom, the heavier the tails of the t-distribution compared with the normal distribution, and the wider the confidence interval. As the sample size increases, the t-distribution approaches the normal distribution.

2.3 Confidence Interval for the Population Variance (Chi-Squared Distribution)

📘 Confidence Interval for the Population Variance

When the population is $ N(\mu, \sigma^2) $, since $ \frac{(n-1)S^2}{\sigma^2} \sim \chi^2_{n-1} $:

$$ \left[ \frac{(n-1)S^2}{\chi^2_{\alpha/2, n-1}}, \frac{(n-1)S^2}{\chi^2_{1-\alpha/2, n-1}} \right] $$

is the $ 1-\alpha $ confidence interval for the population variance $ \sigma^2 $.

💻 Code Example 3: Confidence Interval for the Population Variance

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

# True parameters
mu_true = 50
sigma2_true = 100  # population variance
sigma_true = np.sqrt(sigma2_true)

# Data generation
np.random.seed(42)
n = 20
data = np.random.normal(mu_true, sigma_true, n)

# Sample variance (unbiased)
s2 = np.var(data, ddof=1)

# 95% confidence interval
alpha = 0.05
chi2_lower = stats.chi2.ppf(alpha/2, df=n-1)
chi2_upper = stats.chi2.ppf(1 - alpha/2, df=n-1)

ci_lower = (n-1) * s2 / chi2_upper
ci_upper = (n-1) * s2 / chi2_lower

print("=== Confidence Interval for the Population Variance ===")
print(f"Sample size: {n}")
print(f"Sample variance: {s2:.2f}")
print(f"True population variance: {sigma2_true}")
print(f"\nChi-squared quantiles:")
print(f"  χ²_{{{alpha/2}, {n-1}}} = {chi2_lower:.4f}")
print(f"  χ²_{{{1-alpha/2}, {n-1}}} = {chi2_upper:.4f}")
print(f"\n95% confidence interval for the population variance σ²: [{ci_lower:.2f}, {ci_upper:.2f}]")
print(f"95% confidence interval for the population standard deviation σ: [{np.sqrt(ci_lower):.2f}, {np.sqrt(ci_upper):.2f}]")
print(f"True population variance is inside the interval: {ci_lower <= sigma2_true <= ci_upper}")

# Verify coverage probability via simulation
n_simulations = 1000
coverage_count = 0

np.random.seed(123)
for _ in range(n_simulations):
    sample = np.random.normal(mu_true, sigma_true, n)
    s2_sample = np.var(sample, ddof=1)
    ci_l = (n-1) * s2_sample / chi2_upper
    ci_u = (n-1) * s2_sample / chi2_lower
    if ci_l <= sigma2_true <= ci_u:
        coverage_count += 1

coverage_rate = coverage_count / n_simulations
print(f"\nSimulation ({n_simulations} runs):")
print(f"Coverage probability: {coverage_rate:.3f} (theoretical value: 0.95)")

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

# Chi-squared distribution
x = np.linspace(0, 40, 300)
axes[0].plot(x, stats.chi2.pdf(x, df=n-1), 'b-', linewidth=2,
             label=f'χ² distribution (df={n-1})')
axes[0].axvline(chi2_lower, color='red', linestyle='--', linewidth=2,
                label=f'Lower quantile: {chi2_lower:.2f}')
axes[0].axvline(chi2_upper, color='green', linestyle='--', linewidth=2,
                label=f'Upper quantile: {chi2_upper:.2f}')
axes[0].fill_between(x, 0, stats.chi2.pdf(x, df=n-1),
                      where=(x >= chi2_lower) & (x <= chi2_upper),
                      alpha=0.3, color='yellow', label='95% region')
axes[0].set_xlabel('χ² value')
axes[0].set_ylabel('Probability density')
axes[0].set_title(f'Chi-Squared Distribution (df={n-1})')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Distribution of the sample variance
sample_variances = []
np.random.seed(42)
for _ in range(5000):
    sample = np.random.normal(mu_true, sigma_true, n)
    sample_variances.append(np.var(sample, ddof=1))

axes[1].hist(sample_variances, bins=50, density=True, alpha=0.7,
             color='skyblue', edgecolor='black', label='Distribution of sample variance')
axes[1].axvline(sigma2_true, color='red', linestyle='--', linewidth=2,
                label=f'True population variance: {sigma2_true}')
axes[1].axvline(ci_lower, color='green', linestyle=':', linewidth=2,
                label=f'Lower CI bound: {ci_lower:.1f}')
axes[1].axvline(ci_upper, color='green', linestyle=':', linewidth=2,
                label=f'Upper CI bound: {ci_upper:.1f}')
axes[1].axvspan(ci_lower, ci_upper, alpha=0.2, color='green')
axes[1].set_xlabel('Variance')
axes[1].set_ylabel('Density')
axes[1].set_title('Distribution of the Sample Variance and the Confidence Interval')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

2.4 The Two-Sample Problem and Confidence Interval for the Difference

📘 Confidence Interval for the Difference of Two Population Means

Equal-variance case (using the pooled variance):

$$ (\bar{X}_1 - \bar{X}_2) \pm t_{\alpha/2, n_1+n_2-2} \cdot S_p \sqrt{\frac{1}{n_1} + \frac{1}{n_2}} $$

where $ S_p^2 = \frac{(n_1-1)S_1^2 + (n_2-1)S_2^2}{n_1+n_2-2} $ is the pooled variance.

Unequal-variance case (Welch's t-test):

$$ (\bar{X}_1 - \bar{X}_2) \pm t_{\alpha/2, \nu} \sqrt{\frac{S_1^2}{n_1} + \frac{S_2^2}{n_2}} $$

The degrees of freedom $ \nu $ are computed with the Welch–Satterthwaite approximation.

💻 Code Example 4: Two-Sample t Confidence Interval

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

# Strength data for two manufacturing processes (simulation)
np.random.seed(42)

# Process A
mu_A = 450
sigma_A = 25
n_A = 20
data_A = np.random.normal(mu_A, sigma_A, n_A)

# Process B
mu_B = 470
sigma_B = 30
n_B = 25
data_B = np.random.normal(mu_B, sigma_B, n_B)

# Sample statistics
x_bar_A = np.mean(data_A)
x_bar_B = np.mean(data_B)
s_A = np.std(data_A, ddof=1)
s_B = np.std(data_B, ddof=1)

diff = x_bar_B - x_bar_A
true_diff = mu_B - mu_A

print("=== Two-Sample Problem: Confidence Interval for the Difference of Means ===")
print(f"Process A: n={n_A}, mean={x_bar_A:.2f}, SD={s_A:.2f}")
print(f"Process B: n={n_B}, mean={x_bar_B:.2f}, SD={s_B:.2f}")
print(f"Difference of means: {diff:.2f} (true difference: {true_diff})")

# Test for equality of variances (F-test)
f_stat = s_A**2 / s_B**2 if s_A > s_B else s_B**2 / s_A**2
f_pvalue = 2 * min(stats.f.cdf(s_A**2/s_B**2, n_A-1, n_B-1),
                   1 - stats.f.cdf(s_A**2/s_B**2, n_A-1, n_B-1))
print(f"\nTest for equality of variances (F-test):")
print(f"  F statistic: {f_stat:.4f}, p-value: {f_pvalue:.4f}")

# Method 1: assume equal variances (pooled variance)
alpha = 0.05
s_pooled = np.sqrt(((n_A-1)*s_A**2 + (n_B-1)*s_B**2) / (n_A + n_B - 2))
se_pooled = s_pooled * np.sqrt(1/n_A + 1/n_B)
df_pooled = n_A + n_B - 2
t_crit_pooled = stats.t.ppf(1 - alpha/2, df=df_pooled)
ci_lower_pooled = diff - t_crit_pooled * se_pooled
ci_upper_pooled = diff + t_crit_pooled * se_pooled

print(f"\n[Method 1] Assuming equal variances:")
print(f"  Pooled standard deviation: {s_pooled:.2f}")
print(f"  Degrees of freedom: {df_pooled}")
print(f"  95% confidence interval: [{ci_lower_pooled:.2f}, {ci_upper_pooled:.2f}]")

# Method 2: unequal variances (Welch's method)
se_welch = np.sqrt(s_A**2/n_A + s_B**2/n_B)
df_welch = (s_A**2/n_A + s_B**2/n_B)**2 / \
           (s_A**4/(n_A**2*(n_A-1)) + s_B**4/(n_B**2*(n_B-1)))
t_crit_welch = stats.t.ppf(1 - alpha/2, df=df_welch)
ci_lower_welch = diff - t_crit_welch * se_welch
ci_upper_welch = diff + t_crit_welch * se_welch

print(f"\n[Method 2] Unequal variances (Welch):")
print(f"  Degrees of freedom (Welch–Satterthwaite): {df_welch:.2f}")
print(f"  95% confidence interval: [{ci_lower_welch:.2f}, {ci_upper_welch:.2f}]")

# Using SciPy's function
t_stat, p_value = stats.ttest_ind(data_B, data_A, equal_var=False)
ci_scipy = stats.ttest_ind(data_B, data_A, equal_var=False,
                            alternative='two-sided').confidence_interval(0.95)
print(f"\n[SciPy] Unequal-variance t-test:")
print(f"  95% confidence interval: [{ci_scipy.low:.2f}, {ci_scipy.high:.2f}]")

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

# Distribution of the data
axes[0].boxplot([data_A, data_B], labels=['Process A', 'Process B'],
                patch_artist=True,
                boxprops=dict(facecolor='skyblue', alpha=0.7),
                medianprops=dict(color='red', linewidth=2))
axes[0].scatter(np.ones(n_A)*1 + np.random.normal(0, 0.05, n_A),
                data_A, alpha=0.5, color='blue', s=30)
axes[0].scatter(np.ones(n_B)*2 + np.random.normal(0, 0.05, n_B),
                data_B, alpha=0.5, color='green', s=30)
axes[0].set_ylabel('Strength [MPa]')
axes[0].set_title('Strength Data for the Two Processes')
axes[0].grid(True, alpha=0.3, axis='y')

# Confidence interval for the difference of means
methods = ['Equal-variance\nassumption', 'Welch\n(unequal variance)']
diffs = [diff, diff]
errors_lower = [diff - ci_lower_pooled, diff - ci_lower_welch]
errors_upper = [ci_upper_pooled - diff, ci_upper_welch - diff]

axes[1].errorbar(methods, diffs,
                 yerr=[errors_lower, errors_upper],
                 fmt='o', markersize=10, capsize=10, capthick=2,
                 elinewidth=2, color='blue', label='95% confidence interval')
axes[1].axhline(true_diff, color='red', linestyle='--',
                linewidth=2, label=f'True difference: {true_diff}')
axes[1].axhline(0, color='gray', linestyle=':', linewidth=1,
                label='No difference')
axes[1].set_ylabel('Difference of means [MPa]')
axes[1].set_title('95% Confidence Interval for the Difference of Means')
axes[1].legend()
axes[1].grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

2.5 Confidence Interval for a Proportion (Binomial Distribution)

📘 Confidence Interval for a Population Proportion

Normal approximation (when $ np \geq 5 $ and $ n(1-p) \geq 5 $):

$$ \hat{p} \pm z_{\alpha/2} \sqrt{\frac{\hat{p}(1-\hat{p})}{n}} $$

Wilson interval (applicable even for small samples):

$$ \frac{\hat{p} + \frac{z^2}{2n} \pm z\sqrt{\frac{\hat{p}(1-\hat{p})}{n} + \frac{z^2}{4n^2}}}{1 + \frac{z^2}{n}} $$

💻 Code Example 5: Confidence Interval for a Population Proportion

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.stats.proportion import proportion_confint

# Data (example: survey of a product defect rate)
n = 200  # number of products inspected
x = 15   # number of defective products
p_hat = x / n  # sample proportion

print("=== Confidence Interval for a Population Proportion ===")
print(f"Sample size: {n}")
print(f"Number of defectives: {x}")
print(f"Sample proportion: {p_hat:.4f} ({p_hat*100:.2f}%)")

# Method 1: normal approximation (Wald interval)
alpha = 0.05
z = stats.norm.ppf(1 - alpha/2)
se = np.sqrt(p_hat * (1 - p_hat) / n)
ci_lower_wald = p_hat - z * se
ci_upper_wald = p_hat + z * se

print(f"\n[Wald interval] (normal approximation):")
print(f"  95% confidence interval: [{ci_lower_wald:.4f}, {ci_upper_wald:.4f}]")
print(f"               [{ci_lower_wald*100:.2f}%, {ci_upper_wald*100:.2f}%]")

# Method 2: Wilson interval
term1 = p_hat + z**2 / (2*n)
term2 = z * np.sqrt(p_hat*(1-p_hat)/n + z**2/(4*n**2))
denominator = 1 + z**2/n
ci_lower_wilson = (term1 - term2) / denominator
ci_upper_wilson = (term1 + term2) / denominator

print(f"\n[Wilson interval]:")
print(f"  95% confidence interval: [{ci_lower_wilson:.4f}, {ci_upper_wilson:.4f}]")
print(f"               [{ci_lower_wilson*100:.2f}%, {ci_upper_wilson*100:.2f}%]")

# Method 3: Clopper-Pearson interval (exact)
ci_lower_cp, ci_upper_cp = proportion_confint(x, n, alpha=alpha, method='beta')

print(f"\n[Clopper-Pearson interval] (exact):")
print(f"  95% confidence interval: [{ci_lower_cp:.4f}, {ci_upper_cp:.4f}]")
print(f"               [{ci_lower_cp*100:.2f}%, {ci_upper_cp*100:.2f}%]")

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

# Comparison of confidence intervals across methods
methods = ['Wald\n(normal approx.)', 'Wilson', 'Clopper-\nPearson']
lowers = [ci_lower_wald, ci_lower_wilson, ci_lower_cp]
uppers = [ci_upper_wald, ci_upper_wilson, ci_upper_cp]
centers = [p_hat, p_hat, p_hat]

for i, (method, lower, upper) in enumerate(zip(methods, lowers, uppers)):
    axes[0].plot([lower, upper], [i, i], 'o-', linewidth=3, markersize=8,
                 label=method)

axes[0].axvline(p_hat, color='red', linestyle='--', linewidth=2,
                label=f'Sample proportion: {p_hat:.3f}')
axes[0].set_yticks(range(len(methods)))
axes[0].set_yticklabels(methods)
axes[0].set_xlabel('Proportion')
axes[0].set_title('95% Confidence Interval for the Population Proportion (Method Comparison)')
axes[0].grid(True, alpha=0.3, axis='x')
axes[0].legend()

# Relationship between sample size and confidence-interval width
sample_sizes = np.logspace(1, 3, 50, dtype=int)
widths_wald = []
widths_wilson = []

p_fixed = 0.075  # fixed true proportion

for ns in sample_sizes:
    x_sim = int(p_fixed * ns)
    p_sim = x_sim / ns

    # Wald
    se_sim = np.sqrt(p_sim * (1-p_sim) / ns)
    w_wald = 2 * z * se_sim
    widths_wald.append(w_wald)

    # Wilson
    t1 = p_sim + z**2/(2*ns)
    t2 = z * np.sqrt(p_sim*(1-p_sim)/ns + z**2/(4*ns**2))
    denom = 1 + z**2/ns
    w_wilson = 2 * t2 / denom
    widths_wilson.append(w_wilson)

axes[1].plot(sample_sizes, widths_wald, 'b-', linewidth=2, label='Wald interval')
axes[1].plot(sample_sizes, widths_wilson, 'g--', linewidth=2, label='Wilson interval')
axes[1].set_xlabel('Sample size')
axes[1].set_ylabel('Confidence-interval width')
axes[1].set_title(f'Sample Size vs. Confidence-Interval Width (p={p_fixed})')
axes[1].set_xscale('log')
axes[1].set_yscale('log')
axes[1].legend()
axes[1].grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.show()

2.6 Bootstrap Confidence Intervals

📘 The Bootstrap Method

When the theoretical distribution is unknown, a confidence interval can be constructed by resampling:

Generate bootstrap samples from the original data by sampling with replacement.
Compute the statistic for each bootstrap sample.
Obtain quantiles from the distribution of the statistic and construct the confidence interval.

Percentile method: Take the $ \alpha/2 $ quantile and the $ 1-\alpha/2 $ quantile of the bootstrap distribution as the interval endpoints.

💻 Code Example 6: Bootstrap Confidence Interval

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

# Data generation (skewed distribution)
np.random.seed(42)
# Gamma distribution (right-skewed)
data = np.random.gamma(2, 2, 50)

# Statistic to estimate: the median
observed_median = np.median(data)

print("=== Bootstrap Confidence Interval ===")
print(f"Sample size: {len(data)}")
print(f"Sample median: {observed_median:.2f}")
print(f"Sample mean: {np.mean(data):.2f}")
print(f"Skewness of the data: {stats.skew(data):.2f}")

# Bootstrap resampling
n_bootstrap = 10000
bootstrap_medians = []

np.random.seed(123)
for _ in range(n_bootstrap):
    # Sampling with replacement
    bootstrap_sample = np.random.choice(data, size=len(data), replace=True)
    bootstrap_medians.append(np.median(bootstrap_sample))

bootstrap_medians = np.array(bootstrap_medians)

# Confidence interval by the percentile method
alpha = 0.05
ci_lower_percentile = np.percentile(bootstrap_medians, alpha/2 * 100)
ci_upper_percentile = np.percentile(bootstrap_medians, (1 - alpha/2) * 100)

print(f"\n[Percentile method]:")
print(f"  95% confidence interval: [{ci_lower_percentile:.2f}, {ci_upper_percentile:.2f}]")

# Confidence interval by the BCa method (Bias-Corrected and Accelerated)
# A more accurate method
from scipy.stats import norm

# Compute the bias-correction term
z0 = norm.ppf(np.sum(bootstrap_medians < observed_median) / n_bootstrap)

# Compute the acceleration constant (jackknife method)
jackknife_medians = []
for i in range(len(data)):
    jackknife_sample = np.delete(data, i)
    jackknife_medians.append(np.median(jackknife_sample))

jackknife_medians = np.array(jackknife_medians)
mean_jack = np.mean(jackknife_medians)
numerator = np.sum((mean_jack - jackknife_medians)**3)
denominator = 6 * (np.sum((mean_jack - jackknife_medians)**2))**(3/2)
a = numerator / denominator if denominator != 0 else 0

# Quantiles for the BCa confidence interval
z_alpha = norm.ppf(alpha/2)
z_1alpha = norm.ppf(1 - alpha/2)

p_lower = norm.cdf(z0 + (z0 + z_alpha)/(1 - a*(z0 + z_alpha)))
p_upper = norm.cdf(z0 + (z0 + z_1alpha)/(1 - a*(z0 + z_1alpha)))

ci_lower_bca = np.percentile(bootstrap_medians, p_lower * 100)
ci_upper_bca = np.percentile(bootstrap_medians, p_upper * 100)

print(f"\n[BCa method]:")
print(f"  Bias-correction term z0: {z0:.4f}")
print(f"  Acceleration constant a: {a:.4f}")
print(f"  95% confidence interval: [{ci_lower_bca:.2f}, {ci_upper_bca:.2f}]")

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

# Distribution of the original data
axes[0].hist(data, bins=15, density=True, alpha=0.7,
             color='skyblue', edgecolor='black', label='Data')
axes[0].axvline(observed_median, color='red', linestyle='--',
                linewidth=2, label=f'Sample median: {observed_median:.1f}')
axes[0].axvline(np.mean(data), color='blue', linestyle=':',
                linewidth=2, label=f'Sample mean: {np.mean(data):.1f}')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Density')
axes[0].set_title('Distribution of the Original Data (right-skewed)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Bootstrap distribution
axes[1].hist(bootstrap_medians, bins=50, density=True, alpha=0.7,
             color='lightgreen', edgecolor='black', label='Bootstrap distribution')
axes[1].axvline(observed_median, color='red', linestyle='--',
                linewidth=2, label=f'Observed median: {observed_median:.1f}')
axes[1].axvline(ci_lower_percentile, color='blue', linestyle=':',
                linewidth=2, label='Percentile method')
axes[1].axvline(ci_upper_percentile, color='blue', linestyle=':',
                linewidth=2)
axes[1].axvline(ci_lower_bca, color='purple', linestyle='--',
                linewidth=2, label='BCa method')
axes[1].axvline(ci_upper_bca, color='purple', linestyle='--',
                linewidth=2)
axes[1].set_xlabel('Median')
axes[1].set_ylabel('Density')
axes[1].set_title(f'Bootstrap Distribution ({n_bootstrap} runs)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Comparison with the normal approximation
print(f"\n[Reference] Confidence interval for the median by normal approximation:")
# Approximate standard error of the median: SE ≈ 1.253 * σ / sqrt(n)
se_median_approx = 1.253 * np.std(data) / np.sqrt(len(data))
ci_lower_normal = observed_median - 1.96 * se_median_approx
ci_upper_normal = observed_median + 1.96 * se_median_approx
print(f"  95% confidence interval: [{ci_lower_normal:.2f}, {ci_upper_normal:.2f}]")
print(f"  (may be inaccurate for a skewed distribution)")

📌 Key Points
The bootstrap method is a powerful technique that can construct confidence intervals without assuming the shape of the distribution. It is especially useful for statistics with complicated theoretical distributions, such as the median or skewness. The BCa method is more accurate than the percentile method, correcting for bias and skewness.

2.7 Interval Estimation for Material Strength Data

💻 Code Example 7: Comprehensive Interval Estimation for Material Strength Data

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from statsmodels.stats.proportion import proportion_confint

# Material strength data (simulation of experimental data)
np.random.seed(42)
n_specimens = 30
true_mean = 450  # MPa
true_std = 30    # MPa

# Simulate with a Weibull distribution (brittle material)
k_weibull = 15
lambda_weibull = 470
strength_data = np.random.weibull(k_weibull, n_specimens) * lambda_weibull

print("=== Interval Estimation for Material Strength Data ===")
print(f"Number of specimens: {n_specimens}")
print(f"Mean strength: {np.mean(strength_data):.2f} MPa")
print(f"Standard deviation: {np.std(strength_data, ddof=1):.2f} MPa")
print(f"Minimum: {np.min(strength_data):.2f} MPa")
print(f"Maximum: {np.max(strength_data):.2f} MPa")

# 1. 95% confidence interval for the mean strength (t-interval)
mean_strength = np.mean(strength_data)
std_strength = np.std(strength_data, ddof=1)
se_mean = std_strength / np.sqrt(n_specimens)
alpha = 0.05
t_crit = stats.t.ppf(1 - alpha/2, df=n_specimens-1)
ci_mean_lower = mean_strength - t_crit * se_mean
ci_mean_upper = mean_strength + t_crit * se_mean

print(f"\n[95% confidence interval for the mean strength]:")
print(f"  [{ci_mean_lower:.2f}, {ci_mean_upper:.2f}] MPa")

# 2. 95% confidence interval for the standard deviation (chi-squared distribution)
chi2_lower = stats.chi2.ppf(alpha/2, df=n_specimens-1)
chi2_upper = stats.chi2.ppf(1-alpha/2, df=n_specimens-1)
ci_std_lower = np.sqrt((n_specimens-1) * std_strength**2 / chi2_upper)
ci_std_upper = np.sqrt((n_specimens-1) * std_strength**2 / chi2_lower)

print(f"\n[95% confidence interval for the standard deviation]:")
print(f"  [{ci_std_lower:.2f}, {ci_std_upper:.2f}] MPa")

# 3. Confidence interval for the coefficient of variation (CV) (bootstrap)
cv_observed = std_strength / mean_strength

n_bootstrap = 5000
cv_bootstrap = []
np.random.seed(123)
for _ in range(n_bootstrap):
    boot_sample = np.random.choice(strength_data, size=n_specimens, replace=True)
    cv_bootstrap.append(np.std(boot_sample, ddof=1) / np.mean(boot_sample))

ci_cv_lower = np.percentile(cv_bootstrap, alpha/2 * 100)
ci_cv_upper = np.percentile(cv_bootstrap, (1-alpha/2) * 100)

print(f"\n[95% confidence interval for the coefficient of variation (CV)]:")
print(f"  Observed CV: {cv_observed:.4f} ({cv_observed*100:.2f}%)")
print(f"  95% confidence interval: [{ci_cv_lower:.4f}, {ci_cv_upper:.4f}]")
print(f"               [{ci_cv_lower*100:.2f}%, {ci_cv_upper*100:.2f}%]")

# 4. Confidence interval for the failure probability
# Estimate the probability of failure at or below 400 MPa
threshold = 400  # MPa
n_failed = np.sum(strength_data <= threshold)
p_failure = n_failed / n_specimens

print(f"\n[95% confidence interval for the failure probability] (at or below {threshold} MPa):")
print(f"  Number of failed specimens: {n_failed}/{n_specimens}")
print(f"  Sample failure probability: {p_failure:.4f} ({p_failure*100:.2f}%)")

if n_failed > 0:
    ci_fail_lower, ci_fail_upper = proportion_confint(n_failed, n_specimens,
                                                       alpha=alpha, method='wilson')
    print(f"  Wilson confidence interval: [{ci_fail_lower:.4f}, {ci_fail_upper:.4f}]")
    print(f"                   [{ci_fail_lower*100:.2f}%, {ci_fail_upper*100:.2f}%]")
else:
    print(f"  (No failures: only an upper confidence limit can be computed)")

# 5. Confidence interval for specific percentiles (bootstrap)
# Confidence interval for the 5th percentile (lowest 5% of strength)
percentile_5 = np.percentile(strength_data, 5)
percentile_95 = np.percentile(strength_data, 95)

p5_bootstrap = []
p95_bootstrap = []
for _ in range(n_bootstrap):
    boot_sample = np.random.choice(strength_data, size=n_specimens, replace=True)
    p5_bootstrap.append(np.percentile(boot_sample, 5))
    p95_bootstrap.append(np.percentile(boot_sample, 95))

ci_p5_lower = np.percentile(p5_bootstrap, alpha/2 * 100)
ci_p5_upper = np.percentile(p5_bootstrap, (1-alpha/2) * 100)
ci_p95_lower = np.percentile(p95_bootstrap, alpha/2 * 100)
ci_p95_upper = np.percentile(p95_bootstrap, (1-alpha/2) * 100)

print(f"\n[95% confidence interval for the 5th percentile]:")
print(f"  Observed value: {percentile_5:.2f} MPa")
print(f"  95% confidence interval: [{ci_p5_lower:.2f}, {ci_p5_upper:.2f}] MPa")

print(f"\n[95% confidence interval for the 95th percentile]:")
print(f"  Observed value: {percentile_95:.2f} MPa")
print(f"  95% confidence interval: [{ci_p95_lower:.2f}, {ci_p95_upper:.2f}] MPa")

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

# 1. Data distribution and confidence interval for the mean
ax1 = fig.add_subplot(gs[0, :])
ax1.hist(strength_data, bins=12, density=True, alpha=0.7,
         color='skyblue', edgecolor='black', label='Experimental data')
ax1.axvline(mean_strength, color='blue', linestyle='--', linewidth=2,
            label=f'Mean: {mean_strength:.1f} MPa')
ax1.axvspan(ci_mean_lower, ci_mean_upper, alpha=0.3, color='blue',
            label=f'95% CI for the mean: [{ci_mean_lower:.1f}, {ci_mean_upper:.1f}]')
ax1.axvline(threshold, color='red', linestyle=':', linewidth=2,
            label=f'Design criterion: {threshold} MPa')
ax1.set_xlabel('Strength [MPa]')
ax1.set_ylabel('Density')
ax1.set_title('Material Strength Data and Statistical Interval Estimation')
ax1.legend()
ax1.grid(True, alpha=0.3)

# 2. Bootstrap distribution of the CV
ax2 = fig.add_subplot(gs[1, 0])
ax2.hist(cv_bootstrap, bins=40, density=True, alpha=0.7,
         color='lightgreen', edgecolor='black')
ax2.axvline(cv_observed, color='red', linestyle='--', linewidth=2,
            label=f'Observed CV: {cv_observed:.3f}')
ax2.axvline(ci_cv_lower, color='blue', linestyle=':', linewidth=2)
ax2.axvline(ci_cv_upper, color='blue', linestyle=':', linewidth=2,
            label=f'95% CI: [{ci_cv_lower:.3f}, {ci_cv_upper:.3f}]')
ax2.set_xlabel('Coefficient of variation (CV)')
ax2.set_ylabel('Density')
ax2.set_title('Bootstrap Distribution of the Coefficient of Variation')
ax2.legend()
ax2.grid(True, alpha=0.3)

# 3. Confidence intervals for percentiles
ax3 = fig.add_subplot(gs[1, 1])
percentiles = [5, 25, 50, 75, 95]
obs_values = [np.percentile(strength_data, p) for p in percentiles]
ci_lowers = [ci_p5_lower, 0, 0, 0, ci_p95_lower]
ci_uppers = [ci_p5_upper, 0, 0, 0, ci_p95_upper]

# Show confidence intervals only for the 5th and 95th percentiles
ax3.errorbar([0, 4], [obs_values[0], obs_values[4]],
             yerr=[[obs_values[0]-ci_lowers[0], obs_values[4]-ci_lowers[4]],
                   [ci_uppers[0]-obs_values[0], ci_uppers[4]-obs_values[4]]],
             fmt='o', markersize=10, capsize=8, capthick=2,
             elinewidth=2, color='blue', label='95% confidence interval')
ax3.plot(range(5), obs_values, 'ro-', markersize=8, linewidth=2,
         label='Observed values')
ax3.set_xticks(range(5))
ax3.set_xticklabels([f'{p}%' for p in percentiles])
ax3.set_ylabel('Strength [MPa]')
ax3.set_title('Percentile Values and Confidence Intervals')
ax3.legend()
ax3.grid(True, alpha=0.3, axis='y')

# 4. Comparison of confidence-interval widths
ax4 = fig.add_subplot(gs[2, :])
metrics = ['Mean', 'Std. dev.', 'CV', '5th pct.', '95th pct.']
observed = [mean_strength, std_strength, cv_observed, percentile_5, percentile_95]
lower_ci = [ci_mean_lower, ci_std_lower, ci_cv_lower, ci_p5_lower, ci_p95_lower]
upper_ci = [ci_mean_upper, ci_std_upper, ci_cv_upper, ci_p5_upper, ci_p95_upper]

# Normalize for comparison
relative_widths = [(u-l)/o for o, l, u in zip(observed, lower_ci, upper_ci)]

ax4.bar(metrics, relative_widths, color='skyblue', edgecolor='black', alpha=0.7)
ax4.set_ylabel('Relative interval width (width / estimate)')
ax4.set_title('Relative Confidence-Interval Width of Each Statistic (Comparison of Uncertainty)')
ax4.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.show()

# Practical summary
print("\n=== Practical Interpretation ===")
print(f"✓ The mean strength lies within {ci_mean_lower:.1f}–{ci_mean_upper:.1f} MPa with 95% confidence")
print(f"✓ The data spread (standard deviation) is {ci_std_lower:.1f}–{ci_std_upper:.1f} MPa")
print(f"✓ The coefficient of variation is {ci_cv_lower*100:.1f}%–{ci_cv_upper*100:.1f}% (indicator of quality stability)")
print(f"✓ The lower-5% material strength is {ci_p5_lower:.1f}–{ci_p5_upper:.1f} MPa (reference for the design allowable value)")

📝 Exercises

Compute how the width of a 95% confidence interval at the same confidence level changes between a sample size of 10 and a sample size of 100.
Check how the width of the confidence interval changes as you vary the confidence level to 90%, 95%, and 99%.
Implement a confidence interval for paired two samples (paired t-test).
Construct a confidence interval for the sample mean with the bootstrap method and compare it with the theoretical t-interval.

Summary

Confidence intervals are an important technique for quantifying the uncertainty of a parameter.
Use the z-distribution when the population variance is known, and the t-distribution when it is unknown.
The t-distribution is essential for small samples; for large samples it approaches the normal distribution.
Use the chi-squared distribution for a confidence interval on the population variance.
In the two-sample problem, choose the method according to the equality of variances (pooled variance vs. Welch's method).
The Wilson interval is recommended for a confidence interval on a population proportion (applicable even for small samples).
The bootstrap method is powerful and requires no distributional assumptions.
In materials science, confidence intervals for multiple statistics are evaluated together.