バギングの理論(分散削減の数学的証明)を解説

機械学習において、単一のモデルでは高いバリアンス(分散)を持つ不安定な予測器(例: 決定木)が多く存在します。バギング(Bootstrap Aggregating)は、Breiman(1996)によって提案された、ブートストラップサンプリングと予測の平均化により、モデルの分散を系統的に削減するアンサンブル手法です。

バギングの効果を正確に理解するには、バイアス-バリアンス分解と、相関のある確率変数の平均の分散に関する数学的知識が必要です。本記事では、これらの理論を1行ずつ丁寧に導出し、OOB(Out-of-Bag)推定の理論も含めて解説します。最後にPythonで決定木のバギングをスクラッチ実装し、分散削減を実験的に検証します。

本記事の内容

  • アンサンブル学習の基本思想
  • バイアス-バリアンス分解の完全な導出
  • バギングのアルゴリズム
  • バギングによる分散削減の数学的証明(独立の場合と相関ありの場合)
  • OOB(Out-of-Bag)推定の理論
  • サブサンプリングとの比較
  • Pythonでの実装と実験的検証

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

アンサンブル学習の基本思想

アンサンブル学習とは、複数のモデル(弱学習器)を組み合わせることで、単一モデルよりも優れた性能を得る手法の総称です。基本的なアイデアは「三人寄れば文殊の知恵」と同じで、個々の不完全な予測を集約することで、全体として良い予測を得ます。

アンサンブル手法は大きく分けて2つのカテゴリがあります。

  1. 並列型(Parallel methods): 各モデルを独立に学習し、予測を平均化または多数決する。バギング、ランダムフォレストがこれに該当します。
  2. 逐次型(Sequential methods): 前のモデルの誤りを次のモデルが補正するように逐次的に学習する。ブースティング(AdaBoost, Gradient Boosting)がこれに該当します。

バギングは並列型に属し、分散の削減を主な目的としています。

バイアス-バリアンス分解

設定

データ生成過程を $y = f(\bm{x}) + \epsilon$ とします。ここで $f(\bm{x})$ は真の関数、$\epsilon$ は平均0、分散 $\sigma^2$ のノイズで、$\bm{x}$ と独立とします。

訓練データ $\mathcal{D} = \{(\bm{x}_i, y_i)\}_{i=1}^n$ は確率的に生成されるため、$\mathcal{D}$ に依存して学習される予測関数 $\hat{f}(\bm{x}; \mathcal{D})$ も確率変数です。以降、簡単のため $\hat{f}(\bm{x})$ と書きます。

期待予測誤差の分解

テスト点 $\bm{x}$ における期待予測誤差(MSE)を計算します。期待値は $\mathcal{D}$ の分布と $\epsilon$ の分布について取ります。

$$ \text{EPE}(\bm{x}) = E_{\mathcal{D}, \epsilon}\left[(y – \hat{f}(\bm{x}))^2\right] $$

まず $y = f(\bm{x}) + \epsilon$ を代入します。

$$ \text{EPE}(\bm{x}) = E_{\mathcal{D}, \epsilon}\left[(f(\bm{x}) + \epsilon – \hat{f}(\bm{x}))^2\right] $$

$\epsilon$ と $\hat{f}(\bm{x})$ は独立($\epsilon$ はテストデータのノイズ、$\hat{f}$ は訓練データに依存)なので、展開します。

$$ \begin{align} \text{EPE}(\bm{x}) &= E\left[(f(\bm{x}) – \hat{f}(\bm{x}))^2 + 2\epsilon(f(\bm{x}) – \hat{f}(\bm{x})) + \epsilon^2\right] \\ &= E_{\mathcal{D}}\left[(f(\bm{x}) – \hat{f}(\bm{x}))^2\right] + 2 E[\epsilon] \cdot E_{\mathcal{D}}[f(\bm{x}) – \hat{f}(\bm{x})] + E[\epsilon^2] \end{align} $$

$E[\epsilon] = 0$ なので中間項は消え、$E[\epsilon^2] = \text{Var}[\epsilon] = \sigma^2$ なので、

$$ \text{EPE}(\bm{x}) = E_{\mathcal{D}}\left[(f(\bm{x}) – \hat{f}(\bm{x}))^2\right] + \sigma^2 $$

