総合演習：予測モデル構築レポート

import numpy as np
import pandas as pd
from datetime import datetime, timedelta
from sklearn.model_selection import train_test_split

np.random.seed(42)
n = 6000

# === 問題定義 ===
problem_def = """
=== NetShop 離反予測 問題定義書 ===
ビジネス目標: 顧客離反を事前に予測し、防止施策の効果を最大化する
ML問題: 30日以内の離反を予測する二値分類問題
予測対象: is_churned（1: 離反, 0: 継続）
評価指標: F1スコア（主）、再現率、AUC-ROC（副）
成功基準: F1スコア 0.75以上、再現率 0.80以上
予測頻度: 日次バッチ処理
活用方法: 離反リスクの高い顧客にクーポン/フォローメール配信
"""
print(problem_def)

# === データ生成 ===
base_date = datetime(2026, 3, 1)

df = pd.DataFrame({
    'customer_id': range(1, n + 1),
    'purchase_count': np.random.poisson(8, n),
    'avg_amount': np.random.lognormal(8, 0.8, n).astype(int),
    'days_inactive': np.random.exponential(25, n).astype(int),
    'page_views': np.random.poisson(30, n),
    'support_tickets': np.random.poisson(1.5, n),
    'membership_months': np.random.exponential(18, n).astype(int),
    'coupon_usage': np.random.poisson(2, n),
    'favorite_count': np.random.poisson(5, n),
    'device_type': np.random.choice(['PC', 'スマホ', 'タブレット'], n, p=[0.3, 0.55, 0.15]),
    'membership_rank': np.random.choice(['ブロンズ', 'シルバー', 'ゴールド', 'プラチナ'], n, p=[0.4, 0.3, 0.2, 0.1]),
    'payment_method': np.random.choice(['クレジット', '銀行振込', 'コンビニ', 'PayPay'], n, p=[0.4, 0.2, 0.15, 0.25]),
    'prefecture_region': np.random.choice(['関東', '関西', '中部', '九州', 'その他'], n, p=[0.35, 0.25, 0.15, 0.1, 0.15]),
})

# 離反ラベル
churn_prob = 1 / (1 + np.exp(
    -(df['days_inactive'] - 35) / 12
    + df['purchase_count'] * 0.08
    - df['support_tickets'] * 0.12
    + df['coupon_usage'] * 0.05
))
df['is_churned'] = (np.random.random(n) < churn_prob).astype(int)

# 欠損値
df.loc[np.random.random(n) < 0.03, 'avg_amount'] = np.nan
df.loc[np.random.random(n) < 0.02, 'page_views'] = np.nan
df.loc[np.random.random(n) < 0.01, 'favorite_count'] = np.nan

# === EDA ===
print("=== 基本情報 ===")
print(f"データ件数: {len(df)}")
print(f"離反率: {df['is_churned'].mean():.1%}")
print(f"\n欠損値:\n{df.isnull().sum()[df.isnull().sum() > 0]}")
print(f"\n基本統計量:\n{df.describe()}")

# === データ分割 ===
feature_cols = [c for c in df.columns if c not in ['customer_id', 'is_churned']]
X = df[feature_cols]
y = df['is_churned']

X_temp, X_test, y_temp, y_test = train_test_split(X, y, test_size=0.2, random_state=42, stratify=y)
X_train, X_val, y_train, y_val = train_test_split(X_temp, y_temp, test_size=0.25, random_state=42, stratify=y_temp)

print(f"\n=== データ分割 ===")
print(f"学習: {len(X_train)} ({y_train.mean():.1%})")
print(f"検証: {len(X_val)} ({y_val.mean():.1%})")
print(f"テスト: {len(X_test)} ({y_test.mean():.1%})")

from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder, OrdinalEncoder
from sklearn.impute import SimpleImputer
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier, VotingClassifier, StackingClassifier
import lightgbm as lgb
from sklearn.model_selection import cross_val_score
from sklearn.metrics import f1_score, classification_report
import optuna

