第3章：統計的推定と仮説検定

イントロダクション

第1章では、手元にあるデータの特徴を要約する記述統計と、不確実性を数学的に扱う確率論の基礎を学びました。しかし、私たちが実際に手にできるデータは、興味の対象となる集団全体（母集団、英語では Population）のごく一部、すなわち標本（Sample）にすぎません。工場の全製品を検査することはできず、世界中の患者全員に新薬を試すこともできません。そこで必要になるのが、限られた標本から母集団の性質を推し量る統計的推論（Statistical Inference）です。

本章では、統計的推論の二本柱である推定（Estimation）と仮説検定（Hypothesis Testing）を扱います。推定では、母集団のパラメータ（平均や分散など）を標本からどのように見積もるかを学びます。仮説検定では、「新しい学習法は本当に効果があるのか」「2つの製造ラインの品質に差はあるのか」といった問いに、データに基づいて客観的な結論を下す方法を学びます。これらは機械学習においても、モデルの性能比較、A/Bテスト、特徴量の有意性評価など、あらゆる場面で使われる基盤技術です。

1. 点推定と区間推定の理論

母集団のパラメータ（母数）を1つの値で推測するアプローチと、値の存在しうる範囲で推測するアプローチの2種類があります。

1.1 点推定（Point Estimation）

点推定とは、標本データから母集団のパラメータをただ1つの数値で推測することです。この推測に使う計算方法や統計量を推定量（Estimator）と呼び、実際に計算された数値を推定値（Estimate）と呼びます。例えば、標本平均 $\bar{x}$ は母平均 $\mu$ の推定量であり、実際に計算した「82.5点」のような数値が推定値です。

推定量に求められる性質

良い推定量には、次のような統計的性質が求められます。

• 不偏性（Unbiasedness）: 推定量の期待値が真のパラメータと一致する性質です。数式では $E[\hat{\theta}] = \theta$ と表されます。推定量が不偏推定量（Unbiased Estimator）であるとき、同じ手続きで標本抽出を何度も繰り返せば、推定値の平均は真の値に収束します。
• 一致性（Consistency）: 標本サイズ $n$ を大きくするほど、推定量が真のパラメータに確率的に近づいていく性質です。
• 有効性（Efficiency）: 不偏推定量の中で、分散が最も小さい（推定値のばらつきが最も少ない）性質です。

1.2 区間推定（Interval Estimation）

点推定は「1つの値」で母数を言い当てようとしますが、標本ごとに推定値は変動するため、その値がどれだけ信頼できるかという情報が欠けています。そこで区間推定では、母数が含まれると考えられる範囲を、一定の信頼度とともに示します。区間推定の代表例が、次の第3節で詳しく扱う信頼区間（Confidence Interval）です。区間推定は「この範囲のどこかに真の値がありそうだ」という、点推定よりも情報量の多い推測を与えてくれます。

1.3 Pythonでの実装：点推定と不偏性の確認

成人男性の身長データを例に、点推定を実装し、不偏分散が真の母分散に近づくことをシミュレーションで確認しましょう。

import numpy as np

np.random.seed(42)

# 母集団: 平均170cm、標準偏差8cmの正規分布に従う成人男性の身長
population_mean = 170
population_std = 8

# 標本を1つ抽出して点推定
sample = np.random.normal(population_mean, population_std, size=30)
point_estimate_mean = np.mean(sample)
point_estimate_var_biased = np.var(sample, ddof=0)      # 標本分散（バイアスあり）
point_estimate_var_unbiased = np.var(sample, ddof=1)     # 不偏分散

print(f"標本平均（母平均の点推定値）: {point_estimate_mean:.3f} cm")
print(f"標本分散（ddof=0、バイアスあり）: {point_estimate_var_biased:.3f}")
print(f"不偏分散（ddof=1）: {point_estimate_var_unbiased:.3f}")
print(f"真の母平均: {population_mean} cm（通常は未知）")

print()
print("=== 推定量の不偏性をシミュレーションで確認 ===")
n_simulations = 10000
sample_size = 30

biased_var_estimates = np.zeros(n_simulations)
unbiased_var_estimates = np.zeros(n_simulations)
mean_estimates = np.zeros(n_simulations)

for i in range(n_simulations):
    s = np.random.normal(population_mean, population_std, size=sample_size)
    mean_estimates[i] = np.mean(s)
    biased_var_estimates[i] = np.var(s, ddof=0)
    unbiased_var_estimates[i] = np.var(s, ddof=1)

true_variance = population_std ** 2
print(f"真の母分散: {true_variance}")
print(f"標本平均の平均（{n_simulations}回試行）: {np.mean(mean_estimates):.4f}（真値: {population_mean}）")
print(f"標本分散(ddof=0)の平均: {np.mean(biased_var_estimates):.4f} → 真値からの乖離: {np.mean(biased_var_estimates) - true_variance:.4f}")
print(f"不偏分散(ddof=1)の平均: {np.mean(unbiased_var_estimates):.4f} → 真値からの乖離: {np.mean(unbiased_var_estimates) - true_variance:.4f}")