第1項のさらなる分解

$E_{\mathcal{D}}[\hat{f}(\bm{x})]$ を $\bar{f}(\bm{x})$ と略記します。これは「$\hat{f}(\bm{x})$ の、異なる訓練データにわたる平均的な予測」です。

$f(\bm{x}) – \hat{f}(\bm{x})$ に $\bar{f}(\bm{x})$ を加えて引きます。

$$ f(\bm{x}) – \hat{f}(\bm{x}) = \underbrace{(f(\bm{x}) – \bar{f}(\bm{x}))}_{\text{バイアス成分}} + \underbrace{(\bar{f}(\bm{x}) – \hat{f}(\bm{x}))}_{\text{バリアンス成分}} $$

二乗して期待値を取ります。

$$ \begin{align} E_{\mathcal{D}}\left[(f(\bm{x}) – \hat{f}(\bm{x}))^2\right] &= E_{\mathcal{D}}\left[\left((f(\bm{x}) – \bar{f}(\bm{x})) + (\bar{f}(\bm{x}) – \hat{f}(\bm{x}))\right)^2\right] \\ &= (f(\bm{x}) – \bar{f}(\bm{x}))^2 + 2(f(\bm{x}) – \bar{f}(\bm{x})) E_{\mathcal{D}}[\bar{f}(\bm{x}) – \hat{f}(\bm{x})] + E_{\mathcal{D}}\left[(\bar{f}(\bm{x}) – \hat{f}(\bm{x}))^2\right] \end{align} $$

ここで交差項の期待値を計算します。

$$ E_{\mathcal{D}}[\bar{f}(\bm{x}) – \hat{f}(\bm{x})] = \bar{f}(\bm{x}) – E_{\mathcal{D}}[\hat{f}(\bm{x})] = \bar{f}(\bm{x}) – \bar{f}(\bm{x}) = 0 $$

したがって交差項は消え、

$$ E_{\mathcal{D}}\left[(f(\bm{x}) – \hat{f}(\bm{x}))^2\right] = \underbrace{(f(\bm{x}) – \bar{f}(\bm{x}))^2}_{\text{Bias}^2} + \underbrace{E_{\mathcal{D}}\left[(\hat{f}(\bm{x}) – \bar{f}(\bm{x}))^2\right]}_{\text{Variance}} $$

バイアス-バリアンス分解の完成

以上をまとめると、

$$ \boxed{\text{EPE}(\bm{x}) = \text{Bias}^2[\hat{f}(\bm{x})] + \text{Var}[\hat{f}(\bm{x})] + \sigma^2} $$

各項の意味は以下の通りです。

定義 意味
$\text{Bias}^2$ $(f(\bm{x}) – E[\hat{f}(\bm{x})])^2$ 平均的な予測と真の値のずれ
$\text{Var}$ $E[(\hat{f}(\bm{x}) – E[\hat{f}(\bm{x})])^2]$ 訓練データの変動による予測の揺らぎ
$\sigma^2$ $\text{Var}[\epsilon]$ データ固有のノイズ(削減不可能)

バギングの目標は、Bias をほぼ変えずに Variance を削減することです。

バギングのアルゴリズム

ブートストラップサンプリング

元の訓練データ $\mathcal{D} = \{(\bm{x}_i, y_i)\}_{i=1}^n$ から、復元抽出で $n$ 個のサンプルを選びます。この操作をブートストラップサンプリングと呼び、得られたデータセットをブートストラップサンプル $\mathcal{D}_b^*$ と呼びます。

$\mathcal{D}_b^*$ には同じサンプルが複数回含まれることがあり、逆に一度も選ばれないサンプルも存在します。

バギングの手順

  1. $b = 1, 2, \dots, B$ について: – $\mathcal{D}$ からブートストラップサンプル $\mathcal{D}_b^*$ を生成 – $\mathcal{D}_b^*$ でモデル $\hat{f}_b$ を学習
  2. 予測を平均化: – 回帰: $\hat{f}_{\text{bag}}(\bm{x}) = \frac{1}{B}\sum_{b=1}^{B} \hat{f}_b(\bm{x})$ – 分類: $\hat{y}_{\text{bag}}(\bm{x}) = \text{majority vote}\left(\hat{f}_1(\bm{x}), \dots, \hat{f}_B(\bm{x})\right)$

