Chapter 2: Law of Large Numbers and Central Limit Theorem

Law of Large Numbers and Central Limit Theorem

📖 Reading Time: 20-25 minutes 📊 Difficulty: Beginner 💻 Code Examples: 0 📝 Exercises: 0

🌐 EN | 🇯🇵 JP | Last sync: 2025-11-16

Fundamentals of Mathematics Dojo > Probability Theory and Stochastic Processes > Chapter 2

2.1 Weak Law of Large Numbers

📊 Theorem: Weak Law of Large Numbers
For independent and identically distributed (i.i.d.) random variables \(X_1, X_2, \ldots, X_n\) with expectation \(E[X_i] = \mu\) and variance \(Var(X_i) = \sigma^2 < \infty\), the sample mean \(\bar{X}n = \frac{1}{n}\sum{i=1}^n X_i\) converges in probability to \(\mu\): \[\lim_{n \to \infty} P(|\bar{X}_n - \mu| > \epsilon) = 0 \quad \text{for any } \epsilon > 0\] Proof (using Chebyshev’s inequality): \[P(|\bar{X}_n - \mu| > \epsilon) \leq \frac{Var(\bar{X}_n)}{\epsilon^2} = \frac{\sigma^2}{n\epsilon^2} \to 0 \quad (n \to \infty)\]

💻 Code Example 1: Simulation of Weak Law of Large Numbers

import numpy as np import matplotlib.pyplot as plt from scipy import stats # Parameter settings np.random.seed(42) mu = 5 # True expected value sigma = 2 # True standard deviation max_n = 10000 epsilon = 0.5 # Tolerance # Sampling (from normal distribution) samples = np.random.normal(mu, sigma, max_n) # Calculate cumulative mean cumulative_mean = np.cumsum(samples) / np.arange(1, max_n + 1) # Estimate probability of |X̄_n - μ| > ε n_values = np.logspace(1, 4, 100).astype(int) n_values = np.unique(n_values) # Remove duplicates prob_exceed = [] for n in n_values: # Multiple simulations n_trials = 1000 exceed_count = 0 for _ in range(n_trials): sample_mean = np.random.normal(mu, sigma, n).mean() if abs(sample_mean - mu) > epsilon: exceed_count += 1 prob_exceed.append(exceed_count / n_trials) # Theoretical value (upper bound from Chebyshev’s inequality) theoretical_bound = sigma2 / (n_values * epsilon2) # Visualization fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # (1) Convergence of sample mean axes[0, 0].plot(cumulative_mean, color=‘#667eea’, linewidth=1.5, alpha=0.8) axes[0, 0].axhline(y=mu, color=‘red’, linestyle=’—’, linewidth=2, label=f’True mean μ={mu}’) axes[0, 0].axhline(y=mu+epsilon, color=‘orange’, linestyle=’:’, alpha=0.7, label=f’μ±ε (ε={epsilon})’) axes[0, 0].axhline(y=mu-epsilon, color=‘orange’, linestyle=’:’, alpha=0.7) axes[0, 0].set_xlabel(‘Sample size n’, fontsize=11) axes[0, 0].set_ylabel(‘Sample mean X̄ₙ’, fontsize=11) axes[0, 0].set_title(‘Convergence of Sample Mean’, fontsize=12, fontweight=‘bold’) axes[0, 0].legend() axes[0, 0].grid(alpha=0.3) axes[0, 0].set_xlim([0, max_n]) # (2) Convergence on log scale axes[0, 1].semilogx(cumulative_mean, color=‘#667eea’, linewidth=1.5, alpha=0.8) axes[0, 1].axhline(y=mu, color=‘red’, linestyle=’—’, linewidth=2) axes[0, 1].axhline(y=mu+epsilon, color=‘orange’, linestyle=’:’, alpha=0.7) axes[0, 1].axhline(y=mu-epsilon, color=‘orange’, linestyle=’:’, alpha=0.7) axes[0, 1].set_xlabel(‘Sample size n (log)’, fontsize=11) axes[0, 1].set_ylabel(‘Sample mean X̄ₙ’, fontsize=11) axes[0, 1].set_title(‘Convergence of Sample Mean (Log Scale)’, fontsize=12, fontweight=‘bold’) axes[0, 1].grid(alpha=0.3, which=‘both’) # (3) Estimate of P(|X̄ₙ - μ| > ε) axes[1, 0].loglog(n_values, prob_exceed, ‘o-’, color=‘#667eea’, linewidth=2, markersize=4, label=‘Simulation’) axes[1, 0].loglog(n_values, theoretical_bound, ’—’, color=‘red’, linewidth=2, label=‘Chebyshev upper bound’) axes[1, 0].set_xlabel(‘Sample size n’, fontsize=11) axes[1, 0].set_ylabel(‘P(|X̄ₙ - μ| > ε)’, fontsize=11) axes[1, 0].set_title(‘Verification of Weak Law of Large Numbers’, fontsize=12, fontweight=‘bold’) axes[1, 0].legend() axes[1, 0].grid(alpha=0.3, which=‘both’) # (4) Comparison across different distributions distributions = { ‘Normal’: lambda n: np.random.normal(mu, sigma, n), ‘Uniform’: lambda n: np.random.uniform(mu-sigmanp.sqrt(3), mu+sigmanp.sqrt(3), n), ‘Exponential’: lambda n: np.random.exponential(mu, n), } for name, dist_func in distributions.items(): samples_dist = dist_func(max_n) cumulative_mean_dist = np.cumsum(samples_dist) / np.arange(1, max_n + 1) axes[1, 1].plot(cumulative_mean_dist, linewidth=1.5, alpha=0.7, label=name) axes[1, 1].axhline(y=mu, color=‘red’, linestyle=’—’, linewidth=2, label=f’μ={mu}’) axes[1, 1].set_xlabel(‘Sample size n’, fontsize=11) axes[1, 1].set_ylabel(‘Sample mean X̄ₙ’, fontsize=11) axes[1, 1].set_title(‘Law of Large Numbers for Different Distributions’, fontsize=12, fontweight=‘bold’) axes[1, 1].legend(fontsize=9) axes[1, 1].grid(alpha=0.3) axes[1, 1].set_xlim([0, 5000]) plt.tight_layout() plt.show() print(“Verification of Weak Law of Large Numbers:”) print(f”True expected value: μ = {mu}”) print(f”Sample mean at n=100: {cumulative_mean[99]:.4f}”) print(f”Sample mean at n=1000: {cumulative_mean[999]:.4f}”) print(f”Sample mean at n=10000: {cumulative_mean[9999]:.4f}”)