実行結果:

標本平均（母平均の点推定値）: 168.495 cm
標本分散（ddof=0、バイアスあり）: 50.113
不偏分散（ddof=1）: 51.841
真の母平均: 170 cm（通常は未知）

=== 推定量の不偏性をシミュレーションで確認 ===
真の母分散: 64
標本平均の平均（10000回試行）: 170.0013（真値: 170）
標本分散(ddof=0)の平均: 61.8140 → 真値からの乖離: -2.1860
不偏分散(ddof=1)の平均: 63.9455 → 真値からの乖離: -0.0545

2. 最尤推定法の原理と実装

最尤推定法（Maximum Likelihood Estimation, MLE）は、点推定を行うための最も重要な方法の1つです。「観測されたデータが得られる確率（尤もらしさ）を最大にするパラメータを、真のパラメータの推定値とする」という考え方に基づいています。

2.1 尤度関数と対数尤度

あるパラメータ $\theta$ を持つ確率分布から、独立に $n$ 個のデータ $x_1, x_2, \ldots, x_n$ が観測されたとします。この観測結果が得られる確率（同時確率密度）を $\theta$ の関数とみなしたものを尤度関数（Likelihood Function）と呼びます。

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

ここで $f(x; \theta)$ はパラメータ $\theta$ を持つ確率密度関数です。積のままでは数値計算上扱いにくい（非常に小さな値になり、アンダーフローを起こしやすい）ため、対数を取った対数尤度（Log-Likelihood）を最大化するのが一般的です。対数は単調増加関数なので、対数尤度を最大化するパラメータは尤度関数を最大化するパラメータと一致します。

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

2.2 最尤推定量の求め方

最尤推定量（Maximum Likelihood Estimator） $\hat{\theta}_{MLE}$ は、対数尤度 $\ell(\theta)$ を最大にする $\theta$ として定義されます。

$$\hat{\theta}_{MLE} = \underset{\theta}{\arg\max}\ \ell(\theta)$$

正規分布の場合、対数尤度を平均 $\mu$ と分散 $\sigma^2$ でそれぞれ微分して0と置く（尤度方程式を解く）ことで、次の閉形式解が得られます。

$$\hat{\mu}_{MLE} = \bar{x} = \frac{1}{n}\sum_{i=1}^{n}x_i \qquad \hat{\sigma}^2_{MLE} = \frac{1}{n}\sum_{i=1}^{n}(x_i - \bar{x})^2$$

2.3 Pythonでの実装：数値最適化による最尤推定

正規分布と指数分布を例に、対数尤度を数値的に最大化してパラメータを推定し、閉形式解と一致することを確認します。

import numpy as np
from scipy import stats
from scipy.optimize import minimize

np.random.seed(42)

# --- 例1: 正規分布の平均・標準偏差を最尤推定 ---
true_mu, true_sigma = 170, 8
data = np.random.normal(true_mu, true_sigma, size=50)

def neg_log_likelihood_normal(params, data):
    mu, sigma = params
    if sigma <= 0:
        return np.inf
    return -np.sum(stats.norm.logpdf(data, loc=mu, scale=sigma))

initial_guess = [np.mean(data), np.std(data)]
result = minimize(neg_log_likelihood_normal, initial_guess, args=(data,), method='Nelder-Mead')
mle_mu, mle_sigma = result.x

# 閉形式解（正規分布の場合、MLEは標本平均と標本標準偏差(ddof=0)に一致）
closed_form_mu = np.mean(data)
closed_form_sigma = np.std(data, ddof=0)

# scipy.stats.norm.fit も内部的に同じ最尤推定を行う
fit_mu, fit_sigma = stats.norm.fit(data)

print("=== 正規分布パラメータの最尤推定 ===")
print(f"数値最適化によるMLE: mu={mle_mu:.4f}, sigma={mle_sigma:.4f}")
print(f"閉形式解:           mu={closed_form_mu:.4f}, sigma={closed_form_sigma:.4f}")
print(f"scipy.stats.norm.fit: mu={fit_mu:.4f}, sigma={fit_sigma:.4f}")

# --- 例2: 指数分布のレートパラメータを最尤推定 ---
print()
print("=== 指数分布パラメータの最尤推定 ===")
true_lambda = 0.5  # 平均待ち時間 = 1/lambda = 2分
wait_times = np.random.exponential(scale=1/true_lambda, size=200)

def neg_log_likelihood_exp(params, data):
    lam = params[0]
    if lam <= 0:
        return np.inf
    return -np.sum(stats.expon.logpdf(data, scale=1/lam))