バギングによる分散削減の証明

独立同分布の場合

まず理想的な状況として、$B$ 個のモデル $\hat{f}_1, \dots, \hat{f}_B$ が互いに独立で同一の分布に従う場合を考えます。

各モデルの期待値と分散は共通で、

$$ E[\hat{f}_b(\bm{x})] = \mu, \quad \text{Var}[\hat{f}_b(\bm{x})] = \sigma_f^2 \quad (b = 1, \dots, B) $$

バギング予測の平均は $\hat{f}_{\text{bag}}(\bm{x}) = \frac{1}{B}\sum_{b=1}^{B} \hat{f}_b(\bm{x})$ です。

期待値(バイアスは変化しない):

$$ \begin{align} E[\hat{f}_{\text{bag}}(\bm{x})] &= E\left[\frac{1}{B}\sum_{b=1}^{B} \hat{f}_b(\bm{x})\right] \\ &= \frac{1}{B}\sum_{b=1}^{B} E[\hat{f}_b(\bm{x})] \\ &= \frac{1}{B} \cdot B \cdot \mu \\ &= \mu \end{align} $$

バギング予測の期待値は個々のモデルの期待値と同じです。したがって、バイアスは変化しません。

分散(独立の場合 $1/B$ に削減):

$$ \begin{align} \text{Var}[\hat{f}_{\text{bag}}(\bm{x})] &= \text{Var}\left[\frac{1}{B}\sum_{b=1}^{B} \hat{f}_b(\bm{x})\right] \\ &= \frac{1}{B^2} \text{Var}\left[\sum_{b=1}^{B} \hat{f}_b(\bm{x})\right] \\ &= \frac{1}{B^2} \sum_{b=1}^{B} \text{Var}[\hat{f}_b(\bm{x})] \quad (\because \text{独立なので共分散が0}) \\ &= \frac{1}{B^2} \cdot B \cdot \sigma_f^2 \\ &= \frac{\sigma_f^2}{B} \end{align} $$

$B$ 個の独立なモデルの平均は、分散を $1/B$ に削減します。$B \to \infty$ で分散は0に収束します。

相関がある場合(現実のバギング)

しかし現実のバギングでは、各モデルは同じ訓練データ $\mathcal{D}$ からブートストラップサンプリングで得た $\mathcal{D}_b^*$ で学習するため、$\hat{f}_b$ 間に相関が存在します。

任意の2つのモデル $\hat{f}_b, \hat{f}_{b’}$ ($b \neq b’$) の相関係数を $\rho$ とします。

$$ \rho = \text{Corr}[\hat{f}_b(\bm{x}), \hat{f}_{b’}(\bm{x})] = \frac{\text{Cov}[\hat{f}_b(\bm{x}), \hat{f}_{b’}(\bm{x})]}{\sigma_f^2} $$

つまり $\text{Cov}[\hat{f}_b(\bm{x}), \hat{f}_{b’}(\bm{x})] = \rho \sigma_f^2$ です。

バギング予測の分散を計算します。

$$ \text{Var}[\hat{f}_{\text{bag}}(\bm{x})] = \frac{1}{B^2} \text{Var}\left[\sum_{b=1}^{B} \hat{f}_b(\bm{x})\right] $$

和の分散を展開します。

$$ \text{Var}\left[\sum_{b=1}^{B} \hat{f}_b\right] = \sum_{b=1}^{B} \text{Var}[\hat{f}_b] + \sum_{b=1}^{B} \sum_{b’ \neq b} \text{Cov}[\hat{f}_b, \hat{f}_{b’}] $$

第1項は $B \sigma_f^2$ です。第2項は $B(B-1)$ 個の共分散の和なので、

$$ \sum_{b=1}^{B} \sum_{b’ \neq b} \text{Cov}[\hat{f}_b, \hat{f}_{b’}] = B(B-1) \rho \sigma_f^2 $$

合計して $B^2$ で割ります。

