第5章：機械学習への応用

イントロダクション

ここまでの4つの章で、記述統計と確率の基礎(第1章)、確率分布(第2章)、統計的推定と仮説検定(第3章)、そしてベイズ統計(第4章)を学んできました。本章はシリーズの最終章として、これらの統計学の知識が機械学習の具体的なアルゴリズムの中でどのように使われているかを、実装を通して確認していきます。

機械学習のアルゴリズムの多くは、統計学の概念の上に成り立っています。線形回帰(Linear Regression)は最小二乗法という統計的な最適化手法そのものであり、ロジスティック回帰(Logistic Regression)は最尤推定という統計的推定の枠組みで学習されます。ナイーブベイズ分類器(Naive Bayes Classifier)は第1章で学んだベイズの定理を分類問題に直接応用したものであり、ガウス過程(Gaussian Process)は確率分布とベイズ統計の考え方を回帰問題における不確実性の定量化へと発展させたものです。さらに、学習したモデルの性能を比較するモデル評価や、施策の効果を検証するA/Bテストにも、第3章で学んだ統計的検定の考え方が直接応用されます。

本章を通じて、「統計学を学ぶこと」と「機械学習を使いこなすこと」が地続きであることを実感していただければと思います。

本章で扱うアルゴリズムと、それぞれの統計的な基盤を一覧にまとめました。

1. 線形回帰と最小二乗法の統計的解釈

1.1 最小二乗法とは

最小二乗法(Least Squares Method、Ordinary Least Squares: OLS)は、観測値と予測値の差である残差(Residual)の二乗和を最小にするようにモデルのパラメータを決定する手法です。単回帰モデルは次のように表されます。

$$y_i = \beta_0 + \beta_1 x_i + \varepsilon_i$$

ここで、$\beta_0$は切片、$\beta_1$は傾き、$\varepsilon_i$は誤差項です。最小二乗法では、次の残差二乗和(RSS: Residual Sum of Squares)を最小にする$\beta_0, \beta_1$を求めます。

$$\text{RSS}(\beta_0, \beta_1) = \sum_{i=1}^{n}(y_i - \beta_0 - \beta_1 x_i)^2$$

この最適化問題は解析的に解くことができ、scikit-learnのLinearRegressionは内部でこの解析解(または効率的な数値計算)を用いてパラメータを推定します。

1.2 統計的な仮定と係数の不確実性

線形回帰を単なる「最適化」ではなく「統計的推論」として扱うためには、誤差項$\varepsilon_i$についていくつかの仮定を置きます。

• 線形性: $y$と$x$の関係が線形パラメータで表現できる
• 独立性: 誤差$\varepsilon_i$が互いに独立である
• 等分散性(Homoscedasticity): 誤差の分散が$x$の値によらず一定である
• 正規性: 誤差が正規分布に従う(仮説検定や信頼区間の構成に必要)

これらの仮定のもとでは、最小二乗推定量は不偏性を持ち、線形不偏推定量の中で分散が最小になる(ガウス・マルコフの定理)ことが知られています。さらに正規性の仮定を置くと、係数$\hat{\beta}_1$の標準誤差$\text{SE}(\hat{\beta}_1)$を用いて、第3章で学んだのと同様の枠組みでt検定や信頼区間を構成できます。

「傾き$\beta_1$が0である」という帰無仮説に対する検定統計量は次の通りです。

$$t = \frac{\hat{\beta}_1}{\text{SE}(\hat{\beta}_1)}$$

この$t$は自由度$n-p-1$(サンプルサイズ$n$から、切片を含むパラメータ数$p+1$を引いたもの)のt分布に従います。p値が小さければ、説明変数$x$が目的変数$y$の予測に統計的に有意に寄与していると判断できます。

1.3 Pythonでの実装

学習時間から試験得点を予測する単回帰モデルを、統計的な解釈まで含めて実装してみましょう。

import numpy as np
from scipy import stats
from sklearn.linear_model import LinearRegression

# 実験データを生成: 学習時間(時間) と 試験得点(点)
np.random.seed(42)
n = 30
study_hours = np.random.uniform(1, 10, n)
true_intercept, true_slope = 50, 4.5
noise = np.random.normal(0, 6, n)
scores = true_intercept + true_slope * study_hours + noise

X = study_hours.reshape(-1, 1)
y = scores

# scikit-learnで最小二乗法によるフィット
model = LinearRegression()
model.fit(X, y)

