Claude Code × データサイエンス — pandas・scikit-learn・機械学習プロジェクトをAIペアプログラミングで加速する

# src/feature_engineering.py
import pandas as pd
import numpy as np
from sklearn.preprocessing import LabelEncoder
 
RANDOM_STATE = 42
 
def engineer_features(df: pd.DataFrame) -> pd.DataFrame:
    """
    特徴量エンジニアリングのメインパイプライン。
    元データを変更しないよう .copy() してから処理する。
    """
    df = df.copy()
 
    # ① 日付特徴量（未来の情報を使わないためデータリークなし）
    df["purchase_date"] = pd.to_datetime(df["purchase_date"])
    df["day_of_week"] = df["purchase_date"].dt.dayofweek
    df["month"] = df["purchase_date"].dt.month
    df["quarter"] = df["purchase_date"].dt.quarter
    df["is_weekend"] = (df["day_of_week"] >= 5).astype(int)
    df["days_since_epoch"] = (
        df["purchase_date"] - pd.Timestamp("2020-01-01")
    ).dt.days
 
    # ② 価格・数量の組み合わせ
    df["total_amount"] = df["price"] * df["quantity"]
    df["log_price"] = np.log1p(df["price"])
    df["log_quantity"] = np.log1p(df["quantity"])
 
    # ③ 顧客年齢のビニング
    df["age_group"] = pd.cut(
        df["customer_age"],
        bins=[0, 25, 35, 50, 65, 100],
        labels=["under25", "25-35", "35-50", "50-65", "65+"],
        right=False
    ).astype(str)
 
    # ④ カテゴリのLabel Encoding
    for col in ["product_category", "age_group"]:
        le = LabelEncoder()
        df[f"{col}_encoded"] = le.fit_transform(df[col].astype(str))
 
    df = df.drop(columns=["purchase_date", "product_category", "age_group"])
    return df
 
def validate_no_data_leakage(
    df_train: pd.DataFrame,
    df_test: pd.DataFrame,
    date_col: str = "purchase_date"
) -> bool:
    """
    データリークチェック:
    テストデータの日付がすべてトレーニングデータより後であることを確認する
    """
    if date_col not in df_train.columns:
        print(f"警告: {date_col} が見つかりません。スキップします。")
        return True
 
    train_max = pd.to_datetime(df_train[date_col]).max()
    test_min = pd.to_datetime(df_test[date_col]).min()
 
    if test_min <= train_max:
        print(f"⚠️  データリーク検出: "
              f"train最大日({train_max}) >= test最小日({test_min})")
        return False
 
    print(f"✅ データリークなし: "
          f"train最大={train_max}, test最小={test_min}")
    return True

# src/model_comparison.py
import time
import pandas as pd
import numpy as np
from pathlib import Path
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor
from sklearn.model_selection import cross_val_score, KFold
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
 
try:
    from xgboost import XGBRegressor
    XGBOOST_AVAILABLE = True
except ImportError:
    XGBOOST_AVAILABLE = False
    print("注意: XGBoostがインストールされていません。比較から除外します。")
 
RANDOM_STATE = 42
N_SPLITS = 5
 
