Chapter 1: Foundations of Estimation Theory | Inferential Statistics and Bayesian Statistics

1.1 Overview of Estimation Theory

Estimation Theory is the foundation of statistics for inferring the parameters of a population from sample data. In materials science, it is essential in situations where the properties of an entire material (mean strength, variance) are estimated from a small number of specimen data.

📘 Basic Concepts of Estimation Theory

Population: The entire set of objects under investigation

Sample: A subset of data drawn from the population

Estimator: A function computed from the sample that gives an estimate of a parameter

Estimate: The specific numerical value obtained by substituting the sample data into the estimator

For example, when the sample mean $ \bar{X} = \frac{1}{n}\sum_{i=1}^n X_i $ is used as an estimator of the population mean $ \mu $, the value $ \bar{x} $ computed from the actual data $ x_1, \ldots, x_n $ is the estimate.

1.2 Point Estimation and Properties of Estimators

1.2.1 Unbiasedness

📘 Definition of an Unbiased Estimator

An estimator $ \hat{\theta} $ is an unbiased estimator of the parameter $ \theta $ when:

$$ E[\hat{\theta}] = \theta $$

That is, the expected value of the estimator coincides with the true parameter.

💻 Code Example 1: Verifying the unbiasedness of the sample mean and sample variance

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

# True population parameters
mu_true = 50  # population mean
sigma_true = 10  # population standard deviation

# Simulation settings
n_samples = 30  # sample size
n_simulations = 10000  # number of simulations

# Record the distributions of the estimators
sample_means = []
sample_vars_biased = []  # biased (divide by n)
sample_vars_unbiased = []  # unbiased estimator (divide by n-1)

np.random.seed(42)
for _ in range(n_simulations):
    # Draw a sample
    sample = np.random.normal(mu_true, sigma_true, n_samples)

    # Sample mean
    sample_means.append(np.mean(sample))

    # Sample variance (biased)
    sample_vars_biased.append(np.var(sample, ddof=0))

    # Sample variance (unbiased estimator)
    sample_vars_unbiased.append(np.var(sample, ddof=1))

# Verify the results
print("=== Verification of unbiasedness ===")
print(f"Population mean: {mu_true}")
print(f"Expected value of sample mean: {np.mean(sample_means):.4f}")
print(f"Standard error of sample mean: {np.std(sample_means):.4f}")
print()
print(f"Population variance: {sigma_true**2}")
print(f"Expected value of biased sample variance: {np.mean(sample_vars_biased):.4f}")
print(f"Expected value of unbiased sample variance: {np.mean(sample_vars_unbiased):.4f}")

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

# Distribution of the sample mean
axes[0].hist(sample_means, bins=50, density=True, alpha=0.7,
             color='skyblue', edgecolor='black')
axes[0].axvline(mu_true, color='red', linestyle='--', linewidth=2,
                label=f'True population mean: {mu_true}')
axes[0].axvline(np.mean(sample_means), color='blue', linestyle='--',
                linewidth=2, label=f'Expected value of estimator: {np.mean(sample_means):.2f}')