$$ \begin{align} \text{Var}[\hat{f}_{\text{bag}}(\bm{x})] &= \frac{1}{B^2} \left(B \sigma_f^2 + B(B-1) \rho \sigma_f^2\right) \\ &= \frac{\sigma_f^2}{B^2} \left(B + B(B-1)\rho\right) \\ &= \frac{\sigma_f^2}{B^2} \cdot B \left(1 + (B-1)\rho\right) \\ &= \frac{\sigma_f^2}{B} \left(1 + (B-1)\rho\right) \end{align} $$

整理すると、

$$ \begin{align} \text{Var}[\hat{f}_{\text{bag}}(\bm{x})] &= \frac{\sigma_f^2}{B} + \frac{(B-1)\rho \sigma_f^2}{B} \\ &= \frac{\sigma_f^2 (1 – \rho)}{B} + \rho \sigma_f^2 \end{align} $$

最後の変形を確認します。

$$ \frac{\sigma_f^2}{B} + \frac{(B-1)\rho \sigma_f^2}{B} = \frac{\sigma_f^2 + (B-1)\rho \sigma_f^2}{B} = \frac{\sigma_f^2(1 + B\rho – \rho)}{B} = \frac{\sigma_f^2(1-\rho)}{B} + \rho\sigma_f^2 $$

$$ \boxed{\text{Var}[\hat{f}_{\text{bag}}(\bm{x})] = \rho \sigma_f^2 + \frac{(1-\rho)\sigma_f^2}{B}} $$

この公式から重要な知見が得られます。

  1. 第2項 $\frac{(1-\rho)\sigma_f^2}{B}$: $B$ を増やすことで0に近づけられる項です。
  2. 第1項 $\rho \sigma_f^2$: $B$ を増やしても削減できない項です。モデル間の相関 $\rho$ に比例します。
  3. $\rho = 0$(完全に独立)のとき、$\text{Var} = \sigma_f^2 / B$ となり、理想的な $1/B$ の分散削減が実現します。
  4. $\rho = 1$(完全に相関)のとき、$\text{Var} = \sigma_f^2$ となり、バギングの効果はありません。
  5. したがって、モデル間の相関 $\rho$ を下げることが分散削減の鍵 です。ランダムフォレストが特徴量のランダム選択を行うのは、まさにこの $\rho$ を下げるためです。

OOB(Out-of-Bag)推定の理論

OOBサンプルとは

ブートストラップサンプル $\mathcal{D}_b^*$ は復元抽出で $n$ 個を選ぶため、元のデータの一部は $\mathcal{D}_b^*$ に含まれません。$\mathcal{D}_b^*$ に含まれなかったサンプルをOOB(Out-of-Bag)サンプルと呼びます。

各サンプルがOOBになる確率

元データの特定のサンプル $(\bm{x}_i, y_i)$ が1回の抽出で選ばれない確率は $1 – 1/n$ です。$n$ 回の復元抽出(ブートストラップ1回分)で1度も選ばれない確率は、

$$ P(\text{OOB}) = \left(1 – \frac{1}{n}\right)^n $$

$n \to \infty$ での極限を求めます。対数を取ると、

$$ \ln\left(1 – \frac{1}{n}\right)^n = n \ln\left(1 – \frac{1}{n}\right) $$

$\ln(1-x) \approx -x – x^2/2 – \cdots$ のテイラー展開を $x = 1/n$ に適用すると、

$$ \begin{align} n \ln\left(1 – \frac{1}{n}\right) &= n \left(-\frac{1}{n} – \frac{1}{2n^2} – \frac{1}{3n^3} – \cdots\right) \\ &= -1 – \frac{1}{2n} – \frac{1}{3n^2} – \cdots \\ &\to -1 \quad (n \to \infty) \end{align} $$

したがって、

$$ \boxed{\left(1 – \frac{1}{n}\right)^n \to e^{-1} \approx 0.368} $$

つまり、各ブートストラップサンプルにおいて、元データの約36.8%のサンプルがOOBとなります。逆に言えば、各ブートストラップサンプルには元データの約63.2%のユニークなサンプルが含まれます。

OOB誤差推定

OOBサンプルの性質を利用すると、交差検証を行わずにモデルの汎化性能を推定できます。