# === 特徴量生成 ===
def engineer_features(df):
    df = df.copy()
    df['purchase_per_view'] = df['purchase_count'] / df['page_views'].replace(0, 1)
    df['support_per_purchase'] = df['support_tickets'] / df['purchase_count'].replace(0, 1)
    df['coupon_per_purchase'] = df['coupon_usage'] / df['purchase_count'].replace(0, 1)
    df['log_amount'] = np.log1p(df['avg_amount'])
    df['is_long_inactive'] = (df['days_inactive'] > 30).astype(int)
    df['purchase_amount_interaction'] = df['purchase_count'] * df['avg_amount']
    df['engagement_score'] = df['page_views'] + df['favorite_count'].fillna(0) + df['coupon_usage']
    return df

X_train_fe = engineer_features(X_train)
X_val_fe = engineer_features(X_val)
X_test_fe = engineer_features(X_test)

# カラム分類
numeric_features = ['purchase_count', 'avg_amount', 'days_inactive', 'page_views',
                    'support_tickets', 'membership_months', 'coupon_usage', 'favorite_count',
                    'purchase_per_view', 'support_per_purchase', 'coupon_per_purchase',
                    'log_amount', 'is_long_inactive', 'purchase_amount_interaction', 'engagement_score']
onehot_features = ['device_type', 'payment_method', 'prefecture_region']
ordinal_features = ['membership_rank']

preprocessor = ColumnTransformer([
    ('num', Pipeline([('imputer', SimpleImputer(strategy='median')), ('scaler', StandardScaler())]), numeric_features),
    ('onehot', Pipeline([('imputer', SimpleImputer(strategy='most_frequent')), ('encoder', OneHotEncoder(drop='first', sparse_output=False))]), onehot_features),
    ('ordinal', Pipeline([('imputer', SimpleImputer(strategy='most_frequent')), ('encoder', OrdinalEncoder(categories=[['ブロンズ', 'シルバー', 'ゴールド', 'プラチナ']]))]), ordinal_features),
])

X_train_proc = preprocessor.fit_transform(X_train_fe)
X_val_proc = preprocessor.transform(X_val_fe)
X_test_proc = preprocessor.transform(X_test_fe)

# === モデル構築と比較 ===
improvement_log = []

# ベースライン
lr = LogisticRegression(max_iter=1000, random_state=42)
lr.fit(X_train_proc, y_train)
f1_lr = f1_score(y_test, lr.predict(X_test_proc))
improvement_log.append(('ロジスティック回帰', f1_lr))
print(f"ベースライン F1: {f1_lr:.3f}")

# LightGBM (デフォルト)
lgbm_default = lgb.LGBMClassifier(n_estimators=200, random_state=42, verbose=-1)
lgbm_default.fit(X_train_proc, y_train)
f1_lgbm = f1_score(y_test, lgbm_default.predict(X_test_proc))
improvement_log.append(('LightGBM (デフォルト)', f1_lgbm))
print(f"LightGBM F1: {f1_lgbm:.3f}")

# Optuna 最適化
optuna.logging.set_verbosity(optuna.logging.WARNING)
def objective(trial):
    params = {
        'n_estimators': trial.suggest_int('n_estimators', 100, 400),
        'max_depth': trial.suggest_int('max_depth', 3, 10),
        'learning_rate': trial.suggest_float('learning_rate', 0.01, 0.2, log=True),
        'num_leaves': trial.suggest_int('num_leaves', 10, 60),
        'subsample': trial.suggest_float('subsample', 0.6, 1.0),
        'colsample_bytree': trial.suggest_float('colsample_bytree', 0.6, 1.0),
    }
    model = lgb.LGBMClassifier(**params, random_state=42, verbose=-1)
    scores = cross_val_score(model, X_train_proc, y_train, cv=5, scoring='f1')
    return scores.mean()

study = optuna.create_study(direction='maximize')
study.optimize(objective, n_trials=30)

lgbm_tuned = lgb.LGBMClassifier(**study.best_params, random_state=42, verbose=-1)
lgbm_tuned.fit(X_train_proc, y_train)
f1_tuned = f1_score(y_test, lgbm_tuned.predict(X_test_proc))
improvement_log.append(('LightGBM (最適化)', f1_tuned))
print(f"LightGBM (最適化) F1: {f1_tuned:.3f}")