intercept = model.intercept_
slope = model.coef_[0]
y_pred = model.predict(X)

# 決定係数R^2
r_squared = model.score(X, y)

# 残差の分析
residuals = y - y_pred
n_samples, n_features = X.shape
dof = n_samples - n_features - 1  # 自由度(切片+傾きの2パラメータを推定)
residual_std_error = np.sqrt(np.sum(residuals**2) / dof)

# 係数の標準誤差(計画行列を使った標準的な公式)
X_design = np.column_stack([np.ones(n_samples), X])
cov_matrix = residual_std_error**2 * np.linalg.inv(X_design.T @ X_design)
se_intercept = np.sqrt(cov_matrix[0, 0])
se_slope = np.sqrt(cov_matrix[1, 1])

# 傾きに対する95%信頼区間とt検定(帰無仮説: 傾き=0)
t_critical = stats.t.ppf(0.975, dof)
ci_slope = (slope - t_critical * se_slope, slope + t_critical * se_slope)
t_stat_slope = slope / se_slope
p_value_slope = 2 * (1 - stats.t.cdf(abs(t_stat_slope), dof))

print(f"切片 (Intercept): {intercept:.4f}")
print(f"傾き (Slope): {slope:.4f}")
print(f"決定係数 R^2: {r_squared:.4f}")
print(f"残差標準誤差: {residual_std_error:.4f}")
print(f"\n傾きの標準誤差: {se_slope:.4f}")
print(f"傾きの95%信頼区間: [{ci_slope[0]:.4f}, {ci_slope[1]:.4f}]")
print(f"傾きのt統計量: {t_stat_slope:.4f}")
print(f"傾きのp値: {p_value_slope:.2e}")

実行結果:

切片 (Intercept): 53.5025
傾き (Slope): 3.5713
決定係数 R^2: 0.7671
残差標準誤差: 5.0899

傾きの標準誤差: 0.3719
傾きの95%信頼区間: [2.8095, 4.3330]
傾きのt統計量: 9.6036
傾きのp値: 2.33e-10

2. ロジスティック回帰と最尤推定

2.1 なぜ線形回帰では不十分なのか

目的変数が「合格/不合格」のような2値である分類問題に、そのまま線形回帰を適用すると、予測値が0未満や1超になってしまい、確率として解釈できません。そこでロジスティック回帰(Logistic Regression)では、線形結合をシグモイド関数(Sigmoid Function)に通すことで、出力を$[0,1]$の範囲に収めます。

$$p_i = P(y_i=1|x_i) = \frac{1}{1+e^{-(\beta_0+\beta_1 x_i)}}$$

2.2 最尤推定によるパラメータ推定

ロジスティック回帰のパラメータは、最小二乗法ではなく最尤推定(Maximum Likelihood Estimation, MLE)によって求められます。目的変数$y_i$がベルヌーイ分布に従うと仮定すると、尤度関数は次のように書けます。

$$L(\beta) = \prod_{i=1}^{n} p_i^{y_i}(1-p_i)^{1-y_i}$$

対数を取ると、扱いやすい対数尤度(Log-Likelihood)になります。

$$\ell(\beta) = \sum_{i=1}^{n}\left[y_i \log p_i + (1-y_i)\log(1-p_i)\right]$$

この対数尤度を最大にする$\beta$を求める問題には解析解が存在しないため、scikit-learnのLogisticRegressionは勾配ベースの数値最適化(既定ではL-BFGS法)を用いて反復的に$\beta$を推定します。この対数尤度の符号を反転させたものは対数損失(Log Loss、交差エントロピー損失)と呼ばれ、多くの分類モデルの損失関数として使われています。

2.3 係数の解釈:対数オッズとオッズ比

ロジスティック回帰の係数$\beta_1$は、説明変数$x$が1単位増えたときの対数オッズ(Log-Odds)の変化量として解釈できます。$\exp(\beta_1)$を取るとオッズ比(Odds Ratio)になり、「$x$が1単位増えると、成功のオッズが何倍になるか」という直感的な形で解釈できます。

2.4 Pythonでの実装

import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import log_loss

# 実験データを生成: 勉強時間から合否を予測する2値分類問題
np.random.seed(42)
n = 100
study_hours = np.random.uniform(0, 10, n)
# 真のロジスティックモデル: P(合格=1) = sigmoid(-4 + 1.0 * study_hours)
true_logit = -4 + 1.0 * study_hours
true_prob = 1 / (1 + np.exp(-true_logit))
pass_fail = (np.random.uniform(0, 1, n) < true_prob).astype(int)