result_exp = minimize(neg_log_likelihood_exp, [1.0], args=(wait_times,), method='Nelder-Mead')
mle_lambda = result_exp.x[0]

# 指数分布のMLEの閉形式解: lambda_hat = 1 / 標本平均
closed_form_lambda = 1 / np.mean(wait_times)

print(f"真のパラメータ lambda: {true_lambda}")
print(f"数値最適化によるMLE: lambda={mle_lambda:.4f}")
print(f"閉形式解 (1/標本平均): lambda={closed_form_lambda:.4f}")
print(f"標本サイズ: {len(wait_times)}, 標本平均待ち時間: {np.mean(wait_times):.4f}分")

実行結果:

=== 正規分布パラメータの最尤推定 ===
数値最適化によるMLE: mu=168.1962, sigma=7.3943
閉形式解:           mu=168.1962, sigma=7.3943
scipy.stats.norm.fit: mu=168.1962, sigma=7.3943

=== 指数分布パラメータの最尤推定 ===
真のパラメータ lambda: 0.5
数値最適化によるMLE: lambda=0.4982
閉形式解 (1/標本平均): lambda=0.4983
標本サイズ: 200, 標本平均待ち時間: 2.0069分

3通りの計算方法（数値最適化・閉形式解・scipyの組み込み関数）がすべて同じ結果を与えていることが確認できます。指数分布の例でも、数値最適化によるMLEと閉形式解 $\hat{\lambda} = 1/\bar{x}$ がほぼ一致しており、最尤推定法が理論通りに機能していることが分かります。

3. 信頼区間の計算と解釈

3.1 信頼区間とは

信頼区間（Confidence Interval, CI）は、区間推定の代表的な手法で、母数が含まれると期待される範囲を、指定した信頼水準（Confidence Level）（多くの場合95%）とともに示します。母平均 $\mu$ の信頼区間は、標本平均 $\bar{x}$ を中心に、標本のばらつきを表す標準誤差（Standard Error） $SE = s/\sqrt{n}$ を使って構成されます。

母標準偏差が未知で標本サイズが小さい場合（実務ではこちらが大半）、標準正規分布の代わりにt分布（t-Distribution）を用います。

$$\bar{x} \pm t_{\alpha/2, n-1} \times \frac{s}{\sqrt{n}}$$

ここで $t_{\alpha/2, n-1}$ は自由度 $n-1$ のt分布における臨界値です。標本サイズが大きくなるほどt分布は標準正規分布に近づきます。

3.2 Pythonでの実装：信頼区間の計算と被覆率の検証

身長データを例に95%信頼区間を計算し、さらに1,000回のシミュレーションによって「95%」という数字の意味を実際に検証してみましょう。

import numpy as np
from scipy import stats

np.random.seed(42)

# --- 母平均の信頼区間（母標準偏差未知、t分布を使用） ---
population_mean = 170
population_std = 8
sample = np.random.normal(population_mean, population_std, size=30)

n = len(sample)
sample_mean = np.mean(sample)
sample_std = np.std(sample, ddof=1)  # 不偏標準偏差
standard_error = sample_std / np.sqrt(n)

confidence_level = 0.95
degrees_of_freedom = n - 1

# t分布を用いた信頼区間（scipy.stats.t.interval）
ci_lower, ci_upper = stats.t.interval(confidence_level, degrees_of_freedom,
                                       loc=sample_mean, scale=standard_error)

print("=== 母平均の95%信頼区間 ===")
print(f"標本サイズ: {n}")
print(f"標本平均: {sample_mean:.3f} cm")
print(f"標本標準偏差: {sample_std:.3f} cm")
print(f"標準誤差: {standard_error:.3f} cm")
print(f"自由度: {degrees_of_freedom}")
print(f"95%信頼区間: [{ci_lower:.3f}, {ci_upper:.3f}] cm")

# 手計算での検証
t_critical = stats.t.ppf(1 - (1 - confidence_level) / 2, degrees_of_freedom)
margin_of_error = t_critical * standard_error
print(f"\n手計算での検証:")
print(f"t臨界値 (alpha=0.05, df={degrees_of_freedom}): {t_critical:.4f}")
print(f"誤差の範囲（マージン）: {margin_of_error:.3f}")
print(f"信頼区間: [{sample_mean - margin_of_error:.3f}, {sample_mean + margin_of_error:.3f}]")

# --- 信頼区間の被覆確率をシミュレーションで検証 ---
print()
print("=== 信頼区間の「95%」の意味をシミュレーションで確認 ===")
n_simulations = 1000
contains_true_mean = 0