Verification of Weak Law of Large Numbers: True expected value: μ = 5 Sample mean at n=100: 5.0234 Sample mean at n=1000: 4.9876 Sample mean at n=10000: 5.0012

2.2 Strong Law of Large Numbers

📊 Theorem: Strong Law of Large Numbers
For independent and identically distributed random variables \(X_1, X_2, \ldots\), if \(E[|X_i|] < \infty\), then the sample mean converges almost surely to the expectation: \[P\left(\lim_{n \to \infty} \bar{X}_n = \mu\right) = 1\] Strong law vs Weak law:

Almost sure convergence is a stronger condition than convergence in probability.

💻 Code Example 2: Demonstration of Strong Law of Large Numbers

Simulation with multiple independent paths np.random.seed(123) n_paths = 50 # Number of paths max_n = 5000 mu = 10 sigma = 3 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) # (1) Sample mean trajectories for multiple paths for i in range(n_paths): samples = np.random.normal(mu, sigma, max_n) cumulative_mean = np.cumsum(samples) / np.arange(1, max_n + 1) axes[0, 0].plot(cumulative_mean, alpha=0.3, linewidth=0.8, color=‘#667eea’) axes[0, 0].axhline(y=mu, color=‘red’, linestyle=’—’, linewidth=2.5, label=f’μ={mu}’) axes[0, 0].set_xlabel(‘Sample size n’, fontsize=11) axes[0, 0].set_ylabel(‘Sample mean X̄ₙ’, fontsize=11) axes[0, 0].set_title(f’{n_paths} Independent Paths (Strong Law of Large Numbers)’, fontsize=12, fontweight=‘bold’) axes[0, 0].legend() axes[0, 0].grid(alpha=0.3) axes[0, 0].set_ylim([mu-2, mu+2]) # (2) Decay of maximum deviation n_checkpoints = np.array([10, 50, 100, 500, 1000, 2000, 5000]) max_deviations = [] for n in n_checkpoints: deviations = [] for _ in range(1000): samples = np.random.normal(mu, sigma, n) sample_mean = samples.mean() deviations.append(abs(sample_mean - mu)) max_deviations.append(np.max(deviations)) axes[0, 1].loglog(n_checkpoints, max_deviations, ‘o-’, color=‘#764ba2’, linewidth=2, markersize=8, label=‘Maximum deviation’) axes[0, 1].loglog(n_checkpoints, sigma / np.sqrt(n_checkpoints), ’—’, color=‘red’, linewidth=2, label=‘σ/√n (theoretical)’) axes[0, 1].set_xlabel(‘Sample size n’, fontsize=11) axes[0, 1].set_ylabel(‘max |X̄ₙ - μ|’, fontsize=11) axes[0, 1].set_title(‘Decay of Maximum Deviation’, fontsize=12, fontweight=‘bold’) axes[0, 1].legend() axes[0, 1].grid(alpha=0.3, which=‘both’) # (3) Confidence bands n_trials = 500 all_paths = [] for _ in range(n_trials): samples = np.random.normal(mu, sigma, max_n) cumulative_mean = np.cumsum(samples) / np.arange(1, max_n + 1) all_paths.append(cumulative_mean) all_paths = np.array(all_paths) percentile_5 = np.percentile(all_paths, 5, axis=0) percentile_95 = np.percentile(all_paths, 95, axis=0) median = np.percentile(all_paths, 50, axis=0) x_axis = np.arange(1, max_n + 1) axes[1, 0].fill_between(x_axis, percentile_5, percentile_95, alpha=0.3, color=‘#667eea’, label=‘90% confidence band’) axes[1, 0].plot(median, color=‘#667eea’, linewidth=2, label=‘Median’) axes[1, 0].axhline(y=mu, color=‘red’, linestyle=’—’, linewidth=2, label=f’μ={mu}’) axes[1, 0].set_xlabel(‘Sample size n’, fontsize=11) axes[1, 0].set_ylabel(‘Sample mean X̄ₙ’, fontsize=11) axes[1, 0].set_title(‘Confidence Band of Convergence (90%)’, fontsize=12, fontweight=‘bold’) axes[1, 0].legend() axes[1, 0].grid(alpha=0.3) axes[1, 0].set_ylim([mu-1, mu+1]) # (4) Decay of variance variance_n = np.var(all_paths, axis=0) theoretical_var = sigma**2 / x_axis axes[1, 1].loglog(x_axis[::10], variance_n[::10], ‘o’, color=‘#667eea’, markersize=3, alpha=0.6, label=‘Empirical variance’) axes[1, 1].loglog(x_axis, theoretical_var, ’—’, color=‘red’, linewidth=2, label=‘σ²/n (theoretical)’) axes[1, 1].set_xlabel(‘Sample size n’, fontsize=11) axes[1, 1].set_ylabel(‘Var(X̄ₙ)’, fontsize=11) axes[1, 1].set_title(‘Decay of Sample Mean Variance’, fontsize=12, fontweight=‘bold’) axes[1, 1].legend() axes[1, 1].grid(alpha=0.3, which=‘both’) plt.tight_layout() plt.show() print(“Demonstration of Strong Law of Large Numbers:”) print(f”Range of sample means at n=5000 over 500 independent trials:”) print(f” Minimum: {all_paths[:, -1].min():.4f}”) print(f” Maximum: {all_paths[:, -1].max():.4f}”) print(f” Median: {np.median(all_paths[:, -1]):.4f}”) print(f” Standard deviation: {np.std(all_paths[:, -1]):.4f}“)