def compare_models(
    X: pd.DataFrame,
    y: pd.Series,
    output_dir: str = "reports"
) -> pd.DataFrame:
    """
    複数モデルを5-fold CVで比較し、結果をCSVに保存する
    """
    Path(output_dir).mkdir(parents=True, exist_ok=True)
 
    models = {
        "LinearRegression": Pipeline([
            ("scaler", StandardScaler()),
            ("model", LinearRegression())
        ]),
        "Ridge": Pipeline([
            ("scaler", StandardScaler()),
            ("model", Ridge(alpha=1.0, random_state=RANDOM_STATE))
        ]),
        "RandomForest": RandomForestRegressor(
            n_estimators=100, random_state=RANDOM_STATE, n_jobs=-1
        ),
        "GradientBoosting": GradientBoostingRegressor(
            n_estimators=100, random_state=RANDOM_STATE
        ),
    }
    if XGBOOST_AVAILABLE:
        models["XGBoost"] = XGBRegressor(
            n_estimators=100, random_state=RANDOM_STATE,
            n_jobs=-1, verbosity=0
        )
 
    kf = KFold(n_splits=N_SPLITS, shuffle=True, random_state=RANDOM_STATE)
    results = []
 
    for model_name, model in models.items():
        start = time.time()
 
        # RMSE（neg_mean_squared_errorの符号を反転して平方根）
        neg_mse = cross_val_score(
            model, X, y, cv=kf,
            scoring="neg_mean_squared_error", n_jobs=-1
        )
        rmse_scores = np.sqrt(-neg_mse)
 
        r2_scores = cross_val_score(
            model, X, y, cv=kf, scoring="r2", n_jobs=-1
        )
        elapsed = time.time() - start
 
        results.append({
            "model": model_name,
            "rmse_mean": rmse_scores.mean(),
            "rmse_std": rmse_scores.std(),
            "r2_mean": r2_scores.mean(),
            "r2_std": r2_scores.std(),
            "time_sec": elapsed
        })
        print(
            f"{model_name:20s} | "
            f"RMSE: {rmse_scores.mean():.4f} ± {rmse_scores.std():.4f} | "
            f"R²: {r2_scores.mean():.4f} | {elapsed:.1f}s"
        )
 
    results_df = pd.DataFrame(results).sort_values("rmse_mean")
    results_df.to_csv(f"{output_dir}/model_comparison.csv", index=False)
    print(f"\n[OK] Results saved: {output_dir}/model_comparison.csv")
    return results_df
 
# 期待する出力例（データ次第で変わります）:
# LinearRegression     | RMSE: 12.3456 ± 0.4321 | R²: 0.7234 | 0.2s
# Ridge               | RMSE: 12.1234 ± 0.3456 | R²: 0.7456 | 0.3s
# RandomForest        | RMSE:  8.9876 ± 0.2345 | R²: 0.8901 | 12.4s
# GradientBoosting    | RMSE:  8.4567 ± 0.2123 | R²: 0.9012 | 18.7s
# XGBoost             | RMSE:  8.1234 ± 0.1987 | R²: 0.9123 | 5.6s

# src/hyperparameter_tuning.py
import optuna
import numpy as np
import pandas as pd
import joblib
from pathlib import Path
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.model_selection import cross_val_score, KFold
 
RANDOM_STATE = 42
N_SPLITS = 5
N_TRIALS = 50
MODELS_DIR = Path("models")
MODELS_DIR.mkdir(exist_ok=True)
 
def objective(trial: optuna.Trial, X: pd.DataFrame, y: pd.Series) -> float:
    """Optunaの目的関数（最小化: CV-RMSE）"""
    params = {
        "n_estimators": trial.suggest_int("n_estimators", 50, 500),
        "learning_rate": trial.suggest_float("learning_rate", 0.01, 0.3, log=True),
        "max_depth": trial.suggest_int("max_depth", 3, 8),
        "min_samples_split": trial.suggest_int("min_samples_split", 2, 20),
        "min_samples_leaf": trial.suggest_int("min_samples_leaf", 1, 10),
        "subsample": trial.suggest_float("subsample", 0.6, 1.0),
        "max_features": trial.suggest_categorical(
            "max_features", ["sqrt", "log2", None]
        ),
        "random_state": RANDOM_STATE
    }
    model = GradientBoostingRegressor(**params)
    kf = KFold(n_splits=N_SPLITS, shuffle=True, random_state=RANDOM_STATE)
    neg_mse = cross_val_score(
        model, X, y, cv=kf,
        scoring="neg_mean_squared_error", n_jobs=-1
    )
    return np.sqrt(-neg_mse.mean())
 