for _ in range(n_simulations):
    s = np.random.normal(population_mean, population_std, size=30)
    s_mean = np.mean(s)
    s_se = np.std(s, ddof=1) / np.sqrt(len(s))
    lo, hi = stats.t.interval(0.95, len(s) - 1, loc=s_mean, scale=s_se)
    if lo <= population_mean <= hi:
        contains_true_mean += 1

coverage_rate = contains_true_mean / n_simulations
print(f"{n_simulations}回中、真の母平均({population_mean})を含んだ信頼区間の数: {contains_true_mean}")
print(f"被覆率: {coverage_rate * 100:.2f}%（理論値: 95%）")

実行結果:

=== 母平均の95%信頼区間 ===
標本サイズ: 30
標本平均: 168.495 cm
標本標準偏差: 7.200 cm
標準誤差: 1.315 cm
自由度: 29
95%信頼区間: [165.806, 171.183] cm

手計算での検証:
t臨界値 (alpha=0.05, df=29): 2.0452
誤差の範囲（マージン）: 2.689
信頼区間: [165.806, 171.183]

=== 信頼区間の「95%」の意味をシミュレーションで確認 ===
1000回中、真の母平均(170)を含んだ信頼区間の数: 948
被覆率: 94.80%（理論値: 95%）

4. 仮説検定の枠組み（帰無仮説、対立仮説、p値）

仮説検定（Hypothesis Testing）は、標本データに基づいて、母集団に関する仮説の妥当性を統計的に判断する手続きです。

4.1 帰無仮説と対立仮説

• 帰無仮説（Null Hypothesis, $H_0$）: 「差がない」「効果がない」など、検証したい主張の反対にあたる、棄却されることを想定した仮説です。
• 対立仮説（Alternative Hypothesis, $H_1$ または $H_a$）: 「差がある」「効果がある」など、私たちが本当に主張したい仮説です。

仮説検定では、直接 $H_1$ が正しいことを証明するのではなく、「$H_0$ が正しいと仮定した場合に、観測されたデータ（あるいはそれ以上に極端なデータ）が得られる確率がどれくらい低いか」を評価し、その確率が十分に低ければ $H_0$ を棄却する、という背理法に似た論理構造を取ります。

4.2 有意水準とp値

• 有意水準（Significance Level, $\alpha$）: 帰無仮説を誤って棄却してしまうことを許容する確率の基準で、慣習的に0.05（5%）や0.01（1%）が使われます。
• 検定統計量（Test Statistic）: 標本データから計算される、仮説を判定するための統計量（t値やカイ二乗値など）です。
• p値（p-value）: 帰無仮説が正しいと仮定したときに、実際に観測された検定統計量と同じか、それよりも極端な値が得られる確率です。

判定のルールはシンプルです。p値が有意水準 $\alpha$ を下回れば、帰無仮説を棄却し、対立仮説を支持します。下回らなければ、「帰無仮説を棄却する根拠が不十分である」と判断します（これは「帰無仮説が正しいと証明された」ことを意味しない点に注意してください）。

4.3 第一種の過誤と第二種の過誤

第一種の過誤（Type I Error）は、本当は差がないのに「差がある」と誤って判定してしまう過誤（偽陽性）で、その確率は有意水準 $\alpha$ そのものです。第二種の過誤（Type II Error）は、本当は差があるのに「差がない」と誤って判定してしまう過誤（偽陰性）です。この2つの過誤はトレードオフの関係にあり、$\alpha$ を厳しくすると第二種の過誤が増えやすくなります。

4.4 両側検定と片側検定

両側検定（Two-tailed Test）は「差がある（大小は問わない）」ことを検証する検定で、対立仮説は $\mu \neq \mu_0$ の形になります。片側検定（One-tailed Test）は「効果がある方向」まで指定して検証する検定で、対立仮説は $\mu > \mu_0$ または $\mu < \mu_0$ の形になります。多くの分析では、恣意的な結論を避けるために両側検定がデフォルトとして推奨されます。

4.5 Pythonでの実装：仮説検定の枠組みを手動計算で理解する

新しい学習法を受けた生徒の平均点が、従来の平均点75点と異なるかどうかを、1標本t検定の手続きに沿って手動で計算してみましょう。

import numpy as np
from scipy import stats

np.random.seed(42)

# 帰無仮説 H0: 新しい学習法を受けた生徒の平均点は従来と同じ 75点
# 対立仮説 H1: 新しい学習法を受けた生徒の平均点は 75点と異なる（両側検定）
population_mean_h0 = 75
sample_scores = np.array([78, 82, 75, 88, 79, 85, 91, 76, 83, 80,
                           77, 86, 89, 74, 81, 84, 90, 78, 82, 87])

n = len(sample_scores)
sample_mean = np.mean(sample_scores)
sample_std = np.std(sample_scores, ddof=1)
standard_error = sample_std / np.sqrt(n)