# Stacking
from xgboost import XGBClassifier
stacking = StackingClassifier(
    estimators=[
        ('lr', LogisticRegression(max_iter=1000, random_state=42)),
        ('rf', RandomForestClassifier(n_estimators=100, random_state=42)),
        ('lgbm', lgb.LGBMClassifier(**study.best_params, random_state=42, verbose=-1)),
    ],
    final_estimator=LogisticRegression(random_state=42),
    cv=5, stack_method='predict_proba',
)
stacking.fit(X_train_proc, y_train)
f1_stacking = f1_score(y_test, stacking.predict(X_test_proc))
improvement_log.append(('Stacking', f1_stacking))
print(f"Stacking F1: {f1_stacking:.3f}")

print(f"\n=== 改善履歴 ===")
for name, score in improvement_log:
    print(f"  {name:25s}: F1 = {score:.3f}")

import shap
from sklearn.metrics import roc_curve, roc_auc_score, confusion_matrix, ConfusionMatrixDisplay
import matplotlib.pyplot as plt

# 最良モデルを選定
best_name, best_f1 = max(improvement_log, key=lambda x: x[1])
print(f"最良モデル: {best_name} (F1={best_f1:.3f})")

# LightGBM(最適化) で分析
y_prob = lgbm_tuned.predict_proba(X_test_proc)[:, 1]

# 閾値最適化
best_threshold = 0.5
best_f1_t = 0
for t in np.arange(0.15, 0.75, 0.02):
    y_pred_t = (y_prob >= t).astype(int)
    f1_t = f1_score(y_test, y_pred_t)
    if f1_t > best_f1_t:
        best_f1_t = f1_t
        best_threshold = t

y_final = (y_prob >= best_threshold).astype(int)
auc = roc_auc_score(y_test, y_prob)

# 混同行列
cm = confusion_matrix(y_test, y_final)
tn, fp, fn, tp = cm.ravel()

# SHAP分析
explainer = shap.TreeExplainer(lgbm_tuned)
shap_values = explainer.shap_values(X_test_proc)
feat_names = preprocessor.get_feature_names_out()
mean_abs_shap = np.abs(shap_values[1]).mean(axis=0)
shap_importance = pd.Series(mean_abs_shap, index=feat_names).sort_values(ascending=False)

# === ビジネスレポート ===
report = f"""
{'='*60}
  NetShop 離反予測モデル ビジネスレポート
{'='*60}

【エグゼクティブサマリー】
  機械学習を用いた顧客離反予測モデルを構築しました。
  最適モデルはLightGBMベースで、F1スコア{best_f1_t:.3f}を達成。
  離反リスクの高い顧客を事前に特定し、施策の効率化が可能です。

【モデル性能】
  F1スコア:   {best_f1_t:.3f}（閾値={best_threshold:.2f}）
  AUC-ROC:   {auc:.3f}
  適合率:    {tp/(tp+fp):.3f}
  再現率:    {tp/(tp+fn):.3f}

  混同行列:
    正しく継続と予測:   {tn}人
    誤って離反と予測:   {fp}人（不要なクーポンコスト）
    離反を見逃し:       {fn}人（顧客損失リスク）
    正しく離反と予測:   {tp}人（施策で防止可能）

【主要な離反要因（SHAP分析）】
"""
for i, (feat, imp) in enumerate(shap_importance.head(5).items(), 1):
    report += f"  {i}. {feat} (重要度: {imp:.4f})\n"

report += f"""
【推奨施策】
  1. 離反確率{best_threshold*100:.0f}%以上: 即座にクーポン配布 + 個別フォロー
     対象: {(y_prob >= best_threshold).sum()}人
  2. 離反確率30-{best_threshold*100:.0f}%: フォローメール送信
     対象: {((y_prob >= 0.3) & (y_prob < best_threshold)).sum()}人
  3. 離反確率30%未満: 通常のマーケティング施策
     対象: {(y_prob < 0.3).sum()}人

【コスト試算】
  クーポンコスト（500円 x {fp}人）: {fp * 500:,}円
  離反防止効果（{tp}人 x 年間価値30,000円）: {tp * 30000:,}円
  期待ROI: {(tp * 30000 - fp * 500) / (fp * 500) * 100:.0f}%

【今後の改善方針】
  1. 実際の施策結果をフィードバックし、モデルを再学習
  2. 時系列特徴量の追加（購入トレンド、季節性）
  3. 施策効果の因果推論分析
  4. リアルタイム予測への移行検討
{'='*60}
"""
print(report)