2.3 Central Limit Theorem

📊 Theorem: Central Limit Theorem (CLT)
For independent and identically distributed random variables \(X_1, X_2, \ldots, X_n\) with expectation \(E[X_i] = \mu\) and variance \(Var(X_i) = \sigma^2 < \infty\), the standardized sample mean converges in distribution to the standard normal distribution: \[Z_n = \frac{\bar{X}n - \mu}{\sigma / \sqrt{n}} = \frac{\sum{i=1}^n X_i - n\mu}{\sigma\sqrt{n}} \xrightarrow{d} N(0, 1)\] That is: \[\lim_{n \to \infty} P(Z_n \leq z) = \Phi(z) = \int_{-\infty}^z \frac{1}{\sqrt{2\pi}} e^{-t^2/2} dt\] Significance: Regardless of the original distribution, the sample mean approaches a normal distribution.

💻 Code Example 3: Visualization of Central Limit Theorem

Examine the distribution of sample means from various distributions np.random.seed(42) # Original distributions distributions = { ‘Uniform U(0,1)’: { ‘generator’: lambda n: np.random.uniform(0, 1, n), ‘mu’: 0.5, ‘sigma’: 1/np.sqrt(12) }, ‘Bernoulli (p=0.3)’: { ‘generator’: lambda n: np.random.binomial(1, 0.3, n), ‘mu’: 0.3, ‘sigma’: np.sqrt(0.3 * 0.7) }, ‘Exponential Exp(1)’: { ‘generator’: lambda n: np.random.exponential(1, n), ‘mu’: 1, ‘sigma’: 1 }, ‘Beta Beta(0.5,0.5)’: { ‘generator’: lambda n: np.random.beta(0.5, 0.5, n), ‘mu’: 0.5, ‘sigma’: 1/(2*np.sqrt(2)) } } sample_sizes = [1, 5, 30, 100] n_trials = 10000 fig, axes = plt.subplots(len(distributions), len(sample_sizes), figsize=(16, 12)) for i, (dist_name, dist_info) in enumerate(distributions.items()): generator = dist_info[‘generator’] mu = dist_info[‘mu’] sigma = dist_info[‘sigma’] for j, n in enumerate(sample_sizes): # Calculate sample means sample_means = [] for _ in range(n_trials): samples = generator(n) sample_means.append(samples.mean()) sample_means = np.array(sample_means) # Standardize z_scores = (sample_means - mu) / (sigma / np.sqrt(n)) # Histogram axes[i, j].hist(z_scores, bins=50, density=True, alpha=0.6, color=‘#667eea’, edgecolor=‘black’, linewidth=0.5) # Overlay standard normal distribution x = np.linspace(-4, 4, 1000) axes[i, j].plot(x, stats.norm.pdf(x), ‘r-’, linewidth=2.5, label=‘N(0,1)’) axes[i, j].set_xlim([-4, 4]) axes[i, j].set_ylim([0, 0.5]) if j == 0: axes[i, j].set_ylabel(dist_name, fontsize=10, fontweight=‘bold’) if i == 0: axes[i, j].set_title(f’n={n}’, fontsize=11, fontweight=‘bold’) if i == len(distributions) - 1: axes[i, j].set_xlabel(‘Z = (X̄ₙ - μ)/(σ/√n)’, fontsize=9) axes[i, j].grid(alpha=0.3) axes[i, j].legend(fontsize=8) plt.suptitle(‘Visualization of Central Limit Theorem: Verification Across Various Distributions’, fontsize=14, fontweight=‘bold’, y=0.995) plt.tight_layout() plt.show() print(“Verification of Central Limit Theorem:”) print(”=“*60) for dist_name, dist_info in distributions.items(): print(f”\n{dist_name}:”) for n in [5, 30, 100]: sample_means = [dist_info’generator’.mean() for _ in range(n_trials)] z_scores = (np.array(sample_means) - dist_info[‘mu’]) / (dist_info[‘sigma’] / np.sqrt(n)) print(f” n={n:3d}: mean={np.mean(z_scores):6.3f}, variance={np.var(z_scores):6.3f}“)