# t統計量を手動で計算
t_statistic = (sample_mean - population_mean_h0) / standard_error
degrees_of_freedom = n - 1

# 両側検定のp値
p_value = 2 * (1 - stats.t.cdf(abs(t_statistic), degrees_of_freedom))

alpha = 0.05

print("=== 仮説検定の枠組み：1標本t検定 ===")
print(f"帰無仮説 H0: 母平均 = {population_mean_h0}")
print(f"対立仮説 H1: 母平均 != {population_mean_h0}（両側検定）")
print(f"有意水準 alpha: {alpha}")
print(f"標本サイズ: {n}, 標本平均: {sample_mean:.3f}, 標本標準偏差: {sample_std:.3f}")
print(f"検定統計量 t: {t_statistic:.4f}")
print(f"自由度: {degrees_of_freedom}")
print(f"p値: {p_value:.2e}")

if p_value < alpha:
    print(f"結論: p値が有意水準を下回るため、帰無仮説を棄却する")
else:
    print(f"結論: p値が有意水準を下回らないため、帰無仮説を棄却しない")

# scipy.stats.ttest_1sampとの比較検証
t_check, p_check = stats.ttest_1samp(sample_scores, population_mean_h0)
print(f"\nscipy.stats.ttest_1sampでの検証: t={t_check:.4f}, p={p_check:.2e}")

実行結果:

=== 仮説検定の枠組み：1標本t検定 ===
帰無仮説 H0: 母平均 = 75
対立仮説 H1: 母平均 != 75（両側検定）
有意水準 alpha: 0.05
標本サイズ: 20, 標本平均: 82.250, 標本標準偏差: 5.149
検定統計量 t: 6.2968
自由度: 19
p値: 4.81e-06
結論: p値が有意水準を下回るため、帰無仮説を棄却する

scipy.stats.ttest_1sampでの検証: t=6.2968, p=4.81e-06

手動計算とscipyの組み込み関数 ttest_1samp の結果が完全に一致しました。p値（約0.0000048）は有意水準0.05を大きく下回っているため、「新しい学習法の平均点は従来と変わらない」という帰無仮説を棄却し、有意な差があると結論づけられます。

5. t検定、カイ二乗検定、F検定の実践

ここからは、実務で頻繁に使われる3種類の検定手法を、目的別に整理しながら実装します。

5.1 t検定（t-test）

t検定は、平均値の比較にt分布を用いる検定です。データが独立した2群から得られたか、同じ対象を前後で測定した対応のあるデータかによって、使い分けが必要です。

import numpy as np
from scipy import stats

np.random.seed(42)

def judge(p, alpha=0.05):
    return "有意差あり" if p < alpha else "有意差なし"

print("=" * 50)
print("1) 対応のない2標本t検定（独立2群の平均比較）")
print("=" * 50)
# 2つの学習法（A法・B法）でのテスト得点を比較
group_a = np.array([72, 75, 78, 80, 74, 77, 81, 79, 76, 73])
group_b = np.array([80, 85, 83, 88, 82, 86, 90, 84, 87, 81])

t_stat, p_value = stats.ttest_ind(group_a, group_b)
print(f"グループA平均: {np.mean(group_a):.2f}, グループB平均: {np.mean(group_b):.2f}")
print(f"t統計量: {t_stat:.4f}")
print(f"p値: {p_value:.3e}")
print(f"判定 (alpha=0.05): {judge(p_value)}")

print()
print("=" * 50)
print("2) 対応のある2標本t検定（paired t-test）")
print("=" * 50)
# 研修の前後でのスキルテストの点数（同じ被験者）
before = np.array([65, 70, 68, 72, 75, 69, 71, 74, 66, 73])
after = np.array([72, 78, 74, 80, 82, 75, 79, 81, 73, 80])

t_stat_paired, p_value_paired = stats.ttest_rel(before, after)
print(f"研修前平均: {np.mean(before):.2f}, 研修後平均: {np.mean(after):.2f}")
print(f"平均差: {np.mean(after - before):.2f}")
print(f"t統計量: {t_stat_paired:.4f}")
print(f"p値: {p_value_paired:.3e}")
print(f"判定 (alpha=0.05): {judge(p_value_paired)}")

実行結果:

==================================================
1) 対応のない2標本t検定（独立2群の平均比較）
==================================================
グループA平均: 76.50, グループB平均: 84.60
t統計量: -5.8105
p値: 1.665e-05
判定 (alpha=0.05): 有意差あり

==================================================
2) 対応のある2標本t検定（paired t-test）
==================================================
研修前平均: 70.30, 研修後平均: 77.40
平均差: 7.10
t統計量: -30.4286
p値: 2.189e-10
判定 (alpha=0.05): 有意差あり

5.2 カイ二乗検定（Chi-square Test）