X = study_hours.reshape(-1, 1)
y = pass_fail

# scikit-learnで最尤推定によるロジスティック回帰
# C=1e6として正則化の影響をほぼ無効化し、純粋な最尤推定に近づける
model = LogisticRegression(C=1e6)
model.fit(X, y)

intercept = model.intercept_[0]
coef = model.coef_[0][0]

# 予測確率と対数尤度(log-likelihood)
proba = model.predict_proba(X)[:, 1]
# log_lossは負の対数尤度の平均なので、n倍して符号を反転させると対数尤度が得られる
log_likelihood = -log_loss(y, proba, normalize=True) * n

accuracy = model.score(X, y)

print(f"切片 (Intercept): {intercept:.4f}")
print(f"係数 (Coefficient): {coef:.4f}")
print(f"対数尤度 (Log-Likelihood): {log_likelihood:.4f}")
print(f"訓練データに対する正解率: {accuracy:.4f}")

# 係数の解釈: 勉強時間が1時間増えるとオッズが何倍になるか
odds_ratio = np.exp(coef)
print(f"\nオッズ比 exp(係数): {odds_ratio:.4f}")
print(f"→ 勉強時間が1時間増えると、合格のオッズは約{odds_ratio:.2f}倍になる")

# 具体的な勉強時間での合格確率を予測
test_hours = np.array([[2], [4], [6], [8]])
test_proba = model.predict_proba(test_hours)[:, 1]
print("\n勉強時間ごとの予測合格確率:")
for h, p in zip(test_hours.ravel(), test_proba):
    print(f"  {h:.0f}時間: {p:.4f}")

実行結果:

切片 (Intercept): -5.2407
係数 (Coefficient): 1.2589
対数尤度 (Log-Likelihood): -24.4132
訓練データに対する正解率: 0.8800

オッズ比 exp(係数): 3.5214
→ 勉強時間が1時間増えると、合格のオッズは約3.52倍になる

勉強時間ごとの予測合格確率:
  2時間: 0.0616
  4時間: 0.4489
  6時間: 0.9099
  8時間: 0.9921

3. ナイーブベイズ分類器の実装

3.1 ベイズの定理を分類問題に応用する

第1章で学んだベイズの定理を思い出しましょう。特徴量$x_1,\ldots,x_n$が観測されたときにクラス$y$である事後確率は、次のように書けます。

$$P(y|x_1,\ldots,x_n) = \frac{P(y)\,P(x_1,\ldots,x_n|y)}{P(x_1,\ldots,x_n)}$$

ここで問題になるのが、特徴量の数が増えると同時分布$P(x_1,\ldots,x_n|y)$の推定に必要なデータ量が指数的に増えてしまうことです。ナイーブベイズ分類器(Naive Bayes Classifier)は、「クラス$y$が与えられたもとで、各特徴量は互いに条件付き独立である」という単純化(ナイーブ、Naive)な仮定を置くことで、この問題を回避します。

$$P(x_1,\ldots,x_n|y) \approx \prod_{i=1}^{n}P(x_i|y)$$

この仮定のもとで、分類はクラスごとの事後確率(正確には、分母$P(x_1,\ldots,x_n)$はどのクラスでも共通なので無視した、事前確率と尤度の積)を比較し、最大のものを選ぶ問題に帰着します。

$$\hat{y} = \arg\max_{y}\; P(y)\prod_{i=1}^{n}P(x_i|y)$$

3.2 ガウシアンナイーブベイズ

特徴量が連続値の場合、各クラス内での各特徴量の分布を正規分布と仮定するのがガウシアンナイーブベイズ(Gaussian Naive Bayes)です。訓練データから、クラス$y$ごとに各特徴量の平均$\mu_{y,i}$と分散$\sigma_{y,i}^2$を推定し、それを正規分布の確率密度関数に代入して尤度$P(x_i|y)$を計算します。

3.3 Pythonでの実装

実測データであるIrisデータセット(あやめの花の測定データ、3品種の分類問題)を使って実装します。

import numpy as np
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB
from sklearn.metrics import accuracy_score

# Irisデータセット(あやめの花の測定データ、3クラス分類)を使用
iris = load_iris()
X, y = iris.data, iris.target
class_names = iris.target_names

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42, stratify=y
)

# ガウシアンナイーブベイズ分類器を学習
# 各クラス内で各特徴量が正規分布に従うと仮定し、特徴量間の独立性を仮定する
gnb = GaussianNB()
gnb.fit(X_train, y_train)

