2.1 大数の弱法則
📊 定理: 大数の弱法則（Weak Law of Large Numbers）
独立同分布（i.i.d.）の確率変数 \(X_1, X_2, \ldots, X_n\) について、期待値 \(E[X_i] = \mu\)、分散 \(Var(X_i) = \sigma^2 < \infty\) とするとき、 標本平均 \(\bar{X}_n = \frac{1}{n}\sum_{i=1}^n X_i\) は \(\mu\) に確率収束します： \[\lim_{n \to \infty} P(|\bar{X}_n - \mu| > \epsilon) = 0 \quad \text{for any } \epsilon > 0\] 証明（チェビシェフの不等式による）: \[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)\]
独立同分布（i.i.d.）の確率変数 \(X_1, X_2, \ldots, X_n\) について、期待値 \(E[X_i] = \mu\)、分散 \(Var(X_i) = \sigma^2 < \infty\) とするとき、 標本平均 \(\bar{X}_n = \frac{1}{n}\sum_{i=1}^n X_i\) は \(\mu\) に確率収束します： \[\lim_{n \to \infty} P(|\bar{X}_n - \mu| > \epsilon) = 0 \quad \text{for any } \epsilon > 0\] 証明（チェビシェフの不等式による）: \[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)\]
💻 コード例1: 大数の弱法則のシミュレーション
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
# パラメータ設定
np.random.seed(42)
mu = 5 # 真の期待値
sigma = 2 # 真の標準偏差
max_n = 10000
epsilon = 0.5 # 許容誤差
# サンプリング（正規分布から）
samples = np.random.normal(mu, sigma, max_n)
# 累積平均の計算
cumulative_mean = np.cumsum(samples) / np.arange(1, max_n + 1)
# |X̄_n - μ| > ε となる確率の推定
n_values = np.logspace(1, 4, 100).astype(int)
n_values = np.unique(n_values) # 重複削除
prob_exceed = []
for n in n_values:
# 複数回のシミュレーション
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_bound = sigma**2 / (n_values * epsilon**2)
# 可視化
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (1) 標本平均の収束
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'真の期待値 μ={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('サンプルサイズ n', fontsize=11)
axes[0, 0].set_ylabel('標本平均 X̄ₙ', fontsize=11)
axes[0, 0].set_title('標本平均の収束', fontsize=12, fontweight='bold')
axes[0, 0].legend()
axes[0, 0].grid(alpha=0.3)
axes[0, 0].set_xlim([0, max_n])
# (2) 対数スケールでの収束
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('サンプルサイズ n (対数)', fontsize=11)
axes[0, 1].set_ylabel('標本平均 X̄ₙ', fontsize=11)
axes[0, 1].set_title('標本平均の収束（対数スケール）', fontsize=12, fontweight='bold')
axes[0, 1].grid(alpha=0.3, which='both')
# (3) P(|X̄ₙ - μ| > ε) の推定
axes[1, 0].loglog(n_values, prob_exceed, 'o-', color='#667eea',
linewidth=2, markersize=4, label='シミュレーション')
axes[1, 0].loglog(n_values, theoretical_bound, '--', color='red',
linewidth=2, label='チェビシェフ上界')
axes[1, 0].set_xlabel('サンプルサイズ n', fontsize=11)
axes[1, 0].set_ylabel('P(|X̄ₙ - μ| > ε)', fontsize=11)
axes[1, 0].set_title('大数の弱法則の検証', fontsize=12, fontweight='bold')
axes[1, 0].legend()
axes[1, 0].grid(alpha=0.3, which='both')
# (4) 異なる分布での比較
distributions = {
'正規分布': lambda n: np.random.normal(mu, sigma, n),
'一様分布': lambda n: np.random.uniform(mu-sigma*np.sqrt(3), mu+sigma*np.sqrt(3), n),
'指数分布': 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('サンプルサイズ n', fontsize=11)
axes[1, 1].set_ylabel('標本平均 X̄ₙ', fontsize=11)
axes[1, 1].set_title('異なる分布での大数の法則', 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("大数の弱法則の検証:")
print(f"真の期待値: μ = {mu}")
print(f"n=100での標本平均: {cumulative_mean[99]:.4f}")
print(f"n=1000での標本平均: {cumulative_mean[999]:.4f}")
print(f"n=10000での標本平均: {cumulative_mean[9999]:.4f}")
大数の弱法則の検証:
真の期待値: μ = 5
n=100での標本平均: 5.0234
n=1000での標本平均: 4.9876
n=10000での標本平均: 5.0012
2.2 大数の強法則
📊 定理: 大数の強法則（Strong Law of Large Numbers）
独立同分布の確率変数 \(X_1, X_2, \ldots\) について、\(E[|X_i|] < \infty\) ならば、 標本平均は期待値に概収束（almost surely convergence）します： \[P\left(\lim_{n \to \infty} \bar{X}_n = \mu\right) = 1\] 強法則 vs 弱法則:
独立同分布の確率変数 \(X_1, X_2, \ldots\) について、\(E[|X_i|] < \infty\) ならば、 標本平均は期待値に概収束（almost surely convergence）します： \[P\left(\lim_{n \to \infty} \bar{X}_n = \mu\right) = 1\] 強法則 vs 弱法則:
- 弱法則: 任意の \(\epsilon\) に対して \(P(|\bar{X}_n - \mu| > \epsilon) \to 0\)（確率収束）
- 強法則: \(P(\bar{X}_n \to \mu) = 1\)（概収束）
💻 コード例2: 大数の強法則の実証
# 複数の独立な経路でのシミュレーション
np.random.seed(123)
n_paths = 50 # 経路数
max_n = 5000
mu = 10
sigma = 3
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (1) 複数経路での標本平均の推移
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('サンプルサイズ n', fontsize=11)
axes[0, 0].set_ylabel('標本平均 X̄ₙ', fontsize=11)
axes[0, 0].set_title(f'{n_paths}個の独立経路（大数の強法則）', 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) 最大偏差の減衰
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='最大偏差')
axes[0, 1].loglog(n_checkpoints, sigma / np.sqrt(n_checkpoints), '--',
color='red', linewidth=2, label='σ/√n（理論）')
axes[0, 1].set_xlabel('サンプルサイズ n', fontsize=11)
axes[0, 1].set_ylabel('max |X̄ₙ - μ|', fontsize=11)
axes[0, 1].set_title('最大偏差の減衰', 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%信頼帯')
axes[1, 0].plot(median, color='#667eea', linewidth=2, label='中央値')
axes[1, 0].axhline(y=mu, color='red', linestyle='--', linewidth=2, label=f'μ={mu}')
axes[1, 0].set_xlabel('サンプルサイズ n', fontsize=11)
axes[1, 0].set_ylabel('標本平均 X̄ₙ', fontsize=11)
axes[1, 0].set_title('収束の信頼帯（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) 分散の減衰
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='実測分散')
axes[1, 1].loglog(x_axis, theoretical_var, '--', color='red',
linewidth=2, label='σ²/n（理論）')
axes[1, 1].set_xlabel('サンプルサイズ n', fontsize=11)
axes[1, 1].set_ylabel('Var(X̄ₙ)', fontsize=11)
axes[1, 1].set_title('標本平均の分散の減衰', fontsize=12, fontweight='bold')
axes[1, 1].legend()
axes[1, 1].grid(alpha=0.3, which='both')
plt.tight_layout()
plt.show()
print("大数の強法則の実証:")
print(f"500回の独立試行における、n=5000での標本平均の範囲:")
print(f" 最小値: {all_paths[:, -1].min():.4f}")
print(f" 最大値: {all_paths[:, -1].max():.4f}")
print(f" 中央値: {np.median(all_paths[:, -1]):.4f}")
print(f" 標準偏差: {np.std(all_paths[:, -1]):.4f}")
2.3 中心極限定理
📊 定理: 中心極限定理（Central Limit Theorem, CLT）
独立同分布の確率変数 \(X_1, X_2, \ldots, X_n\) について、期待値 \(E[X_i] = \mu\)、分散 \(Var(X_i) = \sigma^2 < \infty\) とするとき、 標準化した標本平均は標準正規分布に分布収束します： \[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)\] すなわち： \[\lim_{n \to \infty} P(Z_n \leq z) = \Phi(z) = \int_{-\infty}^z \frac{1}{\sqrt{2\pi}} e^{-t^2/2} dt\] 重要性: 元の分布に関わらず、標本平均は正規分布に近づきます。
独立同分布の確率変数 \(X_1, X_2, \ldots, X_n\) について、期待値 \(E[X_i] = \mu\)、分散 \(Var(X_i) = \sigma^2 < \infty\) とするとき、 標準化した標本平均は標準正規分布に分布収束します： \[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)\] すなわち： \[\lim_{n \to \infty} P(Z_n \leq z) = \Phi(z) = \int_{-\infty}^z \frac{1}{\sqrt{2\pi}} e^{-t^2/2} dt\] 重要性: 元の分布に関わらず、標本平均は正規分布に近づきます。
💻 コード例3: 中心極限定理の可視化
# 様々な分布から標本平均の分布を調べる
np.random.seed(42)
# 元の分布
distributions = {
'一様分布 U(0,1)': {
'generator': lambda n: np.random.uniform(0, 1, n),
'mu': 0.5,
'sigma': 1/np.sqrt(12)
},
'ベルヌーイ (p=0.3)': {
'generator': lambda n: np.random.binomial(1, 0.3, n),
'mu': 0.3,
'sigma': np.sqrt(0.3 * 0.7)
},
'指数分布 Exp(1)': {
'generator': lambda n: np.random.exponential(1, n),
'mu': 1,
'sigma': 1
},
'ベータ分布 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):
# 標本平均を計算
sample_means = []
for _ in range(n_trials):
samples = generator(n)
sample_means.append(samples.mean())
sample_means = np.array(sample_means)
# 標準化
z_scores = (sample_means - mu) / (sigma / np.sqrt(n))
# ヒストグラム
axes[i, j].hist(z_scores, bins=50, density=True, alpha=0.6,
color='#667eea', edgecolor='black', linewidth=0.5)
# 標準正規分布を重ねる
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('中心極限定理の可視化：様々な分布での検証',
fontsize=14, fontweight='bold', y=0.995)
plt.tight_layout()
plt.show()
print("中心極限定理の検証:")
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'](n).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}: 平均={np.mean(z_scores):6.3f}, 分散={np.var(z_scores):6.3f}")
2.4 非正規分布から正規分布への収束
💻 コード例4: 非正規分布から正規分布への収束
# 極端に非正規な分布での中心極限定理
np.random.seed(42)
# 極端に歪んだ分布（対数正規分布）
mu_log = 0
sigma_log = 1
def lognormal_generator(n):
return np.random.lognormal(mu_log, sigma_log, n)
# 対数正規分布の理論統計量
true_mean = np.exp(mu_log + sigma_log**2 / 2)
true_var = (np.exp(sigma_log**2) - 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):
# 標本平均を計算
sample_means = []
for _ in range(n_trials):
samples = lognormal_generator(n)
sample_means.append(samples.mean())
sample_means = np.array(sample_means)
# 標準化
z_scores = (sample_means - true_mean) / (true_sigma / np.sqrt(n))
# Q-Qプロット用のデータ
theoretical_quantiles = stats.norm.ppf(np.linspace(0.01, 0.99, 100))
sample_quantiles = np.percentile(z_scores, np.linspace(1, 99, 100))
# ヒストグラム + Q-Qプロット
ax_main = axes[idx]
# ヒストグラム
ax_main.hist(z_scores, bins=40, density=True, alpha=0.5,
color='#667eea', edgecolor='black', linewidth=0.5)
# 理論的な正規分布
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('標準化された標本平均', fontsize=10)
ax_main.set_ylabel('密度', fontsize=10)
ax_main.legend(fontsize=9)
ax_main.grid(alpha=0.3)
plt.suptitle('対数正規分布（極端に歪んだ分布）での中心極限定理',
fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
# Kolmogorov-Smirnov検定で正規性を検定
print("対数正規分布での中心極限定理の検証:")
print("="*60)
print(f"元の分布の歪度: {stats.lognorm.stats(s=sigma_log, moments='s'):.4f}")
print(f"元の分布の尖度: {stats.lognorm.stats(s=sigma_log, moments='k'):.4f}")
print("\nKolmogorov-Smirnov検定（正規性の検定）:")
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統計量={ks_stat:.4f}, p値={p_value:.4f} "
f"{'（正規分布に従う）' if p_value > 0.05 else '（有意に異なる）'}")
2.5 標本平均と標本分散の分布
📊 定理: 標本統計量の分布
\(X_1, \ldots, X_n \sim N(\mu, \sigma^2)\) （i.i.d.）のとき： 標本平均: \[\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i \sim N\left(\mu, \frac{\sigma^2}{n}\right)\] 標本分散: \[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統計量: \[T = \frac{\bar{X} - \mu}{S / \sqrt{n}} \sim t_{n-1}\]
\(X_1, \ldots, X_n \sim N(\mu, \sigma^2)\) （i.i.d.）のとき： 標本平均: \[\bar{X} = \frac{1}{n}\sum_{i=1}^n X_i \sim N\left(\mu, \frac{\sigma^2}{n}\right)\] 標本分散: \[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統計量: \[T = \frac{\bar{X} - \mu}{S / \sqrt{n}} \sim t_{n-1}\]
💻 コード例5: 標本平均の分布
# 標本平均の分布
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):
# 標本平均を計算
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)
# ヒストグラム
axes[idx].hist(sample_means, bins=50, density=True, alpha=0.6,
color='#667eea', edgecolor='black', linewidth=0.5,
label='シミュレーション')
# 理論分布 N(μ, σ²/n)
x = np.linspace(mu - 4*sigma/np.sqrt(n), mu + 4*sigma/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}, {sigma**2/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('標本平均 X̄', fontsize=11)
axes[idx].set_ylabel('密度', fontsize=11)
axes[idx].legend(fontsize=9)
axes[idx].grid(alpha=0.3)
# 統計量
print(f"n={n:3d}: 理論 σ²/n={sigma**2/n:6.2f}, "
f"実測 Var(X̄)={np.var(sample_means):6.2f}")
plt.suptitle('標本平均の分布', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
💻 コード例6: 標本分散の分布
# 標本分散の分布
np.random.seed(42)
mu = 0
sigma = 5
n = 10 # サンプルサイズ
n_trials = 10000
# 標本分散を計算
sample_variances = []
for _ in range(n_trials):
samples = np.random.normal(mu, sigma, n)
sample_var = np.var(samples, ddof=1) # 不偏分散（n-1で割る）
sample_variances.append(sample_var)
sample_variances = np.array(sample_variances)
# (n-1)S²/σ² の分布（カイ二乗分布）
chi_squared_stat = (n - 1) * sample_variances / sigma**2
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# (1) 標本分散の分布
axes[0].hist(sample_variances, bins=50, density=True, alpha=0.6,
color='#667eea', edgecolor='black', linewidth=0.5,
label='シミュレーション')
# 理論分布（スケールされたカイ二乗分布）
x_var = np.linspace(0, 80, 1000)
theoretical_var_pdf = stats.chi2.pdf((n-1)*x_var/sigma**2, n-1) * (n-1)/sigma**2
axes[0].plot(x_var, theoretical_var_pdf, 'r-', linewidth=2.5,
label='理論分布')
axes[0].axvline(sigma**2, color='green', linestyle='--', linewidth=2,
label=f'σ²={sigma**2}')
axes[0].set_xlabel('標本分散 S²', fontsize=11)
axes[0].set_ylabel('密度', fontsize=11)
axes[0].set_title('標本分散の分布', fontsize=12, fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)
# (2) カイ二乗分布との比較
axes[1].hist(chi_squared_stat, bins=50, density=True, alpha=0.6,
color='#764ba2', edgecolor='black', linewidth=0.5,
label='シミュレーション')
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('密度', fontsize=11)
axes[1].set_title(f'カイ二乗分布 χ²({n-1}) との比較', fontsize=12, fontweight='bold')
axes[1].legend()
axes[1].grid(alpha=0.3)
plt.tight_layout()
plt.show()
print("標本分散の検証:")
print(f"理論: E[S²] = σ² = {sigma**2}")
print(f"実測: E[S²] = {np.mean(sample_variances):.4f}")
print(f"\n理論: E[(n-1)S²/σ²] = n-1 = {n-1}")
print(f"実測: E[(n-1)S²/σ²] = {np.mean(chi_squared_stat):.4f}")
2.6 材料科学データへの応用
💻 コード例7: 材料科学データへの応用
# 実データシミュレーション：材料強度の測定
np.random.seed(42)
# 仮想データ：セラミックスの曲げ強度（MPa）
# 真の分布はワイブル分布（材料強度の典型的な分布）
shape_param = 10 # 形状パラメータ
scale_param = 500 # 特性強度 (MPa)
true_dist = stats.weibull_min(c=shape_param, scale=scale_param)
true_mean = true_dist.mean()
true_std = true_dist.std()
print("材料強度データの統計解析")
print("="*60)
print(f"真の分布: ワイブル分布 (shape={shape_param}, scale={scale_param})")
print(f"真の平均: {true_mean:.2f} MPa")
print(f"真の標準偏差: {true_std:.2f} MPa")
# 実験: 異なるサンプルサイズでの測定
sample_sizes = [3, 5, 10, 30, 50, 100]
n_experiments = 1000 # 実験回数
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()
for idx, n in enumerate(sample_sizes):
# 各サンプルサイズでの標本平均を計算
sample_means = []
for _ in range(n_experiments):
samples = true_dist.rvs(n)
sample_means.append(samples.mean())
sample_means = np.array(sample_means)
# ヒストグラム
axes[idx].hist(sample_means, bins=40, density=True, alpha=0.6,
color='#667eea', edgecolor='black', linewidth=0.5,
label='測定結果')
# 中心極限定理による近似（正規分布）
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による近似')
axes[idx].axvline(true_mean, color='green', linestyle='--',
linewidth=2, label=f'真の平均')
axes[idx].set_title(f'試料数 n={n}', fontsize=12, fontweight='bold')
axes[idx].set_xlabel('平均強度 (MPa)', fontsize=10)
axes[idx].set_ylabel('密度', fontsize=10)
axes[idx].legend(fontsize=8)
axes[idx].grid(alpha=0.3)
# 統計的推論
sem = true_std / np.sqrt(n) # 標準誤差
ci_95 = stats.norm.interval(0.95, loc=true_mean, scale=sem)
print(f"\nn={n:3d}:")
print(f" 標準誤差 (SEM): {sem:.2f} MPa")
print(f" 95%信頼区間: [{ci_95[0]:.2f}, {ci_95[1]:.2f}] MPa")
print(f" 実測平均の標準偏差: {np.std(sample_means):.2f} MPa")
plt.suptitle('材料強度測定における中心極限定理の応用',
fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()
# 実用的な推論
print("\n\n【実用的な推論】")
print("="*60)
n_test = 10 # 実際の試験片数
test_samples = true_dist.rvs(n_test)
test_mean = test_samples.mean()
test_std = test_samples.std(ddof=1)
print(f"試験片数: {n_test}")
print(f"測定された平均強度: {test_mean:.2f} MPa")
print(f"測定された標準偏差: {test_std:.2f} MPa")
# t分布による信頼区間（真の分散が未知）
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%信頼区間（t分布）: [{ci_lower:.2f}, {ci_upper:.2f}] MPa")
print(f"真の平均 {true_mean:.2f} MPa は信頼区間に{'含まれる' if ci_lower <= true_mean <= ci_upper else '含まれない'}")
💡 Note: 材料科学では、限られた試験片数から母集団の特性を推定する必要があります。中心極限定理により、元の分布（ワイブル分布など）が正規分布でなくても、標本平均は正規分布に従うため、t検定や信頼区間の構築が可能になります。
演習問題
📝 演習1: 大数の法則の検証
サイコロを振る実験をシミュレートし、以下を確認せよ：
サイコロを振る実験をシミュレートし、以下を確認せよ：
- 出目の標本平均が理論値3.5に収束することを可視化せよ
- チェビシェフの不等式を用いて、|X̄ₙ - 3.5| > 0.2 となる確率の上界を計算せよ
- 異なるサンプルサイズ（n=10, 100, 1000）での収束速度を比較せよ
📝 演習2: 中心極限定理の応用
一様分布 U(0, 10) から独立に n 個のサンプルを取り、標本平均を計算する：
一様分布 U(0, 10) から独立に n 個のサンプルを取り、標本平均を計算する：
- n=5, 10, 30, 100 の各場合で、標本平均の分布をヒストグラムで可視化せよ
- 標準化した標本平均が標準正規分布に従うことを Q-Q プロットで確認せよ
- Kolmogorov-Smirnov検定を用いて、正規性を統計的に検証せよ
📝 演習3: 信頼区間の構築
正規分布 N(100, 15²) から n=25 のサンプルを取り、以下を計算せよ：
正規分布 N(100, 15²) から n=25 のサンプルを取り、以下を計算せよ：
- 標本平均の95%信頼区間を構築せよ（母分散既知の場合）
- 標本平均と標本分散を計算し、t分布を用いた95%信頼区間を構築せよ（母分散未知の場合）
- 1000回の独立なシミュレーションを行い、信頼区間が真の平均を含む割合を確認せよ