def tune_and_save(
    X_train: pd.DataFrame,
    y_train: pd.Series,
    n_trials: int = N_TRIALS,
    model_path: str = "models/best_model.pkl"
) -> dict:
    """チューニングを実行して最良モデルを保存する"""
    optuna.logging.set_verbosity(optuna.logging.WARNING)
 
    study = optuna.create_study(
        direction="minimize",
        sampler=optuna.samplers.TPESampler(seed=RANDOM_STATE)
    )
    print(f"Optuna チューニング開始: {n_trials} trials...")
    study.optimize(
        lambda trial: objective(trial, X_train, y_train),
        n_trials=n_trials,
        show_progress_bar=True
    )
 
    best_params = study.best_params
    best_rmse = study.best_value
 
    print(f"\n✅ チューニング完了")
    print(f"  Best RMSE: {best_rmse:.4f}")
    print(f"  Best params: {best_params}")
 
    best_model = GradientBoostingRegressor(**best_params, random_state=RANDOM_STATE)
    best_model.fit(X_train, y_train)
    joblib.dump(best_model, model_path)
    print(f"  Model saved: {model_path}")
 
    return {
        "best_params": best_params,
        "best_rmse": best_rmse,
        "model_path": model_path,
        "n_trials": len(study.trials)
    }

# src/model_evaluation.py
import shap
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import joblib
from pathlib import Path
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
 
REPORTS_DIR = Path("reports/evaluation")
REPORTS_DIR.mkdir(parents=True, exist_ok=True)
 
def evaluate_model(
    model_path: str,
    X_test: pd.DataFrame,
    y_test: pd.Series,
    output_dir: str = str(REPORTS_DIR)
) -> dict:
    """保存済みモデルをテストデータで評価してSHAP分析まで実行する"""
    model = joblib.load(model_path)
    predictions = model.predict(X_test)
 
    rmse = np.sqrt(mean_squared_error(y_test, predictions))
    mae = mean_absolute_error(y_test, predictions)
    r2 = r2_score(y_test, predictions)
 
    print(f"=== テストデータ評価結果 ===")
    print(f"  RMSE: {rmse:.4f}")
    print(f"  MAE:  {mae:.4f}")
    print(f"  R²:   {r2:.4f}")
 
    # 残差プロット
    residuals = y_test - predictions
    fig, axes = plt.subplots(1, 2, figsize=(12, 5))
    axes[0].scatter(predictions, residuals, alpha=0.5, s=10)
    axes[0].axhline(y=0, color="r", linestyle="--")
    axes[0].set_xlabel("Predicted")
    axes[0].set_ylabel("Residuals")
    axes[0].set_title("Residual Plot")
 
    min_val = min(y_test.min(), predictions.min())
    max_val = max(y_test.max(), predictions.max())
    axes[1].scatter(y_test, predictions, alpha=0.5, s=10)
    axes[1].plot([min_val, max_val], [min_val, max_val], "r--")
    axes[1].set_xlabel("Actual")
    axes[1].set_ylabel("Predicted")
    axes[1].set_title(f"Actual vs Predicted (R²={r2:.3f})")
 
    plt.tight_layout()
    fig.savefig(f"{output_dir}/residual_plots.png", dpi=150, bbox_inches="tight")
    plt.close(fig)
 
    # SHAP分析（TreeExplainer: Gradient Boosting / Random Forestに適合）
    try:
        explainer = shap.TreeExplainer(model)
        # 大きなデータセットでは200件にサンプリングして高速化
        X_shap = X_test.sample(min(200, len(X_test)), random_state=42)
        shap_values = explainer.shap_values(X_shap)
 
        fig = plt.figure(figsize=(10, 8))
        shap.summary_plot(shap_values, X_shap, show=False, plot_size=(10, 8))
        plt.tight_layout()
        plt.savefig(f"{output_dir}/shap_summary.png",
                    dpi=150, bbox_inches="tight")
        plt.close()
        print("[OK] SHAP summary plot saved.")
 
    except Exception as e:
        print(f"警告: SHAP分析エラー: {e}")
        print("  TreeExplainerが使えない場合はKernelExplainerを検討してください。")
 
    return {"rmse": rmse, "mae": mae, "r2": r2}