2.4 Convergence from Non-Normal to Normal Distributions

💻 Code Example 4: Convergence from Non-Normal to Normal Distributions

Central limit theorem for extremely non-normal distributions np.random.seed(42) # Extremely skewed distribution (log-normal distribution) mu_log = 0 sigma_log = 1 def lognormal_generator(n): return np.random.lognormal(mu_log, sigma_log, n) # Theoretical statistics for log-normal distribution true_mean = np.exp(mu_log + sigma_log2 / 2) true_var = (np.exp(sigma_log2) - 1) * np.exp(2*mu_log + sigma_log**2) true_sigma = np.sqrt(true_var) sample_sizes = [1, 2, 5, 10, 30, 100] n_trials = 5000 fig, axes = plt.subplots(2, 3, figsize=(15, 10)) axes = axes.flatten() for idx, n in enumerate(sample_sizes): # Calculate sample means sample_means = [] for _ in range(n_trials): samples = lognormal_generator(n) sample_means.append(samples.mean()) sample_means = np.array(sample_means) # Standardize z_scores = (sample_means - true_mean) / (true_sigma / np.sqrt(n)) # Data for Q-Q plot theoretical_quantiles = stats.norm.ppf(np.linspace(0.01, 0.99, 100)) sample_quantiles = np.percentile(z_scores, np.linspace(1, 99, 100)) # Histogram + Q-Q plot ax_main = axes[idx] # Histogram ax_main.hist(z_scores, bins=40, density=True, alpha=0.5, color=‘#667eea’, edgecolor=‘black’, linewidth=0.5) # Theoretical normal distribution x = np.linspace(-4, 4, 1000) ax_main.plot(x, stats.norm.pdf(x), ‘r-’, linewidth=2.5, label=‘N(0,1)’) ax_main.set_xlim([-4, 4]) ax_main.set_title(f’n={n}’, fontsize=12, fontweight=‘bold’) ax_main.set_xlabel(‘Standardized sample mean’, fontsize=10) ax_main.set_ylabel(‘Density’, fontsize=10) ax_main.legend(fontsize=9) ax_main.grid(alpha=0.3) plt.suptitle(‘Central Limit Theorem for Log-Normal Distribution (Extremely Skewed)’, fontsize=14, fontweight=‘bold’) plt.tight_layout() plt.show() # Test normality with Kolmogorov-Smirnov test print(“Verification of Central Limit Theorem for Log-Normal Distribution:”) print(”=“*60) print(f”Skewness of original distribution: {stats.lognorm.stats(s=sigma_log, moments=’s’):.4f}”) print(f”Kurtosis of original distribution: {stats.lognorm.stats(s=sigma_log, moments=‘k’):.4f}”) print(“\nKolmogorov-Smirnov test (normality test):”) for n in sample_sizes: sample_means = [lognormal_generator(n).mean() for _ in range(n_trials)] z_scores = (np.array(sample_means) - true_mean) / (true_sigma / np.sqrt(n)) ks_stat, p_value = stats.kstest(z_scores, ‘norm’) print(f” n={n:3d}: KS statistic={ks_stat:.4f}, p-value={p_value:.4f} ” f”{‘(follows normal distribution)’ if p_value > 0.05 else ‘(significantly different)’}“)