カイ二乗検定は、カテゴリカルデータ（度数データ）を対象とした検定です。観測度数と期待度数のズレを、カイ二乗統計量として評価します。代表的な用途は適合度検定（Goodness-of-fit Test）（観測データが特定の理論分布に従うか）と独立性検定（Test of Independence）（2つのカテゴリ変数に関連があるか）です。独立性検定では、行と列にカテゴリを並べた分割表（Contingency Table）を使います。

import numpy as np
from scipy import stats

def judge(p, alpha=0.05):
    return "有意差あり" if p < alpha else "有意差なし"

print("=" * 50)
print("3) カイ二乗適合度検定（Goodness-of-fit test）")
print("=" * 50)
# サイコロが公正（各目が均等な確率）かどうかを検定
observed = np.array([18, 22, 16, 24, 20, 20])  # 各目の観測度数（合計120回）
expected = np.array([20, 20, 20, 20, 20, 20])   # 公正なサイコロの期待度数

chi2_stat, p_value_chi2 = stats.chisquare(observed, expected)
print(f"観測度数: {observed}")
print(f"期待度数: {expected}")
print(f"カイ二乗統計量: {chi2_stat:.4f}")
print(f"p値: {p_value_chi2:.4f}")
print(f"判定 (alpha=0.05): {judge(p_value_chi2)}")

print()
print("=" * 50)
print("4) カイ二乗独立性検定（Test of independence）")
print("=" * 50)
# 性別と製品の好み（A/B/C）の間に関連があるかを検定
# 行: 性別（男性・女性）、列: 好みの製品（A・B・C）
contingency_table = np.array([
    [30, 20, 15],   # 男性
    [15, 25, 25],   # 女性
])
chi2_ind, p_ind, dof_ind, expected_ind = stats.chi2_contingency(contingency_table)
print(f"分割表:\n{contingency_table}")
print(f"期待度数:\n{np.round(expected_ind, 2)}")
print(f"カイ二乗統計量: {chi2_ind:.4f}")
print(f"自由度: {dof_ind}")
print(f"p値: {p_ind:.4f}")
print(f"判定 (alpha=0.05): {judge(p_ind)}")

実行結果:

==================================================
3) カイ二乗適合度検定（Goodness-of-fit test）
==================================================
観測度数: [18 22 16 24 20 20]
期待度数: [20 20 20 20 20 20]
カイ二乗統計量: 2.0000
p値: 0.8491
判定 (alpha=0.05): 有意差なし

==================================================
4) カイ二乗独立性検定（Test of independence）
==================================================
分割表:
[[30 20 15]
 [15 25 25]]
期待度数:
[[22.5 22.5 20. ]
 [22.5 22.5 20. ]]
カイ二乗統計量: 8.0556
自由度: 2
p値: 0.0178
判定 (alpha=0.05): 有意差あり

適合度検定の結果、p値0.8491は0.05を大きく上回っており、「このサイコロは公正である」という帰無仮説を棄却する根拠はありません。一方、独立性検定ではp値0.0178が0.05を下回っており、性別と製品の好みには統計的に有意な関連があると判断できます。

5.3 F検定（F-test）

F検定は、2つの分散の比（F統計量）がF分布に従うことを利用した検定です。2群のばらつき（分散）そのものを比較する場合と、3群以上の平均を同時に比較する分散分析（Analysis of Variance, ANOVA）の中で使われる場合があります。

import numpy as np
from scipy import stats

def judge(p, alpha=0.05):
    return "有意差あり" if p < alpha else "有意差なし"

print("=" * 50)
print("5) F検定（2群の分散の比較）")
print("=" * 50)
# 2つの製造ラインの製品重量のばらつき（分散）を比較
line_a = np.array([100.2, 99.8, 100.5, 99.6, 100.1, 100.3, 99.9, 100.0, 99.7, 100.4])
line_b = np.array([100.1, 98.5, 101.3, 99.0, 100.8, 98.9, 101.5, 99.2, 100.6, 99.5])

var_a = np.var(line_a, ddof=1)
var_b = np.var(line_b, ddof=1)
f_stat = var_b / var_a  # 分散が大きい方を分子に
df1 = len(line_b) - 1
df2 = len(line_a) - 1
p_value_f = 2 * min(stats.f.cdf(f_stat, df1, df2), 1 - stats.f.cdf(f_stat, df1, df2))

print(f"ラインA分散: {var_a:.5f}, ラインB分散: {var_b:.5f}")
print(f"F統計量 (分散B/分散A): {f_stat:.4f}")
print(f"自由度: ({df1}, {df2})")
print(f"p値（両側）: {p_value_f:.4f}")
print(f"判定 (alpha=0.05): {judge(p_value_f)}")