# src/utils.py
import os, random, numpy as np
 
RANDOM_STATE = 42
 
def set_all_seeds(seed: int = RANDOM_STATE) -> None:
    """全ライブラリの乱数シードを固定する"""
    random.seed(seed)
    np.random.seed(seed)
    os.environ["PYTHONHASHSEED"] = str(seed)
    try:
        import torch
        torch.manual_seed(seed)
        if torch.cuda.is_available():
            torch.cuda.manual_seed_all(seed)
    except ImportError:
        pass

# 大規模CSV対応のchunk読み込みパターン
def load_large_csv(
    filepath: str,
    chunksize: int = 100_000,
    target_col: str = "target"
) -> tuple[pd.DataFrame, pd.Series]:
    """大規模CSVをchunkで読み込んでメモリ効率よく処理する"""
    chunks = []
    for i, chunk in enumerate(pd.read_csv(filepath, chunksize=chunksize)):
        chunk = chunk.dropna(subset=[target_col])
        chunk["purchase_date"] = pd.to_datetime(chunk["purchase_date"])
        chunks.append(chunk)
        if (i + 1) % 10 == 0:
            print(f"  処理済み: {(i + 1) * chunksize:,} 行")
 
    df = pd.concat(chunks, ignore_index=True)
    print(f"完了: 合計 {len(df):,} 行読み込み")
    return df.drop(columns=[target_col]), df[target_col]

# src/mlops_utils.py
import json, joblib
from datetime import datetime
from pathlib import Path
from typing import Any
 
MODELS_DIR = Path("models")
REGISTRY_PATH = MODELS_DIR / "registry.json"
MODELS_DIR.mkdir(exist_ok=True)
 
def save_model_with_metadata(
    model: Any,
    metrics: dict,
    params: dict,
    model_name: str
) -> str:
    """モデルをバージョン付きで保存し、レジストリに記録する"""
    timestamp = datetime.now().strftime("%Y%m%d_%H%M%S")
    model_filename = f"{model_name}_{timestamp}.pkl"
    model_path = MODELS_DIR / model_filename
 
    joblib.dump(model, model_path)
 
    metadata = {
        "model_name": model_name,
        "timestamp": timestamp,
        "model_path": str(model_path),
        "metrics": metrics,
        "params": params,
        "model_type": type(model).__name__
    }
 
    registry = []
    if REGISTRY_PATH.exists():
        with open(REGISTRY_PATH) as f:
            registry = json.load(f)
    registry.append(metadata)
    with open(REGISTRY_PATH, "w") as f:
        json.dump(registry, f, indent=2, ensure_ascii=False)
 
    print(f"✅ 保存完了: {model_filename}")
    print(f"   RMSE: {metrics.get('rmse', 'N/A'):.4f}")
    return str(model_path)
 
def get_best_model(model_name: str, metric: str = "rmse") -> dict | None:
    """レジストリから最良モデルの情報を取得する"""
    if not REGISTRY_PATH.exists():
        return None
 
    with open(REGISTRY_PATH) as f:
        registry = json.load(f)
 
    candidates = [r for r in registry if r["model_name"] == model_name]
    if not candidates:
        return None
 
    # RMSEは小さい方が良い、R²は大きい方が良い
    reverse = metric not in ("rmse", "mae")
    return sorted(
        candidates,
        key=lambda x: x["metrics"].get(metric, float("inf")),
        reverse=reverse
    )[0]

# Claude Code へのプロンプト例:
# 「src/feature_engineering.py の engineer_features() に対する pytest を書いてください。
#  次の性質を検証してください:
#  - 出力に purchase_date や生のカテゴリ列が残っていないこと
#  - day_of_week が [0, 6] の範囲に収まること
#  - total_amount == price * quantity が全行で成り立つこと
#  - エンコード列に NaN が無いこと
#  - .copy() により入力 DataFrame が変更されないこと」

# tests/test_feature_engineering.py
import pandas as pd
import numpy as np
import pytest
from src.feature_engineering import engineer_features, validate_no_data_leakage
 