y_pred = gnb.predict(X_test)
accuracy = accuracy_score(y_test, y_pred)

print(f"テストデータでの正解率: {accuracy:.4f}")

# 学習されたクラスごとの事前確率(各クラスの訓練データ中の割合)
print("\nクラスごとの事前確率 P(クラス):")
for name, prior in zip(class_names, gnb.class_prior_):
    print(f"  {name}: {prior:.4f}")

# 1つ目の特徴量(がく片の長さ)についてクラスごとの平均と分散(ガウス分布のパラメータ)
print("\n特徴量0(がく片の長さ)のクラスごとの平均と分散:")
for i, name in enumerate(class_names):
    mean = gnb.theta_[i, 0]
    var = gnb.var_[i, 0]
    print(f"  {name}: 平均={mean:.4f}, 分散={var:.4f}")

# テストサンプル1件について、各クラスに属する事後確率を確認
sample = X_test[0].reshape(1, -1)
posterior = gnb.predict_proba(sample)[0]
true_label = class_names[y_test[0]]

print(f"\nテストサンプル1件の特徴量: {X_test[0]}")
print(f"真のクラス: {true_label}")
print("各クラスへの事後確率 P(クラス|特徴量):")
for name, p in zip(class_names, posterior):
    print(f"  {name}: {p:.4f}")

実行結果:

テストデータでの正解率: 0.9111

クラスごとの事前確率 P(クラス):
  setosa: 0.3333
  versicolor: 0.3333
  virginica: 0.3333

特徴量0(がく片の長さ)のクラスごとの平均と分散:
  setosa: 平均=4.9886, 分散=0.1033
  versicolor: 平均=5.9486, 分散=0.2408
  virginica: 平均=6.6829, 分散=0.4248

テストサンプル1件の特徴量: [7.3 2.9 6.3 1.8]
真のクラス: virginica
各クラスへの事後確率 P(クラス|特徴量):
  setosa: 0.0000
  versicolor: 0.0000
  virginica: 1.0000

4. ガウス過程による予測の不確実性推定

4.1 点推定だけでは足りない場面

これまでの線形回帰やロジスティック回帰は、基本的に「最も確からしい1つの値」を予測する点推定(Point Estimate)のモデルです。しかし実務では、「この予測はどれくらい確信できるのか」という不確実性の情報が重要になる場面が数多くあります。例えば、訓練データが少ない領域や、訓練データの分布から大きく外れた入力に対しては、予測の信頼度は本来低くなるはずです。

ガウス過程(Gaussian Process, GP)は、関数そのものに対する確率分布を定義することで、予測値の平均だけでなく、その不確実性(分散)も自然な形で出力できる回帰手法です。

4.2 ガウス過程の考え方

ガウス過程は、任意の有限個の入力点$x_1,\ldots,x_m$に対する関数値$f(x_1),\ldots,f(x_m)$が、常に多変量正規分布に従うという性質を持つ確率過程として定義されます。