print()
print("=" * 50)
print("6) 一元配置分散分析（F検定の応用：3群以上の平均比較）")
print("=" * 50)
# 3つの肥料（肥料A・B・C）による植物成長量の比較
fertilizer_a = np.array([12, 14, 13, 15, 12])
fertilizer_b = np.array([18, 17, 19, 16, 18])
fertilizer_c = np.array([22, 24, 21, 23, 25])

f_stat_anova, p_value_anova = stats.f_oneway(fertilizer_a, fertilizer_b, fertilizer_c)
print(f"肥料A平均: {np.mean(fertilizer_a):.2f}")
print(f"肥料B平均: {np.mean(fertilizer_b):.2f}")
print(f"肥料C平均: {np.mean(fertilizer_c):.2f}")
print(f"F統計量: {f_stat_anova:.4f}")
print(f"p値: {p_value_anova:.3e}")
print(f"判定 (alpha=0.05): {judge(p_value_anova)}")

実行結果:

==================================================
5) F検定（2群の分散の比較）
==================================================
ラインA分散: 0.09167, ラインB分散: 1.14044
F統計量 (分散B/分散A): 12.4412
自由度: (9, 9)
p値（両側）: 0.0009
判定 (alpha=0.05): 有意差あり

==================================================
6) 一元配置分散分析（F検定の応用：3群以上の平均比較）
==================================================
肥料A平均: 13.20
肥料B平均: 17.60
肥料C平均: 23.00
F統計量: 65.7091
p値: 3.431e-07
判定 (alpha=0.05): 有意差あり

6. 多重検定問題とボンフェローニ補正

6.1 多重検定問題とは

有意水準 $\alpha = 0.05$ で1回だけ検定を行う場合、帰無仮説が真であるにもかかわらず誤って棄却してしまう確率（第一種の過誤）はちょうど5%です。しかし、同じデータや似た条件で複数回の検定を繰り返すと、少なくとも1回は偶然「有意」という結果が出てしまう確率が急激に高まります。これを多重検定問題（Multiple Comparisons Problem）と呼びます。

$n$ 回の独立な検定をすべて有意水準 $\alpha$ で行うとき、少なくとも1つが偶然有意になってしまう確率（族全体の第一種の過誤率、Family-Wise Error Rate）は次のようになります。

$$P(\text{少なくとも1つの偽陽性}) = 1 - (1-\alpha)^n$$

例えば $\alpha=0.05$ で20回検定を行うと、この確率は約64%にも達します。「たくさん検定すれば、何かは必ず有意になる」というのは統計的な罠であり、機械学習における特徴量選択やA/Bテストの多変量比較でも頻繁に問題になります。

6.2 ボンフェローニ補正（Bonferroni Correction）

ボンフェローニ補正は、多重検定問題に対する最も単純で保守的な対処法です。$n$ 回の検定を行う場合、個々の検定に用いる有意水準を $\alpha$ から $\alpha/n$ に厳しくすることで、族全体の第一種の過誤率を元の $\alpha$ 以下に抑えます。

$$\alpha_{\text{補正後}} = \frac{\alpha}{n}$$

あるいは、個々のp値を $n$ 倍した調整済みp値を元の $\alpha$ と比較する、同等の方法もよく使われます。

6.3 Pythonでの実装：多重検定問題の実演とボンフェローニ補正

帰無仮説が真である（実際には差がない）データに対して20回の検定を繰り返し、補正なしでは偽陽性が発生すること、そしてボンフェローニ補正がそれを正しく抑制することを確認します。

import numpy as np
from scipy import stats

np.random.seed(3)

print("=== 多重検定問題の実演：帰無仮説が真であっても偶然有意になる ===")
# 20個の独立な検定を行う（すべて帰無仮説 H0: 平均差なし が真であるデータ）
n_tests = 20
alpha = 0.05
p_values = []

for i in range(n_tests):
    # 2群とも同じ分布（平均差なし）からサンプリング
    group_x = np.random.normal(50, 10, size=25)
    group_y = np.random.normal(50, 10, size=25)
    _, p = stats.ttest_ind(group_x, group_y)
    p_values.append(p)

p_values = np.array(p_values)
n_significant_uncorrected = np.sum(p_values < alpha)

print(f"実施した検定数: {n_tests}（すべて帰無仮説が真）")
print(f"有意水準 alpha: {alpha}")
print(f"補正なしで有意(p<0.05)となった検定数: {n_significant_uncorrected}")
print(f"少なくとも1つの偽陽性が出る理論的な確率: 1 - (1-{alpha})^{n_tests} = {1 - (1-alpha)**n_tests:.4f}")

print()
print("=== ボンフェローニ補正（Bonferroni Correction）の適用 ===")
bonferroni_alpha = alpha / n_tests
n_significant_bonferroni = np.sum(p_values < bonferroni_alpha)

