Pythonでのブートストラップ法の実装

ブートストラップ法は、限られたデータから母集団の統計量を推定する際に強力な手法です。特に、標本サイズが小さい場合や、母集団の分布が不明な場合に有効です。この記事では、ブートストラップ法の基本的な考え方と、Pythonを用いた実装方法を解説します。

ブートストラップ法とは
- ブートストラップ法の利点
- ブートストラップ法の具体例
Pythonによる実装
まとめ

ブートストラップ法とは

ブートストラップ法は、手元にある標本データ（サンプル）を母集団とみなして、そこから復元抽出（重複を許してランダムに抽出）を繰り返すことで、多数の仮想的な標本データセットを生成する手法です。各データセットから統計量（例えば平均値や中央値など）を計算し、その分布を調べることで、元の標本の統計量のばらつきや信頼区間を推定します。

復元抽出とは

復元抽出は、毎回同じ条件でサンプリングを行うために、抽出した要素を元に戻してから次の要素を抽出する方法です。これにより、同じ要素が複数回選ばれる可能性があります。

ブートストラップ法の利点

母集団の分布を仮定しない: 正規分布などの特定の分布を仮定せずに、標本データから直接統計量の分布を推定できます。
実装が容易: 比較的シンプルなアルゴリズムで実装できます。
様々な統計量に適用可能: 平均値、中央値、標準偏差など、様々な統計量の推定に利用できます。

ブートストラップ法の具体例

ある湖に生息する魚の平均体重を知りたいとします。しかし、湖のすべての魚を捕まえて体重を測るのは現実的ではありません。そこで、100匹の魚をランダムに捕まえ（標本抽出）、その体重を測定します。

この100匹の魚の体重データだけを使って、湖全体の魚の平均体重（母集団の平均体重）を推定したい場合に、ブートストラップ法が役立ちます。具体的には以下のような考え方をします。

標本を母集団とみなす: 100匹の魚の体重データを、湖全体の魚の体重分布を代表するものとみなします。
復元抽出: 100匹のデータから、重複を許してランダムに100匹分の体重データを選びます（復元抽出）。これにより、例えば、3番目に捕まえた魚の体重データが2回選ばれたり、50番目の魚の体重データが一度も選ばれない、といったことが起こりえます。
仮想標本の作成: 2.の復元抽出を何度も繰り返します（例えば1000回）。すると、100匹の魚の体重データからなる仮想的な標本データセットが1000個できます。
統計量の計算: 1000個の仮想標本それぞれについて、平均体重を計算します。
分布の確認と推定: 1000個の平均体重の分布を調べることで、元の標本（最初に捕まえた100匹）の平均体重のばらつき具合（どれくらい信頼できるか）や、95%信頼区間（母集団の平均体重が、95%の確率で含まれると推定される範囲）などを知ることができます。

このように、ブートストラップ法は、限られたデータから、あたかも大量のデータがあるかのようにシミュレーションを行うことで、母集団の統計量を推定する強力な手法です。

Pythonによる実装

ここでは、Pythonのnumpyとscipyライブラリを使ってブートストラップ法を実装します。

1. データの準備

例として、scikit-learnに含まれるirisデータセット（アヤメの計測データ）を使用します。

from sklearn.datasets import load_iris
import numpy as np

iris = load_iris()
data = iris.data[:, 0]  # 0列目（がく片の長さ）を使用

ここでは、iris.data[:, 0]により、irisデータセットの4つの特徴量のうち、1列目(インデックスは0)の「がく片の長さ(sepal length)」のデータだけを抽出して、data変数に格納しています。これにより分析対象を「がく片の長さ」に限定します。

2. ブートストラップサンプルの生成

numpy.random.choice()関数を使って、元のデータから復元抽出を行い、ブートストラップサンプルを生成します。

def bootstrap_sample(data, sample_size):
  """
  ブートストラップサンプルを生成する関数

  Args:
    data: 元のデータ (numpy配列)
    sample_size: 生成するサンプルのサイズ

  Returns:
    ブートストラップサンプル (numpy配列)
  """
  return np.random.choice(data, size=sample_size, replace=True)

元のデータと同じサイズの仮想的な標本を生成します。np.random.choice(data, size=sample_size, replace=True)はdataからsample_size個の要素を復元抽出(replace=True)するという処理です。

3. 統計量の計算

各ブートストラップサンプルから、求めたい統計量を計算します。ここでは、平均値を計算する関数を定義します。