2.5 Distribution of Sample Mean and Sample Variance

📊 Theorem: Distributions of Sample Statistics
When \(X_1, \ldots, X_n \sim N(\mu, \sigma^2)\) (i.i.d.): Sample mean: \[\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i \sim N\left(\mu, \frac{\sigma^2}{n}\right)\] Sample variance: \[S^2 = \frac{1}{n-1}\sum_{i=1}^n (X_i - \bar{X})^2, \quad \frac{(n-1)S^2}{\sigma^2} \sim \chi^2_{n-1}\] t-statistic: \[T = \frac{\bar{X} - \mu}{S / \sqrt{n}} \sim t_{n-1}\]

💻 Code Example 5: Distribution of Sample Mean

Distribution of sample mean np.random.seed(42) mu = 50 sigma = 10 sample_sizes = [5, 10, 30, 100] n_trials = 10000 fig, axes = plt.subplots(2, 2, figsize=(14, 10)) axes = axes.flatten() for idx, n in enumerate(sample_sizes): # Calculate sample means sample_means = [] for _ in range(n_trials): samples = np.random.normal(mu, sigma, n) sample_means.append(samples.mean()) sample_means = np.array(sample_means) # Histogram axes[idx].hist(sample_means, bins=50, density=True, alpha=0.6, color=‘#667eea’, edgecolor=‘black’, linewidth=0.5, label=‘Simulation’) # Theoretical distribution N(μ, σ²/n) x = np.linspace(mu - 4sigma/np.sqrt(n), mu + 4sigma/np.sqrt(n), 1000) theoretical_pdf = stats.norm.pdf(x, mu, sigma/np.sqrt(n)) axes[idx].plot(x, theoretical_pdf, ‘r-’, linewidth=2.5, label=f’N({mu}, {sigma2/n:.2f})’) axes[idx].axvline(mu, color=‘green’, linestyle=’—’, linewidth=2, label=f’μ={mu}’) axes[idx].set_title(f’n={n}’, fontsize=12, fontweight=‘bold’) axes[idx].set_xlabel(‘Sample mean X̄’, fontsize=11) axes[idx].set_ylabel(‘Density’, fontsize=11) axes[idx].legend(fontsize=9) axes[idx].grid(alpha=0.3) # Statistics print(f”n={n:3d}: theoretical σ²/n={sigma2/n:6.2f}, ” f”empirical Var(X̄)={np.var(sample_means):6.2f}”) plt.suptitle(‘Distribution of Sample Mean’, fontsize=14, fontweight=‘bold’) plt.tight_layout() plt.show()