@pytest.fixture
def sample_df():
    """最小限の現実的なフィクスチャ — 10行、すべて有効な値"""
    return pd.DataFrame({
        "purchase_date": pd.date_range("2024-01-01", periods=10, freq="D").astype(str),
        "customer_age": [25, 30, 35, 40, 45, 50, 55, 60, 65, 70],
        "product_category": ["electronics", "clothing", "food"] * 3 + ["electronics"],
        "price": [10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0],
        "quantity": [1, 2, 3, 1, 2, 3, 1, 2, 3, 1],
    })
 
def test_output_has_no_raw_columns(sample_df):
    """purchase_date とカテゴリ文字列列は落とされていなければならない"""
    result = engineer_features(sample_df)
    assert "purchase_date" not in result.columns
    assert "product_category" not in result.columns
 
def test_day_of_week_range(sample_df):
    """day_of_week は [0, 6] の整数でなければならない"""
    result = engineer_features(sample_df)
    assert result["day_of_week"].between(0, 6).all(), "day_of_week が範囲外です"
 
def test_total_amount_calculation(sample_df):
    """total_amount は price × quantity と厳密に一致しなければならない"""
    result = engineer_features(sample_df)
    expected = sample_df["price"] * sample_df["quantity"]
    pd.testing.assert_series_equal(
        result["total_amount"].reset_index(drop=True),
        expected.reset_index(drop=True),
        check_names=False
    )
 
def test_input_not_mutated(sample_df):
    """engineer_features() は入力 DataFrame を変更してはならない"""
    original_cols = sample_df.columns.tolist()
    original_shape = sample_df.shape
    _ = engineer_features(sample_df)
    assert sample_df.columns.tolist() == original_cols
    assert sample_df.shape == original_shape
 
def test_no_data_leakage_invalid(sample_df):
    """学習とテストの日付が重なる分割 — リークを検出できなければならない"""
    df_train = sample_df.copy()
    df_test = sample_df.head(3).copy()  # 日付が重複
    assert validate_no_data_leakage(df_train, df_test) is False

# run_pipeline.py
"""
売上予測のエンドツーエンド ML パイプライン。
実行: python run_pipeline.py
"""
import pandas as pd
from sklearn.model_selection import train_test_split
 
from src.utils import set_all_seeds, RANDOM_STATE
from src.eda_report import generate_eda_report
from src.feature_engineering import engineer_features
from src.model_comparison import compare_models
from src.hyperparameter_tuning import tune_and_save
from src.model_evaluation import evaluate_model
from src.mlops_utils import save_model_with_metadata
 
DATA_PATH = "data/raw/sales.csv"
TARGET_COL = "revenue"
MODEL_NAME = "sales_prediction"
 
def main():
    set_all_seeds()
 
    generate_eda_report(DATA_PATH)
 
    # 前処理の前に分割する（リーク防止）
    df = pd.read_csv(DATA_PATH)
    df_train, df_test = train_test_split(
        df, test_size=0.2, random_state=RANDOM_STATE, shuffle=False  # 時系列なのでシャッフルしない
    )
 
    X_train = engineer_features(df_train.drop(columns=[TARGET_COL]))
    y_train = df_train[TARGET_COL].reset_index(drop=True)
    X_test = engineer_features(df_test.drop(columns=[TARGET_COL]))
    y_test = df_test[TARGET_COL].reset_index(drop=True)
 
    compare_models(X_train, y_train)
    tune_result = tune_and_save(X_train, y_train)
    metrics = evaluate_model(tune_result["model_path"], X_test, y_test)
 
    save_model_with_metadata(
        model=__import__("joblib").load(tune_result["model_path"]),
        metrics=metrics,
        params=tune_result["best_params"],
        model_name=MODEL_NAME
    )
 
    print(f"完了  RMSE: {metrics['rmse']:.4f} / R²: {metrics['r2']:.4f}")
 
if __name__ == "__main__":
    main()