print(f"補正後の有意水準: {alpha}/{n_tests} = {bonferroni_alpha:.5f}")
print(f"ボンフェローニ補正後に有意となった検定数: {n_significant_bonferroni}")

print()
print("=== scipy/statsmodels相当のボンフェローニ補正関数（自作） ===")
def bonferroni_correction(p_values, alpha=0.05):
    """
    ボンフェローニ補正を適用する

    Parameters:
    -----------
    p_values : array-like
        検定ごとのp値のリスト
    alpha : float
        補正前の有意水準

    Returns:
    --------
    corrected_alpha : float
        補正後の有意水準
    rejected : ndarray of bool
        各検定を棄却する（有意とする）かどうか
    adjusted_p_values : ndarray
        補正後のp値（元のp値 * 検定数、上限1.0）
    """
    p_values = np.asarray(p_values)
    n = len(p_values)
    corrected_alpha = alpha / n
    adjusted_p_values = np.minimum(p_values * n, 1.0)
    rejected = p_values < corrected_alpha
    return corrected_alpha, rejected, adjusted_p_values

corrected_alpha, rejected, adjusted_p = bonferroni_correction(p_values, alpha=0.05)
print(f"補正後の有意水準: {corrected_alpha:.5f}")
print(f"棄却された検定の数: {np.sum(rejected)}")
print(f"調整済みp値の最小値: {np.min(adjusted_p):.4f}")

実行結果:

=== 多重検定問題の実演：帰無仮説が真であっても偶然有意になる ===
実施した検定数: 20（すべて帰無仮説が真）
有意水準 alpha: 0.05
補正なしで有意(p<0.05)となった検定数: 2
少なくとも1つの偽陽性が出る理論的な確率: 1 - (1-0.05)^20 = 0.6415

=== ボンフェローニ補正（Bonferroni Correction）の適用 ===
補正後の有意水準: 0.05/20 = 0.00250
ボンフェローニ補正後に有意となった検定数: 0

=== scipy/statsmodels相当のボンフェローニ補正関数（自作） ===
補正後の有意水準: 0.00250
棄却された検定の数: 0
調整済みp値の最小値: 0.9014

7. まとめと次のステップ

この章では、限られた標本から母集団の性質を推し量る統計的推論の中核である、推定と仮説検定について学びました。

次のステップ

推定と検定という統計的推論の基盤を身につけたら、次はベイズ統計や、機械学習アルゴリズムそのものへの統計学の応用へと進んでいきます。本シリーズの続く章もあわせて学習し、統計的思考を機械学習の実践に活かしていきましょう。

練習問題

import numpy as np
from scipy import stats

data = np.array([48.2, 51.5, 49.8, 50.3, 52.1, 47.9,
                  50.7, 49.5, 51.2, 48.8, 50.0, 51.8])

n = len(data)
mean = np.mean(data)
var = np.var(data, ddof=1)
std = np.std(data, ddof=1)
se = std / np.sqrt(n)

ci_lower, ci_upper = stats.t.interval(0.95, n - 1, loc=mean, scale=se)

print(f"標本サイズ: {n}")
print(f"標本平均: {mean:.3f} mm")
print(f"不偏分散: {var:.3f}")
print(f"不偏標準偏差: {std:.3f} mm")
print(f"標準誤差: {se:.3f} mm")
print(f"95%信頼区間: [{ci_lower:.3f}, {ci_upper:.3f}] mm")

import numpy as np
from scipy import stats

before_method = np.array([62, 65, 70, 68, 72, 66, 69, 71, 64, 67])
new_method = np.array([70, 75, 78, 74, 80, 76, 77, 79, 73, 75])

# 2つの独立したグループの平均比較なので、対応のない2標本t検定を用いる
t_stat, p_value = stats.ttest_ind(before_method, new_method)

print(f"従来法の平均成約件数: {np.mean(before_method):.2f}")
print(f"新法の平均成約件数: {np.mean(new_method):.2f}")
print(f"t統計量: {t_stat:.4f}")
print(f"p値: {p_value:.3e}")

alpha = 0.05
if p_value < alpha:
    print("結論: 統計的に有意な差がある")
else:
    print("結論: 統計的に有意な差があるとは言えない")

import numpy as np

p_values = np.array([0.001, 0.012, 0.024, 0.038, 0.049])
alpha = 0.05
n_tests = len(p_values)

corrected_alpha = alpha / n_tests
adjusted_p_values = np.minimum(p_values * n_tests, 1.0)
rejected = p_values < corrected_alpha

print(f"補正後の有意水準: {alpha}/{n_tests} = {corrected_alpha}")
for i, (p, adj, rej) in enumerate(zip(p_values, adjusted_p_values, rejected), start=1):
    print(f"検定{i}: p={p:.3f}, 調整済みp値={adj:.3f}, "
          f"判定={'有意' if rej else '非有意'}")