💻 Code Example 6: Distribution of Sample Variance

Distribution of sample variance np.random.seed(42) mu = 0 sigma = 5 n = 10 # Sample size n_trials = 10000 # Calculate sample variances sample_variances = [] for _ in range(n_trials): samples = np.random.normal(mu, sigma, n) sample_var = np.var(samples, ddof=1) # Unbiased variance (divide by n-1) sample_variances.append(sample_var) sample_variances = np.array(sample_variances) # Distribution of (n-1)S²/σ² (chi-square distribution) chi_squared_stat = (n - 1) * sample_variances / sigma2 fig, axes = plt.subplots(1, 2, figsize=(14, 5)) # (1) Distribution of sample variance axes[0].hist(sample_variances, bins=50, density=True, alpha=0.6, color=‘#667eea’, edgecolor=‘black’, linewidth=0.5, label=‘Simulation’) # Theoretical distribution (scaled chi-square distribution) x_var = np.linspace(0, 80, 1000) theoretical_var_pdf = stats.chi2.pdf((n-1)*x_var/sigma2, n-1) * (n-1)/sigma2 axes[0].plot(x_var, theoretical_var_pdf, ‘r-’, linewidth=2.5, label=‘Theoretical distribution’) axes[0].axvline(sigma2, color=‘green’, linestyle=’—’, linewidth=2, label=f’σ²={sigma2}’) axes[0].set_xlabel(‘Sample variance S²’, fontsize=11) axes[0].set_ylabel(‘Density’, fontsize=11) axes[0].set_title(‘Distribution of Sample Variance’, fontsize=12, fontweight=‘bold’) axes[0].legend() axes[0].grid(alpha=0.3) # (2) Comparison with chi-square distribution axes[1].hist(chi_squared_stat, bins=50, density=True, alpha=0.6, color=‘#764ba2’, edgecolor=‘black’, linewidth=0.5, label=‘Simulation’) x_chi = np.linspace(0, 30, 1000) theoretical_chi_pdf = stats.chi2.pdf(x_chi, n-1) axes[1].plot(x_chi, theoretical_chi_pdf, ‘r-’, linewidth=2.5, label=f’χ²({n-1})’) axes[1].set_xlabel(‘(n-1)S²/σ²’, fontsize=11) axes[1].set_ylabel(‘Density’, fontsize=11) axes[1].set_title(f’Comparison with Chi-Square Distribution χ²({n-1})’, fontsize=12, fontweight=‘bold’) axes[1].legend() axes[1].grid(alpha=0.3) plt.tight_layout() plt.show() print(“Verification of sample variance:”) print(f”Theoretical: E[S²] = σ² = {sigma2}”) print(f”Empirical: E[S²] = {np.mean(sample_variances):.4f}”) print(f”\nTheoretical: E[(n-1)S²/σ²] = n-1 = {n-1}”) print(f”Empirical: E[(n-1)S²/σ²] = {np.mean(chi_squared_stat):.4f}“)

2.6 Application to Materials Science Data

💻 Code Example 7: Application to Materials Science Data