def calculate_statistic(sample):
  """
  サンプルから統計量（ここでは平均値）を計算する関数

  Args:
    sample: データサンプル (numpy配列)

  Returns:
    統計量 (float)
  """
  return np.mean(sample)

入力されたsampleの平均値を計算します。np.mean(sample)により、ブートストラップサンプルごとに平均を計算します。

4. ブートストラップ法の実行

ここが、ブートストラップ法の中核となる処理です。ブートストラップ標本から統計量（ここでは平均値）を繰り返し計算し、その結果をリストに格納します。

n_iterations = 1000  # ブートストラップサンプルの生成回数
sample_size = len(data)  # サンプルサイズは元のデータと同じにする

# ブートストラップ統計量を格納するリスト
bootstrap_statistics = []

for _ in range(n_iterations):
  sample = bootstrap_sample(data, sample_size)
  statistic = calculate_statistic(sample)
  bootstrap_statistics.append(statistic)

反復回数とサンプルサイズの設定

n_iterations = 1000: ブートストラップサンプルを生成する回数を1000回に設定します。この回数が多いほど、推定の精度が向上する傾向にありますが、計算時間も長くなります。
sample_size = len(data): ブートストラップサンプル1つあたりのサイズを、元のデータ(data)のサイズと同じにします。つまり、元のデータが150個の要素を持っていれば、ブートストラップサンプルも150個の要素を持つようにします。

ブートストラップサンプルの生成と統計量の計算の繰り返し

ブートストラップサンプルの生成と統計量の計算を1000回繰り返します。bootstrap_sample関数（前に定義）を呼び出し、元のデータ data から復元抽出によって、sample_size個の要素を持つブートストラップサンプル sample を生成します。続けて、calculate_statistic関数（前に定義）を呼び出し、生成されたブートストラップサンプル sample の統計量（ここでは平均値）を計算し、statistic に格納します。

この処理により、bootstrap_statistics リストには、1000個のブートストラップサンプルそれぞれの平均値が格納されます。この1000個の平均値の分布を調べることで、元の標本の平均値のばらつきや信頼区間を推定することができます（これは次のステップの「結果の可視化と分析」で行ういます）。

5. 結果の可視化と分析

ブートストラップ統計量の分布をヒストグラムで可視化し、信頼区間などを計算します。

import matplotlib.pyplot as plt

# ヒストグラムの表示
plt.hist(bootstrap_statistics, bins=30)
plt.xlabel('Bootstrap Mean')
plt.ylabel('Frequency')
plt.title('Distribution of Bootstrap Means')
plt.show()

# 95%信頼区間の計算
confidence_interval = np.percentile(bootstrap_statistics, [2.5, 97.5])
print(f"95%信頼区間: {confidence_interval}")

1000個のブートストラップ平均値の分布をヒストグラムで表示します。plt.hist(bootstrap_statistics, bins=30)によって、bins=30でヒストグラムの棒の数を30に設定します。 * np.percentile(bootstrap_statistics, [2.5, 97.5])では bootstrap_statistics（1000個のブートストラップ平均）の2.5パーセンタイル点と97.5パーセンタイル点を計算しています。これにより、95%信頼区間を求めることができます。

応用: パーセンタイルブートストラップ法

上記の例は、最も基本的なブートストラップ法です。母集団の分布が歪んでいる場合など、より正確な信頼区間を求めるために、パーセンタイルブートストラップ法などの改良された方法があります。

パーセンタイルブートストラップ法では、ブートストラップ統計量の分布から直接パーセンタイル点を求めて信頼区間とします。上記のコード例の最後の部分（信頼区間の計算）は、すでにパーセンタイルブートストラップ法の実装になっています。

# (上記のコードの続き) 95%信頼区間を計算
confidence_interval = np.percentile(bootstrap_statistics, [2.5, 97.5])
print("95% confidence interval:", confidence_interval)

この95%信頼区間とは、同じ方法で推定を100回繰り返すと、そのうち95回はこの区間内に母集団の平均値が含まれることを意味します。

応用: 他の信頼区間の計算方法

scipy.stats.bootstrapを使うとより簡単に信頼区間を計算できます。

from scipy import stats
# axis=0を指定しないとdata全体で一つの統計量を計算してしまうので注意
res = stats.bootstrap((data,), np.mean, confidence_level=0.95,
                        random_state=np.random.default_rng())

print(res.confidence_interval)

stats.bootstrap((data,), np.mean, confidence_level=0.95, random_state=np.random.default_rng())の各値は下記のような意味があります