各サンプル $(\bm{x}_i, y_i)$ に対して:

  1. $(\bm{x}_i, y_i)$ がOOBであるブートストラップサンプルの番号を $\mathcal{B}_i = \{b : (\bm{x}_i, y_i) \notin \mathcal{D}_b^*\}$ とする
  2. OOBモデルによる予測を平均化: $\hat{f}_{\text{OOB}}(\bm{x}_i) = \frac{1}{|\mathcal{B}_i|}\sum_{b \in \mathcal{B}_i} \hat{f}_b(\bm{x}_i)$
  3. OOB誤差: $\text{OOB Error} = \frac{1}{n}\sum_{i=1}^{n} L(y_i, \hat{f}_{\text{OOB}}(\bm{x}_i))$

各 $\hat{f}_b$ は $(\bm{x}_i, y_i)$ を学習に使っていないため、この推定は交差検証と同等のバイアスの少ない汎化誤差推定になります。

OOB推定の計算量は $O(nB)$(各サンプルについてOOBモデルの予測を集計するだけ)であり、$K$-fold交差検証の $O(KnB/K) = O(nB)$ と同程度ですが、追加の学習が不要という利点があります。

サブサンプリングとの比較

ブートストラップ(復元抽出)の代わりに、非復元抽出(サブサンプリング)を使うこともできます。$n$ 個から $m$ 個を非復元で選ぶ方法です。

ブートストラップ サブサンプリング
抽出方法 復元抽出($n$個選択) 非復元抽出($m$個選択)
各サブセットのサイズ $n$(ただしユニーク約$0.632n$) $m$(全てユニーク)
OOBサンプル 約$0.368n$個 $n – m$個
モデル間の相関 比較的高い $m < n$ で低くなる

サブサンプリングで $m$ を小さくすると各モデルの独立性が増し、相関 $\rho$ が下がります。しかし、各モデルが使うデータ量が減るためバイアスが増加します。ブートストラップは各サブセットのサイズを $n$ に保ちながらランダム性を導入するバランスの取れた方法です。

Pythonでの実装

決定木のバギング(スクラッチ実装)

scikit-learnのDecisionTreeRegressorを基本学習器として使い、バギングのアルゴリズムをスクラッチ実装します。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.tree import DecisionTreeRegressor
from sklearn.metrics import mean_squared_error

# ----------------------------
# バギング回帰器(スクラッチ実装)
# ----------------------------
class BaggingRegressor:
    """決定木のバギング(Bootstrap Aggregating)"""

    def __init__(self, n_estimators=50, max_depth=None, random_state=None):
        self.n_estimators = n_estimators
        self.max_depth = max_depth
        self.random_state = random_state
        self.estimators = []
        self.oob_indices = []  # 各推定器のOOBサンプルのインデックス

    def fit(self, X, y):
        """バギングの学習"""
        n_samples = X.shape[0]
        rng = np.random.RandomState(self.random_state)

        self.estimators = []
        self.oob_indices = []

        for b in range(self.n_estimators):
            # ブートストラップサンプリング(復元抽出)
            bootstrap_idx = rng.choice(n_samples, size=n_samples, replace=True)

            # OOBインデックスの計算
            oob_idx = np.setdiff1d(np.arange(n_samples), np.unique(bootstrap_idx))
            self.oob_indices.append(oob_idx)

            # 決定木の学習
            tree = DecisionTreeRegressor(max_depth=self.max_depth,
                                         random_state=rng.randint(0, 10000))
            tree.fit(X[bootstrap_idx], y[bootstrap_idx])
            self.estimators.append(tree)

        return self

    def predict(self, X):
        """予測(全推定器の平均)"""
        predictions = np.array([tree.predict(X) for tree in self.estimators])
        return np.mean(predictions, axis=0)

    def predict_individual(self, X):
        """各推定器の個別予測を返す"""
        return np.array([tree.predict(X) for tree in self.estimators])

    def oob_score(self, X, y):
        """OOBスコア(R^2)の計算"""
        n_samples = X.shape[0]
        oob_preds = np.zeros(n_samples)
        oob_counts = np.zeros(n_samples)

        for b, tree in enumerate(self.estimators):
            oob_idx = self.oob_indices[b]
            if len(oob_idx) > 0:
                oob_preds[oob_idx] += tree.predict(X[oob_idx])
                oob_counts[oob_idx] += 1

        # OOB予測が得られたサンプルのみで評価
        valid = oob_counts > 0
        oob_preds[valid] /= oob_counts[valid]

        ss_res = np.sum((y[valid] - oob_preds[valid]) ** 2)
        ss_tot = np.sum((y[valid] - np.mean(y[valid])) ** 2)

        return 1 - ss_res / ss_tot

