PythonとXGBoostによる睡眠品質分類器の最適化：特徴量エンジニアリングガイド

pip install pandas numpy scikit-learn xgboost scipy

# src/data_simulation.py
import pandas as pd
import numpy as np
from datetime import datetime, timedelta

def simulate_wearable_data(nights=50):
    """
    複数夜の合成ウェアラブルデータセットを生成します。
    
    各夜は1秒間隔で8時間分のデータで構成されます。
    特徴量：heart_rate（bpm）、acceleration（g）。
    """
    data_frames = []
    
    for night in range(nights):
        # この夜のランダムな睡眠品質を割り当て
        quality = np.random.choice(['Poor', 'Average', 'Good'], p=[0.25, 0.45, 0.3])
        start_time = datetime(2025, 1, 1, 22, 0, 0) + timedelta(days=night)
        
        # 睡眠品質に基づく基本パラメータ
        if quality == 'Good':
            hr_base, hr_var = 60, 3
            acc_base, acc_var = 0.01, 0.005
            num_awakenings = np.random.randint(0, 2)
        elif quality == 'Average':
            hr_base, hr_var = 68, 5
            acc_base, acc_var = 0.02, 0.01
            num_awakenings = np.random.randint(2, 5)
        else: # Poor
            hr_base, hr_var = 75, 8
            acc_base, acc_var = 0.04, 0.02
            num_awakenings = np.random.randint(5, 10)
            
        # 8時間分のデータを生成（28800秒）
        timestamps = pd.to_datetime([start_time + timedelta(seconds=i) for i in range(28800)])
        
        # データをシミュレート
        hr = hr_base + np.random.randn(28800) * hr_var
        acc = np.abs(acc_base + np.random.randn(28800) * acc_var)
        
        # 覚醒や不安定さのスパイクを追加
        for _ in range(num_awakenings):
            idx = np.random.randint(0, 28800)
            spike_duration = np.random.randint(60, 300)
            hr[idx:idx+spike_duration] += np.random.uniform(5, 15)
            acc[idx:idx+spike_duration] += np.random.uniform(0.1, 0.5)

        night_df = pd.DataFrame({
            'timestamp': timestamps,
            'heart_rate': hr,
            'acceleration': acc,
            'night_id': night,
            'sleep_quality': quality
        })
        data_frames.append(night_df)
        
    return pd.concat(data_frames, ignore_index=True)

# データを生成して保存
raw_data = simulate_wearable_data(nights=100)
raw_data.to_csv('wearable_sleep_data.csv', index=False)

print("Generated wearable_sleep_data.csv with 100 nights of data.")
print(raw_data.head())

# src/feature_engineering.py
import pandas as pd
import numpy as np
from scipy.stats import iqr

def calculate_hrv(hr_series):
    """心拍数シリーズから心拍変動（RMSSD）を計算します。"""
    # 心拍間隔（ミリ秒単位での心拍間の時間）を計算
    rr_intervals = 60000 / hr_series
    # 連続する差分を計算
    successive_diffs = np.diff(rr_intervals)
    # RMSSD（連続する差分の二乗平均平方根）を計算
    if len(successive_diffs) > 0:
        rmssd = np.sqrt(np.mean(successive_diffs ** 2))
    else:
        rmssd = 0
    return rmssd

def engineer_features(df):
    """
    夜のグループ化された生のウェアラブルデータから特徴量をエンジニアリングします。
    """
    # 夜でグループ化
    grouped = df.groupby('night_id')
    
    feature_list = []
    for night_id, group in grouped:
        features = {'night_id': night_id}
        
        # 基本的な心拍数特徴量
        features['hr_mean'] = group['heart_rate'].mean()
        features['hr_std'] = group['heart_rate'].std()
        features['hr_min'] = group['heart_rate'].min()
        features['hr_max'] = group['heart_rate'].max()
        
        # 心拍変動（HRV）
        features['hrv_rmssd'] = calculate_hrv(group['heart_rate'])
        
        # 加速度特徴量（動き）
        features['acc_mean'] = group['acceleration'].mean()
        features['acc_std'] = group['acceleration'].std()
        features['acc_max'] = group['acceleration'].max()
        features['restless_moments'] = (group['acceleration'] > 0.1).sum() # 高動きをした瞬間をカウント
        
        # 睡眠段階の持続時間（簡略版）
        # 深い睡眠：低心拍、低動き
        deep_sleep_mask = (group['heart_rate'] < (features['hr_mean'] * 0.9)) & (group['acceleration'] < 0.05)
        features['deep_sleep_duration_pct'] = deep_sleep_mask.sum() / len(group)
        
        # 浅い睡眠：中程度のすべて
        light_sleep_mask = ~deep_sleep_mask
        features['light_sleep_duration_pct'] = light_sleep_mask.sum() / len(group)
        
        # ターゲット変数
        features['sleep_quality'] = group['sleep_quality'].iloc[0]
        
        feature_list.append(features)
        
    return pd.DataFrame(feature_list)

# 生データを読み込んで特徴量をエンジニアリング
raw_data = pd.read_csv('wearable_sleep_data.csv')
featured_data = engineer_features(raw_data)
featured_data.to_csv('featured_sleep_data.csv', index=False)

print("Engineered features and saved to featured_sleep_data.csv")
print(featured_data.head())

# src/train_model.py
import pandas as pd
import xgboost as xgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report, accuracy_score
from sklearn.preprocessing import LabelEncoder

# 特徴量データを読み込み
df = pd.read_csv('featured_sleep_data.csv')

# XGBoost用にデータを準備
X = df.drop(['night_id', 'sleep_quality'], axis=1)
y = df['sleep_quality']

# ターゲット変数をエンコード（例：'Poor' -> 0、'Average' -> 1、'Good' -> 2）
le = LabelEncoder()
y_encoded = le.fit_transform(y)

# データを訓練セットとテストセットに分割
X_train, X_test, y_train, y_test = train_test_split(
    X, y_encoded, test_size=0.2, random_state=42, stratify=y_encoded
)

# XGBoost分類器を初期化して訓練
model = xgb.XGBClassifier(
    objective='multi:softmax', # マルチクラス分類用
    num_class=len(le.classes_),
    eval_metric='mlogloss',
    use_label_encoder=False,
    seed=42
)

model.fit(X_train, y_train)

# 予測を行う
y_pred = model.predict(X_test)

# モデルを評価
accuracy = accuracy_score(y_test, y_pred)
print(f"Model Accuracy: {accuracy * 100:.2f}%")

# 詳細な分類レポートを印刷
print("\
Classification Report:")
print(classification_report(y_test, y_pred, target_names=le.classes_))

# train_model.pyの末尾に追加
import matplotlib.pyplot as plt

# 特徴量の重要度をプロット
feature_importances = pd.DataFrame(
    {'feature': X.columns, 'importance': model.feature_importances_}
).sort_values('importance', ascending=False)

print("\
Feature Importances:")
print(feature_importances)

# オプション：プロット
# plt.figure(figsize=(10, 6))
# plt.barh(feature_importances['feature'], feature_importances['importance'])
# plt.xlabel("XGBoost Feature Importance")
# plt.gca().invert_yaxis()
# plt.show()