axes[0].set_xlabel('Sample mean')
axes[0].set_ylabel('Density')
axes[0].set_title('Distribution of the sample mean (verifying unbiasedness)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Distribution of the sample variance
axes[1].hist(sample_vars_biased, bins=50, density=True, alpha=0.5,
             color='orange', edgecolor='black', label='Biased (ddof=0)')
axes[1].hist(sample_vars_unbiased, bins=50, density=True, alpha=0.5,
             color='green', edgecolor='black', label='Unbiased estimator (ddof=1)')
axes[1].axvline(sigma_true**2, color='red', linestyle='--',
                linewidth=2, label=f'True population variance: {sigma_true**2}')
axes[1].set_xlabel('Sample variance')
axes[1].set_ylabel('Density')
axes[1].set_title('Distribution of the sample variance (comparison of unbiasedness)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Comparison with the theoretical standard error
theoretical_se = sigma_true / np.sqrt(n_samples)
print(f"\nTheoretical standard error: {theoretical_se:.4f}")
print(f"Empirical standard error: {np.std(sample_means):.4f}")

📌 Key Point
The sample mean is an unbiased estimator of the population mean, but the sample variance must be divided by $ n-1 $ (Bessel's correction). This corrects for the tendency of the sample variance to underestimate.

1.2.2 Consistency

📘 Definition of a Consistent Estimator

An estimator $ \hat{\theta}_n $ is a consistent estimator when, as the sample size $ n \to \infty $:

$$ \hat{\theta}_n \xrightarrow{P} \theta $$

That is, it converges in probability to the true parameter as the sample size grows.

1.2.3 Efficiency

When multiple unbiased estimators exist, the one with the minimum variance is the most efficient. The Cramér-Rao lower bound gives the theoretical lower limit that the variance of an unbiased estimator can attain.

1.3 Maximum Likelihood Estimation

📘 Principle of Maximum Likelihood Estimation

Given observed data $ x_1, \ldots, x_n $, this is a method that estimates the parameter $ \theta $ under which these data are most likely to occur.

Likelihood Function:

$$ L(\theta; x_1, \ldots, x_n) = \prod_{i=1}^n f(x_i; \theta) $$

Log-Likelihood Function:

$$ \ell(\theta) = \log L(\theta) = \sum_{i=1}^n \log f(x_i; \theta) $$

The maximum likelihood estimator (MLE) is the $ \hat{\theta}_{\text{MLE}} $ that maximizes $ \ell(\theta) $.

💻 Code Example 2: Maximum likelihood estimation of normal distribution parameters

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

# True parameters
mu_true = 100
sigma_true = 15

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

# Log-likelihood function (normal distribution)
def neg_log_likelihood(params, data):
    """Negative log-likelihood (for minimization)"""
    mu, sigma = params
    if sigma <= 0:
        return np.inf
    n = len(data)
    return n/2 * np.log(2*np.pi*sigma**2) + np.sum((data - mu)**2) / (2*sigma**2)

# Maximum likelihood estimation (numerical optimization)
initial_guess = [np.mean(data), np.std(data)]
result = minimize(neg_log_likelihood, initial_guess, args=(data,),
                  method='Nelder-Mead')

mu_mle, sigma_mle = result.x

print("=== Maximum likelihood estimation results ===")
print(f"True population mean: {mu_true}, MLE: {mu_mle:.4f}")
print(f"True population std. dev.: {sigma_true}, MLE: {sigma_mle:.4f}")
print(f"Sample mean (analytical solution): {np.mean(data):.4f}")
print(f"Sample std. dev. (analytical solution, ddof=0): {np.std(data, ddof=0):.4f}")

# Visualization of the likelihood function
mu_grid = np.linspace(90, 110, 100)
sigma_grid = np.linspace(10, 20, 100)
MU, SIGMA = np.meshgrid(mu_grid, sigma_grid)

# Computation of the log-likelihood
log_likelihood = np.zeros_like(MU)
for i in range(len(mu_grid)):
    for j in range(len(sigma_grid)):
        log_likelihood[j, i] = -neg_log_likelihood([MU[j, i], SIGMA[j, i]], data)

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

# Contour plot of the log-likelihood
contour = axes[0].contourf(MU, SIGMA, log_likelihood, levels=30, cmap='viridis')
axes[0].plot(mu_true, sigma_true, 'r*', markersize=15, label='True value')
axes[0].plot(mu_mle, sigma_mle, 'wo', markersize=10, label='MLE')
axes[0].set_xlabel('μ (mean)')
axes[0].set_ylabel('σ (standard deviation)')
axes[0].set_title('Contour plot of the log-likelihood function')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
plt.colorbar(contour, ax=axes[0], label='Log-likelihood')

# Histogram of the data and the estimated distribution
axes[1].hist(data, bins=15, density=True, alpha=0.7,
             color='skyblue', edgecolor='black', label='Data')
x_range = np.linspace(data.min(), data.max(), 200)
axes[1].plot(x_range, stats.norm.pdf(x_range, mu_true, sigma_true),
             'r-', linewidth=2, label=f'True distribution N({mu_true}, {sigma_true}²)')
axes[1].plot(x_range, stats.norm.pdf(x_range, mu_mle, sigma_mle),
             'b--', linewidth=2, label=f'MLE distribution N({mu_mle:.1f}, {sigma_mle:.1f}²)')
axes[1].set_xlabel('Value')
axes[1].set_ylabel('Probability density')
axes[1].set_title('Data and the estimated distribution')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

📌 Key Point
For the normal distribution, the analytical solution of the MLE coincides with the sample mean and the sample variance (divided by $ n $). For general distributions, numerical optimization is required.

💻 Code Example 3: Maximum likelihood estimation of the binomial and Poisson distributions

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

# === Maximum likelihood estimation of the binomial distribution ===
print("=== Maximum likelihood estimation of the binomial distribution ===")
# Data generation (n=10 trials performed 50 times)
n_trials = 10
p_true = 0.3
np.random.seed(42)
data_binomial = np.random.binomial(n_trials, p_true, 50)

# Analytical MLE: p_hat = (total successes) / (total trials)
p_mle = np.sum(data_binomial) / (len(data_binomial) * n_trials)
print(f"True p: {p_true}, MLE: {p_mle:.4f}")

# === Maximum likelihood estimation of the Poisson distribution ===
print("\n=== Maximum likelihood estimation of the Poisson distribution ===")
# Data generation
lambda_true = 5.0
data_poisson = np.random.poisson(lambda_true, 100)

# Analytical MLE: lambda_hat = sample mean
lambda_mle = np.mean(data_poisson)
print(f"True λ: {lambda_true}, MLE: {lambda_mle:.4f}")

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

# Binomial distribution
x_binom = np.arange(0, n_trials+1)
axes[0].hist(data_binomial, bins=np.arange(-0.5, n_trials+1.5, 1),
             density=True, alpha=0.7, color='skyblue',
             edgecolor='black', label='Data')
axes[0].plot(x_binom, stats.binom.pmf(x_binom, n_trials, p_true),
             'ro-', markersize=8, linewidth=2, label=f'True distribution B({n_trials}, {p_true})')
axes[0].plot(x_binom, stats.binom.pmf(x_binom, n_trials, p_mle),
             'b^--', markersize=6, linewidth=2, label=f'MLE distribution B({n_trials}, {p_mle:.2f})')
axes[0].set_xlabel('Number of successes')
axes[0].set_ylabel('Probability mass')
axes[0].set_title('Maximum likelihood estimation of the binomial distribution')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Poisson distribution
x_poisson = np.arange(0, max(data_poisson)+1)
axes[1].hist(data_poisson, bins=np.arange(-0.5, max(data_poisson)+1.5, 1),
             density=True, alpha=0.7, color='lightgreen',
             edgecolor='black', label='Data')
axes[1].plot(x_poisson, stats.poisson.pmf(x_poisson, lambda_true),
             'ro-', markersize=6, linewidth=2, label=f'True distribution Poisson({lambda_true})')
axes[1].plot(x_poisson, stats.poisson.pmf(x_poisson, lambda_mle),
             'b^--', markersize=5, linewidth=2, label=f'MLE distribution Poisson({lambda_mle:.2f})')
axes[1].set_xlabel('Number of occurrences')
axes[1].set_ylabel('Probability mass')
axes[1].set_title('Maximum likelihood estimation of the Poisson distribution')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

1.4 Method of Moments

📘 Principle of the Method of Moments

Parameters are estimated by equating the population moments with the sample moments.

k-th moment:

$$ m_k = E[X^k], \quad \hat{m}_k = \frac{1}{n}\sum_{i=1}^n X_i^k $$

We set up and solve as many moment equations as there are parameters.

💻 Code Example 4: Estimation by the method of moments (gamma distribution)

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

# True parameters (gamma distribution)
# Gamma distribution: Gamma(k, θ), E[X] = kθ, Var[X] = kθ²
k_true = 3.0  # shape parameter
theta_true = 2.0  # scale parameter

# Data generation
np.random.seed(42)
n = 200
data = np.random.gamma(k_true, theta_true, n)

# Estimation by the method of moments
# From E[X] = kθ, Var[X] = kθ²
# m1 = kθ, m2 - m1² = kθ²
m1 = np.mean(data)
m2 = np.mean(data**2)
var_sample = m2 - m1**2

# Simultaneous equations: m1 = kθ, var = kθ²
# Since var/m1 = θ, θ_hat = var/m1
# k_hat = m1/θ_hat = m1²/var
theta_mom = var_sample / m1
k_mom = m1**2 / var_sample

print("=== Estimation by the method of moments (gamma distribution) ===")
print(f"True shape parameter k: {k_true}")
print(f"Method-of-moments estimate k: {k_mom:.4f}")
print(f"True scale parameter θ: {theta_true}")
print(f"Method-of-moments estimate θ: {theta_mom:.4f}")

# Comparison with maximum likelihood estimation (requires numerical solution)
def gamma_mle_equations(params, data):
    """Equations for the MLE of the gamma distribution"""
    k, theta = params
    n = len(data)
    eq1 = n * np.log(theta) + np.sum(np.log(data)) - n * (np.log(k) + stats.digamma(k))
    eq2 = np.sum(data) - n * k * theta
    return [eq1, eq2]

# Use the method-of-moments result as the initial value
initial_guess = [k_mom, theta_mom]
k_mle, theta_mle = fsolve(gamma_mle_equations, initial_guess, args=(data,))

print(f"\nMaximum likelihood estimate k: {k_mle:.4f}")
print(f"Maximum likelihood estimate θ: {theta_mle:.4f}")

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

# Histogram of the data and the estimated distributions
x_range = np.linspace(0, data.max(), 200)
axes[0].hist(data, bins=30, density=True, alpha=0.7,
             color='skyblue', edgecolor='black', label='Data')
axes[0].plot(x_range, stats.gamma.pdf(x_range, k_true, scale=theta_true),
             'r-', linewidth=2.5, label=f'True distribution Γ({k_true}, {theta_true})')
axes[0].plot(x_range, stats.gamma.pdf(x_range, k_mom, scale=theta_mom),
             'g--', linewidth=2, label=f'Method of moments Γ({k_mom:.2f}, {theta_mom:.2f})')
axes[0].plot(x_range, stats.gamma.pdf(x_range, k_mle, scale=theta_mle),
             'b:', linewidth=2, label=f'MLE Γ({k_mle:.2f}, {theta_mle:.2f})')
axes[0].set_xlabel('Value')
axes[0].set_ylabel('Probability density')
axes[0].set_title('Estimation of the gamma distribution (method of moments vs MLE)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Q-Q plot (visual assessment of estimation accuracy)
theoretical_quantiles = np.linspace(0.01, 0.99, 100)
data_sorted = np.sort(data)
empirical_quantiles = np.linspace(0, 1, len(data_sorted))

for method, k, theta, color, linestyle in [
    ('True distribution', k_true, theta_true, 'red', '-'),
    ('Method of moments', k_mom, theta_mom, 'green', '--'),
    ('MLE', k_mle, theta_mle, 'blue', ':')
]:
    theoretical_values = stats.gamma.ppf(theoretical_quantiles, k, scale=theta)
    axes[1].plot(theoretical_values,
                 np.quantile(data, theoretical_quantiles),
                 color=color, linestyle=linestyle, linewidth=2, label=method)

axes[1].plot([0, data.max()], [0, data.max()], 'k--', alpha=0.5, label='y=x')
axes[1].set_xlabel('Theoretical quantiles')
axes[1].set_ylabel('Sample quantiles')
axes[1].set_title('Q-Q plot (comparison of estimation accuracy)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

1.5 Bias-Variance Trade-off of Estimators

📘 Decomposition of the Mean Squared Error (MSE)

The mean squared error of an estimator $ \hat{\theta} $ is:

$$ \text{MSE}(\hat{\theta}) = E[(\hat{\theta} - \theta)^2] = \text{Bias}^2(\hat{\theta}) + \text{Var}(\hat{\theta}) $$

where:

Bias: $ \text{Bias}(\hat{\theta}) = E[\hat{\theta}] - \theta $
Variance: $ \text{Var}(\hat{\theta}) = E[(\hat{\theta} - E[\hat{\theta}])^2] $

Even an unbiased estimator can have a large MSE if its variance is large.

💻 Code Example 5: Visualizing the bias-variance trade-off

import numpy as np
import matplotlib.pyplot as plt

# True population mean and population variance
mu_true = 100
sigma_true = 20

# Sample size
n = 20

# Simulation
np.random.seed(42)
n_simulations = 5000

# Three types of estimators
# 1. Sample mean (unbiased)
# 2. Constant estimator (biased, zero variance)
# 3. Shrinkage estimator (biased, small variance)

estimates_unbiased = []
estimates_constant = []
estimates_shrinkage = []

constant_value = 95  # guess based on prior knowledge
shrinkage_factor = 0.7  # shrinkage coefficient

for _ in range(n_simulations):
    sample = np.random.normal(mu_true, sigma_true, n)

    # Unbiased estimator (sample mean)
    estimates_unbiased.append(np.mean(sample))

    # Constant estimator
    estimates_constant.append(constant_value)

    # Shrinkage estimator: shrink the sample mean toward the constant value
    estimates_shrinkage.append(
        shrinkage_factor * np.mean(sample) + (1 - shrinkage_factor) * constant_value
    )

# Computation of the MSE
def compute_mse_components(estimates, true_value):
    estimates = np.array(estimates)
    bias = np.mean(estimates) - true_value
    variance = np.var(estimates)
    mse = np.mean((estimates - true_value)**2)
    return bias, variance, mse

bias_u, var_u, mse_u = compute_mse_components(estimates_unbiased, mu_true)
bias_c, var_c, mse_c = compute_mse_components(estimates_constant, mu_true)
bias_s, var_s, mse_s = compute_mse_components(estimates_shrinkage, mu_true)

print("=== Bias-variance trade-off ===")
print(f"\nUnbiased estimator (sample mean):")
print(f"  Bias: {bias_u:.4f}, Variance: {var_u:.4f}, MSE: {mse_u:.4f}")
print(f"\nConstant estimator:")
print(f"  Bias: {bias_c:.4f}, Variance: {var_c:.4f}, MSE: {mse_c:.4f}")
print(f"\nShrinkage estimator:")
print(f"  Bias: {bias_s:.4f}, Variance: {var_s:.4f}, MSE: {mse_s:.4f}")

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

# Distributions of the estimators
axes[0, 0].hist(estimates_unbiased, bins=50, density=True, alpha=0.6,
                color='blue', edgecolor='black', label='Unbiased estimator')
axes[0, 0].hist(estimates_shrinkage, bins=50, density=True, alpha=0.6,
                color='green', edgecolor='black', label='Shrinkage estimator')
axes[0, 0].axvline(mu_true, color='red', linestyle='--', linewidth=2,
                   label=f'True value: {mu_true}')
axes[0, 0].axvline(constant_value, color='orange', linestyle='--',
                   linewidth=2, label=f'Constant value: {constant_value}')
axes[0, 0].set_xlabel('Estimate')
axes[0, 0].set_ylabel('Density')
axes[0, 0].set_title('Distributions of the estimators')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Bar chart of the MSE decomposition
methods = ['Unbiased estimator', 'Constant estimator', 'Shrinkage estimator']
biases_sq = [bias_u**2, bias_c**2, bias_s**2]
variances = [var_u, var_c, var_s]
x_pos = np.arange(len(methods))

axes[0, 1].bar(x_pos, biases_sq, width=0.35, label='Bias²',
               color='orange', alpha=0.8)
axes[0, 1].bar(x_pos, variances, width=0.35, bottom=biases_sq,
               label='Variance', color='skyblue', alpha=0.8)
axes[0, 1].set_xticks(x_pos)
axes[0, 1].set_xticklabels(methods, rotation=15, ha='right')
axes[0, 1].set_ylabel('Value')
axes[0, 1].set_title('Decomposition of the MSE (Bias² + Variance)')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3, axis='y')

# Distribution of the squared error
sq_errors_u = (np.array(estimates_unbiased) - mu_true)**2
sq_errors_s = (np.array(estimates_shrinkage) - mu_true)**2

axes[1, 0].hist(sq_errors_u, bins=50, density=True, alpha=0.6,
                color='blue', edgecolor='black', label='Unbiased estimator')
axes[1, 0].hist(sq_errors_s, bins=50, density=True, alpha=0.6,
                color='green', edgecolor='black', label='Shrinkage estimator')
axes[1, 0].axvline(mse_u, color='blue', linestyle='--', linewidth=2)
axes[1, 0].axvline(mse_s, color='green', linestyle='--', linewidth=2)
axes[1, 0].set_xlabel('Squared error')
axes[1, 0].set_ylabel('Density')
axes[1, 0].set_title('Distribution of the squared error')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Relationship between the shrinkage coefficient and the MSE
shrinkage_factors = np.linspace(0, 1, 50)
mse_by_shrinkage = []

for sf in shrinkage_factors:
    est = [sf * e + (1-sf) * constant_value for e in estimates_unbiased]
    _, _, mse = compute_mse_components(est, mu_true)
    mse_by_shrinkage.append(mse)

axes[1, 1].plot(shrinkage_factors, mse_by_shrinkage, 'b-', linewidth=2)
axes[1, 1].axvline(shrinkage_factor, color='green', linestyle='--',
                   linewidth=2, label=f'Chosen coefficient: {shrinkage_factor}')
axes[1, 1].axhline(mse_u, color='blue', linestyle=':', linewidth=2,
                   label=f'MSE of unbiased estimator: {mse_u:.2f}')
optimal_sf = shrinkage_factors[np.argmin(mse_by_shrinkage)]
axes[1, 1].axvline(optimal_sf, color='red', linestyle='--',
                   linewidth=2, label=f'Optimal coefficient: {optimal_sf:.2f}')
axes[1, 1].set_xlabel('Shrinkage coefficient')
axes[1, 1].set_ylabel('MSE')
axes[1, 1].set_title('Relationship between the shrinkage coefficient and the MSE')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\nOptimal shrinkage coefficient: {optimal_sf:.4f}")

📌 Key Point
An unbiased estimator is not always the best. When allowing a little bias can greatly reduce the variance, the MSE becomes smaller (as with the James-Stein estimator, etc.).

1.6 Fisher Information and the Cramér-Rao Lower Bound

📘 Fisher Information

The Fisher information of a parameter $ \theta $ is:

$$ I(\theta) = E\left[\left(\frac{\partial \log f(X;\theta)}{\partial \theta}\right)^2\right] = -E\left[\frac{\partial^2 \log f(X;\theta)}{\partial \theta^2}\right] $$

It represents the amount of information the data carry about the parameter.

📘 Cramér-Rao Lower Bound

The variance of any unbiased estimator $ \hat{\theta} $ satisfies:

$$ \text{Var}(\hat{\theta}) \geq \frac{1}{nI(\theta)} $$

An estimator that attains this lower bound is an efficient estimator.

💻 Code Example 6: Computing the Fisher information and the Cramér-Rao lower bound

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

# Fisher information of the normal distribution N(μ, σ²)
# I(μ) = 1/σ², I(σ²) = 1/(2σ⁴)

def normal_fisher_info_mean(sigma):
    """Fisher information about the mean of the normal distribution"""
    return 1 / sigma**2

def normal_fisher_info_var(sigma):
    """Fisher information about the variance of the normal distribution"""
    return 1 / (2 * sigma**4)

# Parameter settings
mu_true = 50
sigma_true = 10
sample_sizes = np.array([10, 20, 50, 100, 200, 500])

# Simulation
n_simulations = 10000
np.random.seed(42)

results = []
for n in sample_sizes:
    sample_means = []
    sample_vars = []

    for _ in range(n_simulations):
        sample = np.random.normal(mu_true, sigma_true, n)
        sample_means.append(np.mean(sample))
        sample_vars.append(np.var(sample, ddof=1))

    # Empirical variance
    empirical_var_mean = np.var(sample_means)
    empirical_var_var = np.var(sample_vars)

    # Cramér-Rao lower bound
    cr_bound_mean = 1 / (n * normal_fisher_info_mean(sigma_true))
    cr_bound_var = 1 / (n * normal_fisher_info_var(sigma_true))

    results.append({
        'n': n,
        'empirical_var_mean': empirical_var_mean,
        'cr_bound_mean': cr_bound_mean,
        'empirical_var_var': empirical_var_var,
        'cr_bound_var': cr_bound_var
    })

# Display the results
print("=== Fisher information and the Cramér-Rao lower bound ===")
print(f"Normal distribution N({mu_true}, {sigma_true}²)")
print(f"\nFisher information (per sample):")
print(f"  I(μ) = {normal_fisher_info_mean(sigma_true):.6f}")
print(f"  I(σ²) = {normal_fisher_info_var(sigma_true):.10f}")
print()

for r in results:
    print(f"n={r['n']}:")
    print(f"  Estimator of the mean: empirical variance={r['empirical_var_mean']:.4f}, "
          f"CR bound={r['cr_bound_mean']:.4f}, "
          f"efficiency={(r['cr_bound_mean']/r['empirical_var_mean'])*100:.2f}%")
    print(f"  Estimator of the variance: empirical variance={r['empirical_var_var']:.4f}, "
          f"CR bound={r['cr_bound_var']:.4f}, "
          f"efficiency={(r['cr_bound_var']/r['empirical_var_var'])*100:.2f}%")
    print()

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

# Estimator of the mean
empirical_vars_mean = [r['empirical_var_mean'] for r in results]
cr_bounds_mean = [r['cr_bound_mean'] for r in results]

axes[0].plot(sample_sizes, empirical_vars_mean, 'bo-',
             markersize=8, linewidth=2, label='Empirical variance')
axes[0].plot(sample_sizes, cr_bounds_mean, 'r^--',
             markersize=8, linewidth=2, label='Cramér-Rao lower bound')
axes[0].set_xlabel('Sample size n')
axes[0].set_ylabel('Variance')
axes[0].set_title('Variance of the sample mean vs Cramér-Rao lower bound')
axes[0].set_xscale('log')
axes[0].set_yscale('log')
axes[0].legend()
axes[0].grid(True, alpha=0.3, which='both')

# Estimator of the variance
empirical_vars_var = [r['empirical_var_var'] for r in results]
cr_bounds_var = [r['cr_bound_var'] for r in results]

axes[1].plot(sample_sizes, empirical_vars_var, 'bo-',
             markersize=8, linewidth=2, label='Empirical variance')
axes[1].plot(sample_sizes, cr_bounds_var, 'r^--',
             markersize=8, linewidth=2, label='Cramér-Rao lower bound')
axes[1].set_xlabel('Sample size n')
axes[1].set_ylabel('Variance')
axes[1].set_title('Variance of the sample variance vs Cramér-Rao lower bound')
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()

📌 Key Point
For the normal distribution, the sample mean attains the Cramér-Rao lower bound (100% efficient), whereas the sample variance attains the bound only asymptotically.

1.7 Estimation of Material Property Data

💻 Code Example 7: Statistical estimation of material strength data

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

# Generation of experimental data (in practice, obtained from experiments)
# The distribution of material strength is often modeled by a normal or Weibull distribution
np.random.seed(42)

# Scenario 1: modeled by a normal distribution
n_samples = 30
true_mean_strength = 450  # MPa
true_std_strength = 30  # MPa
strength_data_normal = np.random.normal(true_mean_strength, true_std_strength, n_samples)

# Scenario 2: modeled by a Weibull distribution (suitable for brittle materials)
# Weibull distribution: parameters (k: shape, λ: scale)
k_true = 15  # shape parameter (larger means smaller variance)
lambda_true = 470  # scale parameter
strength_data_weibull = np.random.weibull(k_true, n_samples) * lambda_true

print("=== Statistical estimation of material strength data ===")
print("\n[Normal distribution model]")
# Estimation for the normal distribution
mean_est = np.mean(strength_data_normal)
std_est = np.std(strength_data_normal, ddof=1)

print(f"Sample mean (MLE): {mean_est:.2f} MPa (true value: {true_mean_strength} MPa)")
print(f"Sample std. dev.: {std_est:.2f} MPa (true value: {true_std_strength} MPa)")

# 95% confidence interval (covered in detail in the next chapter)
se = std_est / np.sqrt(n_samples)
ci_95 = stats.t.interval(0.95, n_samples-1, loc=mean_est, scale=se)
print(f"95% confidence interval for the mean strength: [{ci_95[0]:.2f}, {ci_95[1]:.2f}] MPa")

print("\n[Weibull distribution model]")
# MLE for the Weibull distribution (using the SciPy function)
# scipy.stats.weibull_min.fit() returns (c, loc, scale)
# c = k (shape), scale = λ (scale)
params = stats.weibull_min.fit(strength_data_weibull, floc=0)
k_est, loc_est, lambda_est = params

print(f"Shape parameter k (MLE): {k_est:.2f} (true value: {k_true})")
print(f"Scale parameter λ (MLE): {lambda_est:.2f} (true value: {lambda_true})")
print(f"Mean strength (Weibull): {lambda_est * stats.gamma(1 + 1/k_est):.2f} MPa")

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

# Normal distribution model
axes[0, 0].hist(strength_data_normal, bins=12, density=True, alpha=0.7,
                color='skyblue', edgecolor='black', label='Experimental data')
x_range = np.linspace(strength_data_normal.min(), strength_data_normal.max(), 200)
axes[0, 0].plot(x_range, stats.norm.pdf(x_range, true_mean_strength, true_std_strength),
                'r-', linewidth=2.5, label=f'True distribution N({true_mean_strength}, {true_std_strength}²)')
axes[0, 0].plot(x_range, stats.norm.pdf(x_range, mean_est, std_est),
                'b--', linewidth=2, label=f'Estimated distribution N({mean_est:.1f}, {std_est:.1f}²)')
axes[0, 0].set_xlabel('Strength [MPa]')
axes[0, 0].set_ylabel('Probability density')
axes[0, 0].set_title('Normal distribution model: estimation of material strength')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# Q-Q plot for the normal distribution
stats.probplot(strength_data_normal, dist="norm", plot=axes[0, 1])
axes[0, 1].set_title('Normal Q-Q plot (checking normality)')
axes[0, 1].grid(True, alpha=0.3)

# Weibull distribution model
axes[1, 0].hist(strength_data_weibull, bins=12, density=True, alpha=0.7,
                color='lightgreen', edgecolor='black', label='Experimental data')
x_range_w = np.linspace(0, strength_data_weibull.max(), 200)
axes[1, 0].plot(x_range_w, stats.weibull_min.pdf(x_range_w, k_true, scale=lambda_true),
                'r-', linewidth=2.5, label=f'True distribution Weibull({k_true}, {lambda_true})')
axes[1, 0].plot(x_range_w, stats.weibull_min.pdf(x_range_w, k_est, scale=lambda_est),
                'b--', linewidth=2, label=f'Estimated distribution Weibull({k_est:.1f}, {lambda_est:.1f})')
axes[1, 0].set_xlabel('Strength [MPa]')
axes[1, 0].set_ylabel('Probability density')
axes[1, 0].set_title('Weibull distribution model: estimation of material strength')
axes[1, 0].legend()
axes[1, 0].grid(True, alpha=0.3)

# Weibull probability plot
sorted_data = np.sort(strength_data_weibull)
n = len(sorted_data)
empirical_cdf = np.arange(1, n+1) / (n+1)
weibull_y = np.log(-np.log(1 - empirical_cdf))
weibull_x = np.log(sorted_data)

axes[1, 1].plot(weibull_x, weibull_y, 'bo', markersize=6, label='Experimental data')
# Theoretical line of the estimated Weibull distribution
x_fit = np.array([weibull_x.min(), weibull_x.max()])
y_fit = k_est * (x_fit - np.log(lambda_est))
axes[1, 1].plot(x_fit, y_fit, 'r-', linewidth=2, label='Estimated Weibull line')
axes[1, 1].set_xlabel('ln(strength)')
axes[1, 1].set_ylabel('ln(-ln(1-F))')
axes[1, 1].set_title('Weibull probability plot (checking goodness of fit)')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Estimation of the failure probability (important in engineering)
print("\n[Estimation of the failure probability]")
critical_strength = 400  # MPa
prob_failure_normal = stats.norm.cdf(critical_strength, mean_est, std_est)
prob_failure_weibull = stats.weibull_min.cdf(critical_strength, k_est, scale=lambda_est)

print(f"Probability of failure at or below 400 MPa:")
print(f"  Normal distribution model: {prob_failure_normal*100:.2f}%")
print(f"  Weibull model: {prob_failure_weibull*100:.2f}%")

📌 Practical Point
In the analysis of material strength data:

Ductile materials (metals, etc.) → normal or log-normal distribution
Brittle materials (ceramics, etc.) → Weibull distribution
Check the goodness of fit of the distribution with Q-Q plots and probability plots
Use the estimation results for failure probability and reliability evaluation

📝 Exercises

Investigate whether the sample median is an unbiased estimator of the population mean (in the case of a normal distribution).
Derive the maximum likelihood estimator of the parameter $ \lambda $ of the exponential distribution $ f(x; \lambda) = \lambda e^{-\lambda x} $.
Show that $ \max(X_1, \ldots, X_n) $ is a consistent estimator of the parameter $ \theta $ of the uniform distribution $ U(0, \theta) $.
Compute the Fisher information about the parameter $ p $ of the Bernoulli distribution.

Summary

Estimation theory is the foundation of statistics for inferring population parameters from a sample
Unbiasedness, consistency, and efficiency are important properties of estimators
Maximum likelihood estimation is the most general and powerful estimation method
The method of moments is easy to compute and useful for obtaining initial estimates
The bias-variance trade-off is an important consideration in the design of estimators
Fisher information and the Cramér-Rao lower bound give the theoretical limits of estimation
In materials science, estimation using the normal and Weibull distributions is important