# ----------------------------
# データ生成: 非線形回帰タスク
# ----------------------------
np.random.seed(42)

# 真の関数
def true_function(x):
    return np.sin(2 * x) + 0.5 * np.cos(5 * x)

n_train = 100
n_test = 200
noise_std = 0.3

# 訓練データ
X_train = np.random.uniform(-3, 3, (n_train, 1))
y_train = true_function(X_train[:, 0]) + noise_std * np.random.randn(n_train)

# テストデータ
X_test = np.linspace(-3, 3, n_test).reshape(-1, 1)
y_test = true_function(X_test[:, 0])

# ----------------------------
# 実験1: バギングによる分散削減の検証
# ----------------------------
n_experiments = 50  # 異なる訓練データで実験を繰り返す
B_values = [1, 5, 10, 25, 50]

# 各推定器数での予測の分散を計算
fig, axes = plt.subplots(2, 3, figsize=(18, 10))
axes = axes.flatten()

for idx, B in enumerate(B_values):
    all_predictions = np.zeros((n_experiments, n_test))

    for exp in range(n_experiments):
        # 毎回異なる訓練データを生成
        X_exp = np.random.uniform(-3, 3, (n_train, 1))
        y_exp = true_function(X_exp[:, 0]) + noise_std * np.random.randn(n_train)

        bag = BaggingRegressor(n_estimators=B, max_depth=None, random_state=exp)
        bag.fit(X_exp, y_exp)
        all_predictions[exp] = bag.predict(X_test)

    # 平均予測と分散
    mean_pred = np.mean(all_predictions, axis=0)
    var_pred = np.var(all_predictions, axis=0)

    ax = axes[idx]
    ax.plot(X_test[:, 0], y_test, 'k-', linewidth=2, label='True function')
    ax.plot(X_test[:, 0], mean_pred, 'r-', linewidth=2, label='Mean prediction')
    ax.fill_between(X_test[:, 0],
                     mean_pred - 2 * np.sqrt(var_pred),
                     mean_pred + 2 * np.sqrt(var_pred),
                     alpha=0.3, color='red', label='$\pm 2\sigma$')
    ax.set_title(f'B = {B}, Mean Var = {np.mean(var_pred):.4f}')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)

# 分散の推移
ax = axes[5]
mean_vars = []
for B in range(1, 51):
    all_preds = np.zeros((n_experiments, n_test))
    for exp in range(n_experiments):
        X_exp = np.random.uniform(-3, 3, (n_train, 1))
        y_exp = true_function(X_exp[:, 0]) + noise_std * np.random.randn(n_train)
        bag = BaggingRegressor(n_estimators=B, max_depth=None, random_state=exp)
        bag.fit(X_exp, y_exp)
        all_preds[exp] = bag.predict(X_test)
    mean_vars.append(np.mean(np.var(all_preds, axis=0)))

ax.plot(range(1, 51), mean_vars, 'b-', linewidth=2)
ax.set_xlabel('Number of Estimators (B)')
ax.set_ylabel('Mean Variance')
ax.set_title('Variance vs Number of Estimators')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# ----------------------------
# 実験2: OOBスコアの検証
# ----------------------------
print("=== OOBスコアとテストスコアの比較 ===")
bag_full = BaggingRegressor(n_estimators=100, max_depth=None, random_state=42)
bag_full.fit(X_train, y_train)

# OOBスコア
oob = bag_full.oob_score(X_train, y_train)

# テストスコア
y_pred_test = bag_full.predict(X_test)
ss_res = np.sum((y_test - y_pred_test) ** 2)
ss_tot = np.sum((y_test - np.mean(y_test)) ** 2)
test_r2 = 1 - ss_res / ss_tot

print(f"OOBスコア (R^2): {oob:.4f}")
print(f"テストスコア (R^2): {test_r2:.4f}")