$$f(x) \sim \mathcal{GP}(m(x), k(x, x'))$$

ここで$m(x)$は平均関数(多くの場合0とする)、$k(x,x')$はカーネル関数(Kernel Function)と呼ばれる共分散関数で、2点$x$と$x'$がどれだけ「似ている」かを表現します。本章では、最も広く使われるRBFカーネル(動径基底関数カーネル、Radial Basis Function Kernel)を使用します。近い入力どうしは強く相関し、遠い入力どうしは相関が弱くなるという滑らかさの仮定をカーネルが表現しています。

観測データ$D$を条件づけると、新しい入力点$x_*$における予測分布もまた正規分布になり、その平均が点推定、標準偏差が予測の不確実性を表します。訓練データに近い$x_*$では標準偏差が小さく(確信度が高く)、訓練データから離れるほど標準偏差が大きくなる(確信度が低くなる)という性質が、ガウス過程の大きな特徴です。

4.3 Pythonでの実装

少数の観測点しかないセンサーデータを想定し、観測点の近くと遠くで予測の不確実性がどう変化するかを確認します。

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

# 少数の観測点しかない状況を再現(センサーの点検データを想定)
np.random.seed(42)
X_train = np.array([[1.0], [2.0], [3.5], [6.0], [7.5]])
def true_function(x):
    return np.sin(x) * 2 + 0.3 * x

y_train = true_function(X_train).ravel() + np.random.normal(0, 0.15, X_train.shape[0])

# カーネル: RBFカーネル(滑らかさを表現)+ WhiteKernel(観測ノイズを表現)
# 探索範囲を適度に制限し、少数データでもハイパーパラメータが妥当な範囲に収まるようにする
kernel = RBF(length_scale=1.5, length_scale_bounds=(0.5, 5.0)) \
    + WhiteKernel(noise_level=0.05, noise_level_bounds=(1e-3, 0.5))

gpr = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, random_state=42)
gpr.fit(X_train, y_train)

print(f"学習後のカーネルパラメータ: {gpr.kernel_}")
print(f"対数周辺尤度: {gpr.log_marginal_likelihood(gpr.kernel_.theta):.4f}")

# 観測点付近と観測点から大きく離れた点で予測の不確実性(標準偏差)を比較
X_test = np.array([[1.5], [4.5], [7.0], [12.0]])
y_mean, y_std = gpr.predict(X_test, return_std=True)

print("\n予測点ごとの平均と標準偏差(不確実性):")
for x, m, s in zip(X_test.ravel(), y_mean, y_std):
    print(f"  x={x:.1f}: 予測平均={m:.4f}, 標準偏差={s:.4f}, "
          f"95%予測区間=[{m - 1.96*s:.4f}, {m + 1.96*s:.4f}]")

print("\n観測データが近くにある点(x=1.5, 4.5)は標準偏差が小さく、")
print("観測データの範囲外に大きく外挿した点(x=12.0)は標準偏差が大きくなる傾向を確認できる。")

実行結果:

学習後のカーネルパラメータ: RBF(length_scale=1.53) + WhiteKernel(noise_level=0.5)
対数周辺尤度: -13.1404

予測点ごとの平均と標準偏差(不確実性):
  x=1.5: 予測平均=1.8016, 標準偏差=0.8472, 95%予測区間=[0.1412, 3.4621]
  x=4.5: 予測平均=0.3935, 標準偏差=0.9498, 95%予測区間=[-1.4681, 2.2551]
  x=7.0: 予測平均=2.5160, 標準偏差=0.8721, 95%予測区間=[0.8067, 4.2253]
  x=12.0: 予測平均=0.0369, 標準偏差=1.2247, 95%予測区間=[-2.3635, 2.4373]

観測データが近くにある点(x=1.5, 4.5)は標準偏差が小さく、
観測データの範囲外に大きく外挿した点(x=12.0)は標準偏差が大きくなる傾向を確認できる。

5. モデル評価と統計的検定

5.1 なぜ1回の評価では不十分なのか

訓練データとテストデータを1回だけ分割してモデルを評価すると、その分割の仕方(たまたま簡単なサンプルがテストデータに多く含まれた、など)によって評価結果が偶然左右されてしまうことがあります。交差検証(Cross Validation)は、データを$k$個に分割し、そのうち1つをテスト用、残りを訓練用として使うことを$k$回繰り返す手法です。これにより、評価結果のばらつき(標準偏差)も含めてモデルの性能を把握できます。

5.2 モデル間の性能差を統計的に検定する

2つのモデルの交差検証スコアを比較する際、単に平均値の大小を見るだけでは、その差が偶然のばらつきによるものなのか、統計的に意味のある差なのかを判断できません。ここで第3章の仮説検定の考え方が役立ちます。同一の分割(fold)で評価した2モデルのスコアは対応関係にあるため、対応のあるt検定(Paired t-test)を用いて、スコア差の平均が0であるという帰無仮説を検定します。

5.3 Pythonでの実装

通常の線形回帰と、正則化を加えたRidge回帰の性能を10分割交差検証で比較し、その差が統計的に有意かどうかを検定します。

import numpy as np
from scipy import stats
from sklearn.datasets import make_regression
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.model_selection import KFold, cross_val_score

# 回帰用の合成データ(特徴量間に多少の相関を持たせる)
X, y = make_regression(
    n_samples=200, n_features=15, n_informative=8,
    noise=15.0, random_state=42
)

# 通常の線形回帰と、正則化ありのRidge回帰を比較する
model_a = LinearRegression()
model_b = Ridge(alpha=10.0)

# 同一の分割(KFold)を両モデルに使うことで、対応のある(paired)比較にする
kf = KFold(n_splits=10, shuffle=True, random_state=42)

scores_a = cross_val_score(model_a, X, y, cv=kf, scoring="neg_mean_squared_error")
scores_b = cross_val_score(model_b, X, y, cv=kf, scoring="neg_mean_squared_error")

mse_a = -scores_a
mse_b = -scores_b

print("=== 10分割交差検証によるMSE(平均二乗誤差) ===")
print(f"線形回帰    : 平均MSE={mse_a.mean():.4f}, 標準偏差={mse_a.std():.4f}")
print(f"Ridge回帰   : 平均MSE={mse_b.mean():.4f}, 標準偏差={mse_b.std():.4f}")

# 対応のあるt検定(paired t-test): 同じ分割で評価した2モデルのスコア差が
# 統計的に有意かどうかを検定する(帰無仮説: 平均差は0)
t_stat, p_value = stats.ttest_rel(mse_a, mse_b)

print(f"\n=== 対応のあるt検定(モデル間のMSE差の検定) ===")
print(f"t統計量: {t_stat:.4f}")
print(f"p値: {p_value:.4f}")

alpha = 0.05
if p_value < alpha:
    print(f"p値 < {alpha} のため、2モデルの性能差は統計的に有意であると判断できる")
else:
    print(f"p値 >= {alpha} のため、2モデルの性能差が統計的に有意であるとは言えない")

実行結果:

=== 10分割交差検証によるMSE(平均二乗誤差) ===
線形回帰    : 平均MSE=244.2298, 標準偏差=87.0243
Ridge回帰   : 平均MSE=357.3553, 標準偏差=96.6626

=== 対応のあるt検定(モデル間のMSE差の検定) ===
t統計量: -2.7153
p値: 0.0238
p値 < 0.05 のため、2モデルの性能差は統計的に有意であると判断できる

6. A/Bテストの統計的手法

6.1 頻度論的アプローチによるA/Bテスト

第4章では、A/Bテストをベイズ統計の枠組みで扱い、「B案がA案より優れている確率」を直接計算しました。本節では対照的に、統計学で古くから使われてきた2標本t検定(Two-Sample t-test)による頻度論的なA/Bテストを実装します。連続値の指標(購入金額、滞在時間など)を2群で比較する際に広く使われる手法です。

2群の分散が等しいと仮定できない場合に頑健なWelchのt検定(Welch's t-test)を使うのが実務上安全な選択です。検定統計量は次の式で計算されます。

$$t = \frac{\bar{x}_A - \bar{x}_B}{\sqrt{s_A^2/n_A + s_B^2/n_B}}$$

ここで$\bar{x}_A, \bar{x}_B$は各群の標本平均、$s_A^2, s_B^2$は各群の標本分散、$n_A, n_B$は各群のサンプルサイズです。自由度はWelch–Satterthwaiteの式という近似式で計算されます。

6.2 効果量:統計的有意性だけでは分からないこと

p値は「差が偶然によるものかどうか」を判断する材料にはなりますが、「差がどれだけ実質的に大きいか」は教えてくれません。サンプルサイズが非常に大きい場合、実務上ほとんど意味のない小さな差でも統計的に有意になってしまうことがあります。そこで効果量(Effect Size)、特にCohenのd(Cohen's d)を併用し、差の実質的な大きさを評価します。

$$d = \frac{\bar{x}_B - \bar{x}_A}{s_{\text{pooled}}}$$

一般的な目安として、$|d|\approx 0.2$は小さい効果、$0.5$は中程度の効果、$0.8$以上は大きい効果とされています(あくまで目安であり、分野や文脈によって解釈は異なります)。

6.3 Pythonでの実装

ECサイトの新レイアウト(B案)が平均購入金額を増やすかどうかを、Welchのt検定で検証します。標準的な公式を手計算で実装し、scipy.statsの結果と照合します。

import numpy as np
from scipy import stats

# A/Bテスト: ECサイトの新レイアウト(B案)が平均購入金額を増やすかを検証する
# A案(既存レイアウト)とB案(新レイアウト)それぞれの購入者の購入金額(円)を観測したと想定
np.random.seed(42)
n_a, n_b = 60, 65
group_a = np.random.normal(loc=4200, scale=900, size=n_a)   # A案: 平均4200円
group_b = np.random.normal(loc=4550, scale=950, size=n_b)   # B案: 平均4550円

def welch_t_test(x, y):
    """
    Welchのt検定(等分散を仮定しない2標本t検定)を手計算で実装する

    Parameters:
    -----------
    x, y : ndarray
        比較する2つの群のデータ

    Returns:
    --------
    t_stat, dof, p_value : t統計量、自由度、両側検定のp値
    """
    mean_x, mean_y = np.mean(x), np.mean(y)
    var_x, var_y = np.var(x, ddof=1), np.var(y, ddof=1)
    n_x, n_y = len(x), len(y)

    se_diff = np.sqrt(var_x / n_x + var_y / n_y)
    t_stat = (mean_x - mean_y) / se_diff

    # Welch-Satterthwaiteの式による自由度の近似
    dof = (var_x / n_x + var_y / n_y) ** 2 / (
        (var_x / n_x) ** 2 / (n_x - 1) + (var_y / n_y) ** 2 / (n_y - 1)
    )

    p_value = 2 * (1 - stats.t.cdf(abs(t_stat), dof))
    return t_stat, dof, p_value

t_stat, dof, p_value = welch_t_test(group_a, group_b)

print("=== 手計算によるWelchのt検定 ===")
print(f"A案: 平均={group_a.mean():.2f}円, 標準偏差={group_a.std(ddof=1):.2f}円, n={n_a}")
print(f"B案: 平均={group_b.mean():.2f}円, 標準偏差={group_b.std(ddof=1):.2f}円, n={n_b}")
print(f"\nt統計量: {t_stat:.4f}")
print(f"自由度(近似): {dof:.2f}")
print(f"p値(両側検定): {p_value:.4f}")

# scipy.stats.ttest_indで検算(equal_var=Falseを指定するとWelchの検定になる)
t_stat_scipy, p_value_scipy = stats.ttest_ind(group_a, group_b, equal_var=False)
print(f"\n=== scipy.stats.ttest_ind による検算 ===")
print(f"t統計量: {t_stat_scipy:.4f}")
print(f"p値: {p_value_scipy:.4f}")

# 効果量(Cohen's d): 統計的有意性だけでなく、差の実質的な大きさを評価する
pooled_std = np.sqrt(((n_a - 1) * group_a.var(ddof=1) + (n_b - 1) * group_b.var(ddof=1)) / (n_a + n_b - 2))
cohens_d = (group_b.mean() - group_a.mean()) / pooled_std

print(f"\n効果量(Cohen's d): {cohens_d:.4f}")

alpha = 0.05
print(f"\n=== 結論 ===")
if p_value < alpha:
    print(f"p値({p_value:.4f}) < 有意水準({alpha}) のため、A案とB案の平均購入金額には")
    print("統計的に有意な差があると判断できる。B案への切り替えを検討する価値がある。")
else:
    print(f"p値({p_value:.4f}) >= 有意水準({alpha}) のため、統計的に有意な差があるとは言えない。")

実行結果:

=== 手計算によるWelchのt検定 ===
A案: 平均=4060.81円, 標準偏差=817.67円, n=60
B案: 平均=4553.63円, 標準偏差=906.21円, n=65

t統計量: -3.1960
自由度(近似): 122.94
p値(両側検定): 0.0018

=== scipy.stats.ttest_ind による検算 ===
t統計量: -3.1960
p値: 0.0018

効果量(Cohen's d): 0.5698

=== 結論 ===
p値(0.0018) < 有意水準(0.05) のため、A案とB案の平均購入金額には
統計的に有意な差があると判断できる。B案への切り替えを検討する価値がある。

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

この章では、これまでのシリーズで学んだ統計学の知識が、機械学習の代表的なアルゴリズムの中でどのように使われているかを、理論と実装の両面から確認しました。

学習目標の振り返り

この章の冒頭で示した学習目標を振り返ってみましょう。

• ✅ 機械学習アルゴリズムの統計的基礎を理解する — 第1節・第2節・第3節で線形回帰、ロジスティック回帰、ナイーブベイズの統計的な成り立ちを確認した
• ✅ 統計学の知識を機械学習に応用できる — 各節でscikit-learnとnumpy/scipyを組み合わせた実装を通じて確認した
• ✅ 予測の不確実性を定量化できる — 第4節でガウス過程による予測分散の計算を実装した
• ✅ モデル評価を統計的に行える — 第5節・第6節で交差検証・対応のあるt検定・2標本t検定・効果量を実装した

シリーズ全体の振り返り

本シリーズ「機械学習のための統計学入門」は、この第5章をもって完結です。第1章の記述統計と確率の基礎から始まり、第2章で確率分布の広がりを、第3章で推定と仮説検定という統計的推論の骨格を、第4章でベイズ統計という現代的な視点を学び、そして本章でそれらすべてが機械学習の実践にどうつながるかを確認してきました。統計学は決して機械学習の「前座」ではなく、モデルの振る舞いを理解し、結果を正しく解釈し、意思決定に説得力を持たせるための土台であり続けます。ここで得た統計的な視点を、今後さらに高度な機械学習の学習・実践にも活かしていってください。

練習問題

import numpy as np
from sklearn.linear_model import LinearRegression

# 広告費(万円)と売上(万円)のデータ
ad_spend = np.array([10, 15, 20, 25, 30, 35, 40, 45, 50, 55]).reshape(-1, 1)
sales = np.array([120, 145, 158, 190, 210, 225, 260, 275, 300, 320])

model = LinearRegression()
model.fit(ad_spend, sales)

r_squared = model.score(ad_spend, sales)
print(f"切片: {model.intercept_:.4f}")
print(f"傾き: {model.coef_[0]:.4f}")
print(f"決定係数 R^2: {r_squared:.4f}")

# 広告費42万円のときの売上予測
pred = model.predict(np.array([[42]]))
print(f"広告費42万円のときの売上予測: {pred[0]:.2f}万円")

# 出力:
# 切片: 74.7394
# 傾き: 4.4788
# 決定係数 R^2: 0.9965
# 広告費42万円のときの売上予測: 262.85万円

import numpy as np
from sklearn.datasets import load_wine
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB
from sklearn.metrics import accuracy_score

# Wineデータセット(ワインの化学成分から3つの品種を分類する)を使用
wine = load_wine()
X, y = wine.data, wine.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=0, stratify=y
)

gnb = GaussianNB()
gnb.fit(X_train, y_train)
y_pred = gnb.predict(X_test)

accuracy = accuracy_score(y_test, y_pred)
print(f"テストデータでの正解率: {accuracy:.4f}")

# 事前確率と、最初のテストサンプルの事後確率を確認
print("\nクラスごとの事前確率:")
for name, prior in zip(wine.target_names, gnb.class_prior_):
    print(f"  {name}: {prior:.4f}")

posterior = gnb.predict_proba(X_test[:1])[0]
posterior_str = ", ".join(f"{p:.4e}" for p in posterior)
print(f"\nテストサンプル1件目の事後確率: [{posterior_str}]")
print(f"予測クラス: {wine.target_names[y_pred[0]]}, 真のクラス: {wine.target_names[y_test[0]]}")

# 出力:
# テストデータでの正解率: 0.9630
#
# クラスごとの事前確率:
#   class_0: 0.3306
#   class_1: 0.4032
#   class_2: 0.2661
#
# テストサンプル1件目の事後確率: [1.0000e+00, 1.0339e-13, 1.9798e-38]
# 予測クラス: class_0, 真のクラス: class_0

import numpy as np
from scipy import stats

# A/Bテスト: メールマガジンの件名を2パターン(A/B)で配信し、開封後の滞在時間(秒)を比較
np.random.seed(1)
group_a = np.random.normal(loc=45, scale=12, size=40)  # A案(既存の件名)
group_b = np.random.normal(loc=52, scale=13, size=42)  # B案(新しい件名)

t_stat, p_value = stats.ttest_ind(group_a, group_b, equal_var=False)

pooled_std = np.sqrt(((len(group_a) - 1) * group_a.var(ddof=1)
                       + (len(group_b) - 1) * group_b.var(ddof=1))
                      / (len(group_a) + len(group_b) - 2))
cohens_d = (group_b.mean() - group_a.mean()) / pooled_std

print(f"A案: 平均滞在時間={group_a.mean():.2f}秒, n={len(group_a)}")
print(f"B案: 平均滞在時間={group_b.mean():.2f}秒, n={len(group_b)}")
print(f"\nt統計量: {t_stat:.4f}")
print(f"p値: {p_value:.4f}")
print(f"効果量(Cohen's d): {cohens_d:.4f}")

alpha = 0.05
if p_value < alpha:
    print(f"\np値 < {alpha} のため、統計的に有意な差があると判断できる")
else:
    print(f"\np値 >= {alpha} のため、統計的に有意な差があるとは言えない")

# 出力:
# A案: 平均滞在時間=43.92秒, n=40
# B案: 平均滞在時間=54.20秒, n=42
#
# t統計量: -3.9320
# p値: 0.0002
# 効果量(Cohen's d): 0.8683
#
# p値 < 0.05 のため、統計的に有意な差があると判断できる