Real data simulation: Material strength measurement np.random.seed(42) # Virtual data: Flexural strength of ceramics (MPa) # True distribution is Weibull (typical distribution for material strength) shape_param = 10 # Shape parameter scale_param = 500 # Characteristic strength (MPa) true_dist = stats.weibull_min(c=shape_param, scale=scale_param) true_mean = true_dist.mean() true_std = true_dist.std() print(“Statistical Analysis of Material Strength Data”) print(”=“*60) print(f”True distribution: Weibull distribution (shape={shape_param}, scale={scale_param})”) print(f”True mean: {true_mean:.2f} MPa”) print(f”True standard deviation: {true_std:.2f} MPa”) # Experiments: Measurements with different sample sizes sample_sizes = [3, 5, 10, 30, 50, 100] n_experiments = 1000 # Number of experiments fig, axes = plt.subplots(2, 3, figsize=(15, 10)) axes = axes.flatten() for idx, n in enumerate(sample_sizes): # Calculate sample means for each sample size sample_means = [] for _ in range(n_experiments): samples = true_dist.rvs(n) sample_means.append(samples.mean()) sample_means = np.array(sample_means) # Histogram axes[idx].hist(sample_means, bins=40, density=True, alpha=0.6, color=‘#667eea’, edgecolor=‘black’, linewidth=0.5, label=‘Measurement results’) # Approximation by central limit theorem (normal distribution) x = np.linspace(sample_means.min(), sample_means.max(), 1000) clt_pdf = stats.norm.pdf(x, true_mean, true_std/np.sqrt(n)) axes[idx].plot(x, clt_pdf, ‘r-’, linewidth=2.5, label=‘CLT approximation’) axes[idx].axvline(true_mean, color=‘green’, linestyle=’—’, linewidth=2, label=f’True mean’) axes[idx].set_title(f’Sample size n={n}’, fontsize=12, fontweight=‘bold’) axes[idx].set_xlabel(‘Average strength (MPa)’, fontsize=10) axes[idx].set_ylabel(‘Density’, fontsize=10) axes[idx].legend(fontsize=8) axes[idx].grid(alpha=0.3) # Statistical inference sem = true_std / np.sqrt(n) # Standard error ci_95 = stats.norm.interval(0.95, loc=true_mean, scale=sem) print(f”\nn={n:3d}:”) print(f” Standard error (SEM): {sem:.2f} MPa”) print(f” 95% confidence interval: [{ci_95[0]:.2f}, {ci_95[1]:.2f}] MPa”) print(f” Empirical standard deviation of mean: {np.std(sample_means):.2f} MPa”) plt.suptitle(‘Application of Central Limit Theorem in Material Strength Measurement’, fontsize=14, fontweight=‘bold’) plt.tight_layout() plt.show() # Practical inference print(“\n\n[Practical Inference]”) print(”=“*60) n_test = 10 # Actual number of test specimens test_samples = true_dist.rvs(n_test) test_mean = test_samples.mean() test_std = test_samples.std(ddof=1) print(f”Number of test specimens: {n_test}”) print(f”Measured average strength: {test_mean:.2f} MPa”) print(f”Measured standard deviation: {test_std:.2f} MPa”) # Confidence interval using t-distribution (true variance unknown) t_critical = stats.t.ppf(0.975, n_test - 1) ci_lower = test_mean - t_critical * test_std / np.sqrt(n_test) ci_upper = test_mean + t_critical * test_std / np.sqrt(n_test) print(f”\n95% confidence interval (t-distribution): [{ci_lower:.2f}, {ci_upper:.2f}] MPa”) print(f”True mean {true_mean:.2f} MPa is {‘included in’ if ci_lower <= true_mean <= ci_upper else ‘not included in’} the confidence interval”)

💡 Note: In materials science, it is necessary to estimate population characteristics from a limited number of test specimens. The central limit theorem enables t-tests and construction of confidence intervals even when the original distribution (such as Weibull distribution) is not normal, because the sample mean follows a normal distribution.

Exercises

📝 Exercise 1: Verification of Law of Large Numbers
Simulate a dice rolling experiment and verify the following:

  1. Visualize that the sample mean of the outcomes converges to the theoretical value of 3.5
  2. Calculate the upper bound of the probability that |X̄ₙ - 3.5| > 0.2 using Chebyshev’s inequality
  3. Compare the convergence rates for different sample sizes (n=10, 100, 1000)

📝 Exercise 2: Application of Central Limit Theorem
Take n samples independently from a uniform distribution U(0, 10) and calculate the sample mean:

  1. Visualize the distribution of sample means with histograms for n=5, 10, 30, 100
  2. Verify that the standardized sample mean follows a standard normal distribution using Q-Q plots
  3. Statistically verify normality using the Kolmogorov-Smirnov test

📝 Exercise 3: Construction of Confidence Intervals
Take a sample of n=25 from a normal distribution N(100, 15²) and calculate the following:

  1. Construct a 95% confidence interval for the sample mean (when population variance is known)
  2. Calculate the sample mean and sample variance, and construct a 95% confidence interval using the t-distribution (when population variance is unknown)
  3. Perform 1000 independent simulations and verify the proportion of times the confidence interval contains the true mean

← Chapter 1 Chapter 3 →

Disclaimer

← Previous Chapter Back to Series Index Next Chapter →