# ----------------------------
# 実験3: OOB確率の検証
# ----------------------------
print("\n=== OOB確率の検証 ===")
n_samples = 1000
n_bootstrap = 10000
oob_counts = np.zeros(n_samples)

rng = np.random.RandomState(42)
for _ in range(n_bootstrap):
    bootstrap_idx = rng.choice(n_samples, size=n_samples, replace=True)
    unique_idx = np.unique(bootstrap_idx)
    oob_idx = np.setdiff1d(np.arange(n_samples), unique_idx)
    oob_counts[oob_idx] += 1

empirical_oob_prob = np.mean(oob_counts) / n_bootstrap
theoretical_oob_prob = (1 - 1/n_samples) ** n_samples

print(f"経験的OOB確率: {empirical_oob_prob:.6f}")
print(f"理論的OOB確率 (1-1/n)^n: {theoretical_oob_prob:.6f}")
print(f"理論的極限値 1/e: {1/np.e:.6f}")

# ----------------------------
# 実験4: 単一決定木 vs バギングの比較
# ----------------------------
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 単一決定木
single_tree = DecisionTreeRegressor(max_depth=None, random_state=42)
single_tree.fit(X_train, y_train)
y_pred_tree = single_tree.predict(X_test)

# バギング(100本)
y_pred_bag = bag_full.predict(X_test)

# 単一決定木のプロット
ax1.scatter(X_train[:, 0], y_train, alpha=0.3, s=20, label='Training data')
ax1.plot(X_test[:, 0], y_test, 'k-', linewidth=2, label='True function')
ax1.plot(X_test[:, 0], y_pred_tree, 'r-', linewidth=1.5, label='Single tree')
mse_tree = mean_squared_error(y_test, y_pred_tree)
ax1.set_title(f'Single Decision Tree (MSE = {mse_tree:.4f})')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.legend()
ax1.grid(True, alpha=0.3)

# バギングのプロット
ax2.scatter(X_train[:, 0], y_train, alpha=0.3, s=20, label='Training data')
ax2.plot(X_test[:, 0], y_test, 'k-', linewidth=2, label='True function')
ax2.plot(X_test[:, 0], y_pred_bag, 'r-', linewidth=1.5, label='Bagging (B=100)')
mse_bag = mean_squared_error(y_test, y_pred_bag)
ax2.set_title(f'Bagging with B=100 (MSE = {mse_bag:.4f})')
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.legend()
ax2.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"\n単一決定木のMSE: {mse_tree:.4f}")
print(f"バギング(B=100)のMSE: {mse_bag:.4f}")
print(f"MSE削減率: {(1 - mse_bag/mse_tree) * 100:.1f}%")

上記のコードを実行すると、以下のことが確認できます。

  1. 分散削減の実験的検証: $B$ を増やすにつれて予測の分散($\pm 2\sigma$ の幅)が減少し、理論通りの結果が得られます。
  2. OOBスコアの信頼性: OOBスコアとテストスコアが近い値になり、OOB推定がモデル評価として妥当であることが確認できます。
  3. OOB確率の検証: 経験的なOOB確率が理論値 $1/e \approx 0.368$ に一致することが確認できます。
  4. 単一決定木 vs バギング: バギングにより予測が滑らかになり、MSEが大幅に削減されます。

まとめ

本記事では、バギング(Bootstrap Aggregating)の理論について解説しました。

  • バイアス-バリアンス分解により、期待予測誤差は $\text{Bias}^2 + \text{Var} + \sigma^2$ に分解されます
  • バギングはバイアスをほぼ変えずに分散を削減するアンサンブル手法です
  • 独立なモデルの平均は分散を $\sigma_f^2 / B$ に削減しますが、現実には相関 $\rho$ があるため $\rho \sigma_f^2 + (1-\rho)\sigma_f^2 / B$ が下限となります
  • 相関 $\rho$ を下げることが分散削減の鍵であり、これがランダムフォレストの動機になります
  • OOB(Out-of-Bag)推定により、交差検証なしで汎化誤差を推定できます。各サンプルがOOBになる確率は $(1-1/n)^n \to 1/e \approx 0.368$ です

次のステップとして、以下の記事も参考にしてください。