因果推論学習：Phase4

# ローカル環境の場合
# pip install scikit-learn xgboost causalml econml shap

# ライブラリのインポート
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split, cross_val_score
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
from sklearn.metrics import accuracy_score, mean_squared_error, roc_auc_score
import warnings
warnings.filterwarnings('ignore')

# 日本語フォント設定
plt.rcParams['font.sans-serif'] = ['Hiragino Sans', 'Yu Gothic', 'Meiryo', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("✅ 環境構築完了！")

# サンプルデータ: 顧客属性と広告への反応
np.random.seed(42)
n_samples = 1000

# 顧客属性
age = np.random.randint(20, 70, n_samples) # 年齢
income = np.random.randint(300, 1500, n_samples) # 年収（万円）
is_urban = np.random.binomial(1, 0.6, n_samples) # 都心在住（1=都心、0=郊外）

# 処置（広告を見た=1、見ていない=0）をランダム割り当て
treatment = np.random.binomial(1, 0.5, n_samples)

# 真の因果効果（CATEが顧客属性に依存）
# 若い人・高所得・都心在住ほど広告効果が高い
true_cate = (
0.1 * (70 - age) + # 若いほど効果大
0.01 * income + # 高所得ほど効果大
5 * is_urban # 都心在住で効果大
)

# 結果変数（問い合わせするか: 1=する、0=しない）
# 基準確率
base_prob = 1 / (1 + np.exp(-(0.01 * age + 0.001 * income - 1)))

# 処置効果を加算
prob = np.where(treatment == 1,
base_prob + true_cate / 100, # 広告ありの場合
base_prob) # 広告なしの場合

# 確率からサンプリング
outcome = np.random.binomial(1, np.clip(prob, 0, 1))

# DataFrameに変換
df = pd.DataFrame({
'age': age,
'income': income,
'is_urban': is_urban,
'treatment': treatment,
'outcome': outcome,
'true_cate': true_cate
})

print(f"\n=== データ概要 ===")
print(f"サンプル数: {len(df)}")
print(f"処置群: {df['treatment'].sum()} 人")
print(f"対照群: {(1-df['treatment']).sum()} 人")
print(f"問い合わせ率（全体）: {df['outcome'].mean():.2%}")
print(f"問い合わせ率（処置群）: {df[df['treatment']==1]['outcome'].mean():.2%}")
print(f"問い合わせ率（対照群）: {df[df['treatment']==0]['outcome'].mean():.2%}")

# Random Forestで問い合わせ予測
# 特徴量とターゲット
X = df[['age', 'income', 'is_urban', 'treatment']]
y = df['outcome']

# 訓練データとテストデータに分割（80:20）
X_train, X_test, y_train, y_test = train_test_split(
X, y, test_size=0.2, random_state=42
)

print(f"\n訓練データ: {len(X_train)} 件")
print(f"テストデータ: {len(X_test)} 件")

# Random Forestモデルの構築
rf_model = RandomForestClassifier(
n_estimators=100, # 決定木の数
max_depth=5, # 木の最大深さ
random_state=42
)

# 訓練
rf_model.fit(X_train, y_train)

# 予測
y_pred = rf_model.predict(X_test)
y_pred_proba = rf_model.predict_proba(X_test)[:, 1] # 確率

# 評価
accuracy = accuracy_score(y_test, y_pred)
auc = roc_auc_score(y_test, y_pred_proba)

print(f"\n=== Random Forest 予測精度 ===")
print(f"正解率（Accuracy）: {accuracy:.2%}")
print(f"AUC: {auc:.3f}")

# Feature Importance（変数重要度）
feature_importance = pd.DataFrame({
'変数': X.columns,
'重要度': rf_model.feature_importances_
}).sort_values('重要度', ascending=False)

print(f"\n=== Feature Importance ===")
print(feature_importance)

# 可視化
plt.figure(figsize=(10, 6))
plt.barh(feature_importance['変数'], feature_importance['重要度'], color='#e91e63')
plt.xlabel('重要度')
plt.title('Random Forest: Feature Importance', fontsize=14, fontweight='bold')
plt.gca().invert_yaxis()
plt.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.show()

# XGBoostモデルの構築
xgb_model = XGBClassifier(
n_estimators=100, # ブースティング回数
max_depth=3, # 木の最大深さ
learning_rate=0.1, # 学習率
random_state=42,
eval_metric='logloss' # 評価指標
)

# 訓練
xgb_model.fit(X_train, y_train)

# 予測
y_pred_xgb = xgb_model.predict(X_test)
y_pred_proba_xgb = xgb_model.predict_proba(X_test)[:, 1]

# 評価
accuracy_xgb = accuracy_score(y_test, y_pred_xgb)
auc_xgb = roc_auc_score(y_test, y_pred_proba_xgb)

print(f"=== XGBoost 予測精度 ===")
print(f"正解率（Accuracy）: {accuracy_xgb:.2%}")
print(f"AUC: {auc_xgb:.3f}")

# Random Forestと比較
print(f"\n=== Random Forest vs XGBoost ===")
print(f"Random Forest AUC: {auc:.3f}")
print(f"XGBoost AUC: {auc_xgb:.3f}")
print(f"改善: {(auc_xgb - auc) / auc * 100:.1f}%")

# Feature Importance
feature_importance_xgb = pd.DataFrame({
'変数': X.columns,
'重要度': xgb_model.feature_importances_
}).sort_values('重要度', ascending=False)

print(f"\n=== XGBoost Feature Importance ===")
print(feature_importance_xgb)

# 可視化: RF vs XGBoost
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

axes[0].barh(feature_importance['変数'], feature_importance['重要度'], color='#e91e63')
axes[0].set_xlabel('重要度')
axes[0].set_title('Random Forest', fontsize=12, fontweight='bold')
axes[0].invert_yaxis()
axes[0].grid(axis='x', alpha=0.3)

axes[1].barh(feature_importance_xgb['変数'], feature_importance_xgb['重要度'], color='#2196f3')
axes[1].set_xlabel('重要度')
axes[1].set_title('XGBoost', fontsize=12, fontweight='bold')
axes[1].invert_yaxis()
axes[1].grid(axis='x', alpha=0.3)

plt.tight_layout()
plt.show()

# S-Learnerの実装
# ステップ1: 処置変数を含めてモデル訓練
X_train_s = X_train.copy()
X_test_s = X_test.copy()

s_learner = RandomForestRegressor(n_estimators=100, max_depth=5, random_state=42)
s_learner.fit(X_train_s, y_train)

# ステップ2: CATE推定
# 各サンプルについて、treatment=1とtreatment=0の両方で予測
X_test_t1 = X_test_s.copy()
X_test_t1['treatment'] = 1 # 処置ありと仮定

X_test_t0 = X_test_s.copy()
X_test_t0['treatment'] = 0 # 処置なしと仮定

# 予測
y_pred_t1 = s_learner.predict(X_test_t1)
y_pred_t0 = s_learner.predict(X_test_t0)

# CATE = E[Y|T=1, X] - E[Y|T=0, X]
cate_s = y_pred_t1 - y_pred_t0

# テストデータに追加
test_results = X_test.copy()
test_results['true_cate'] = df.loc[X_test.index, 'true_cate']
test_results['estimated_cate_s'] = cate_s

# CATEの推定精度
from sklearn.metrics import mean_absolute_error, mean_squared_error

mae_s = mean_absolute_error(test_results['true_cate'], test_results['estimated_cate_s'])
rmse_s = np.sqrt(mean_squared_error(test_results['true_cate'], test_results['estimated_cate_s']))

print(f"\n=== S-Learner CATE推定精度 ===")
print(f"MAE: {mae_s:.2f}")
print(f"RMSE: {rmse_s:.2f}")

# 可視化: 真のCATEと推定CATE
plt.figure(figsize=(10, 6))
plt.scatter(test_results['true_cate'], test_results['estimated_cate_s'],
alpha=0.5, color='#e91e63')
plt.plot([test_results['true_cate'].min(), test_results['true_cate'].max()],
[test_results['true_cate'].min(), test_results['true_cate'].max()],
'r--', linewidth=2, label='理想的な予測')
plt.xlabel('真のCATE')
plt.ylabel('推定CATE（S-Learner）')
plt.title('S-Learner: CATE推定精度', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# セグメント別のCATE
print("\n=== セグメント別 平均CATE ===")
test_results['age_group'] = pd.cut(test_results['age'], bins=[0, 30, 50, 100],
labels=['若年層', '中年層', 'シニア層'])

segment_cate = test_results.groupby('age_group').agg({
'true_cate': 'mean',
'estimated_cate_s': 'mean'
}).round(2)

print(segment_cate)

print("\n✅ S-Learnerの利点: 実装が簡単、既存の機械学習ライブラリをそのまま使える")
print("⚠️ S-Learnerの欠点: 処置群と対照群で別々のパターンを学習できない")

# 処置群と対照群に分割
train_treated = X_train[X_train['treatment'] == 1].drop('treatment', axis=1)
train_control = X_train[X_train['treatment'] == 0].drop('treatment', axis=1)

y_train_treated = y_train[X_train['treatment'] == 1]
y_train_control = y_train[X_train['treatment'] == 0]

# 処置群用モデル
model_treated = RandomForestRegressor(n_estimators=100, max_depth=5, random_state=42)
model_treated.fit(train_treated, y_train_treated)

# 対照群用モデル
model_control = RandomForestRegressor(n_estimators=100, max_depth=5, random_state=42)
model_control.fit(train_control, y_train_control)

# テストデータでCATEを推定
X_test_features = X_test.drop('treatment', axis=1)

y_pred_treated = model_treated.predict(X_test_features)
y_pred_control = model_control.predict(X_test_features)

# CATE = 処置群の予測 - 対照群の予測
cate_t = y_pred_treated - y_pred_control

test_results['estimated_cate_t'] = cate_t

# 推定精度
mae_t = mean_absolute_error(test_results['true_cate'], test_results['estimated_cate_t'])
rmse_t = np.sqrt(mean_squared_error(test_results['true_cate'], test_results['estimated_cate_t']))

print(f"\n=== T-Learner CATE推定精度 ===")
print(f"MAE: {mae_t:.2f}")
print(f"RMSE: {rmse_t:.2f}")

print("\n✅ T-Learnerの利点: 処置群と対照群のパターンを別々に学習できる")
print("⚠️ T-Learnerの欠点: サンプル数が少ない場合、各モデルの精度が低下")

# X-Learner（簡易版の概念説明）
print("\n=== X-Learner（Extended Learner）===")
print("特徴: T-Learnerの改良版、傾向スコアで重み付け")
print("※ 完全な実装はCausalMLライブラリで学習します")

# Meta-Learnersの比較
print("\n=== Meta-Learners比較 ===")
comparison = pd.DataFrame({
'S-Learner': [mae_s, rmse_s],
'T-Learner': [mae_t, rmse_t]
}, index=['MAE', 'RMSE']).round(2)

print(comparison)

# 可視化: S-Learner vs T-Learner
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

axes[0].scatter(test_results['true_cate'], test_results['estimated_cate_s'],
alpha=0.5, color='#e91e63')
axes[0].plot([test_results['true_cate'].min(), test_results['true_cate'].max()],
[test_results['true_cate'].min(), test_results['true_cate'].max()],
'r--', linewidth=2)
axes[0].set_xlabel('真のCATE')
axes[0].set_ylabel('推定CATE')
axes[0].set_title(f'S-Learner (MAE={mae_s:.2f})', fontsize=12, fontweight='bold')
axes[0].grid(alpha=0.3)

axes[1].scatter(test_results['true_cate'], test_results['estimated_cate_t'],
alpha=0.5, color='#2196f3')
axes[1].plot([test_results['true_cate'].min(), test_results['true_cate'].max()],
[test_results['true_cate'].min(), test_results['true_cate'].max()],
'r--', linewidth=2)
axes[1].set_xlabel('真のCATE')
axes[1].set_ylabel('推定CATE')
axes[1].set_title(f'T-Learner (MAE={mae_t:.2f})', fontsize=12, fontweight='bold')
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

print("=== CausalMLでMeta-Learners ===")

# データ準備
X_cm = df[['age', 'income', 'is_urban']].values
treatment_cm = df['treatment'].values
y_cm = df['outcome'].values

# 訓練/テストデータ分割
from sklearn.model_selection import train_test_split
X_train_cm, X_test_cm, t_train, t_test, y_train_cm, y_test_cm = train_test_split(
X_cm, treatment_cm, y_cm, test_size=0.2, random_state=42
)

# ベースモデル（XGBoost）
base_model = XGBRegressor(n_estimators=100, max_depth=3, random_state=42)

# S-Learner
learner_s = BaseSLearner(learner=base_model)
learner_s.fit(X=X_train_cm, treatment=t_train, y=y_train_cm)
cate_pred_s = learner_s.predict(X=X_test_cm)

# T-Learner
learner_t = BaseTLearner(learner=base_model)
learner_t.fit(X=X_train_cm, treatment=t_train, y=y_train_cm)
cate_pred_t = learner_t.predict(X=X_test_cm)

# X-Learner
learner_x = BaseXLearner(learner=base_model)
learner_x.fit(X=X_train_cm, treatment=t_train, y=y_train_cm)
cate_pred_x = learner_x.predict(X=X_test_cm)

# 真のCATEと比較
true_cate_test = df.loc[df.index[len(X_train_cm):], 'true_cate'].values

# 各手法の精度
results_cm = pd.DataFrame({
'S-Learner': [
mean_absolute_error(true_cate_test, cate_pred_s.flatten()),
np.sqrt(mean_squared_error(true_cate_test, cate_pred_s.flatten()))
],
'T-Learner': [
mean_absolute_error(true_cate_test, cate_pred_t.flatten()),
np.sqrt(mean_squared_error(true_cate_test, cate_pred_t.flatten()))
],
'X-Learner': [
mean_absolute_error(true_cate_test, cate_pred_x.flatten()),
np.sqrt(mean_squared_error(true_cate_test, cate_pred_x.flatten()))
]
}, index=['MAE', 'RMSE']).round(3)

print("\n=== CausalML Meta-Learners 精度比較 ===")
print(results_cm)

# 最良モデル
best_model = results_cm.loc['MAE'].idxmin()
print(f"\n✅ 最良モデル: {best_model}")

# 可視化
fig, axes = plt.subplots(1, 3, figsize=(18, 5))

for i, (learner_name, cate_pred) in enumerate([
('S-Learner', cate_pred_s),
('T-Learner', cate_pred_t),
('X-Learner', cate_pred_x)
]):
axes[i].scatter(true_cate_test, cate_pred.flatten(), alpha=0.5)
axes[i].plot([true_cate_test.min(), true_cate_test.max()],
[true_cate_test.min(), true_cate_test.max()],
'r--', linewidth=2)
axes[i].set_xlabel('真のCATE')
axes[i].set_ylabel('推定CATE')
axes[i].set_title(f'{learner_name}\nMAE={results_cm.loc["MAE", learner_name]:.3f}',
fontsize=12, fontweight='bold')
axes[i].grid(alpha=0.3)

plt.tight_layout()
plt.show()

print("=== アップリフト評価指標 ===")

# Gainカーブ（アップリフトカーブ）の作成
# 推定CATEの高い順に顧客をターゲティングした場合の累積効果

# T-Learnerの結果を使用
cate_pred_flat = cate_pred_t.flatten()

# Gainカーブのプロット
fig, ax = plt.subplots(figsize=(10, 6))

plot_gain(y_test_cm, cate_pred_flat, t_test,
plot_type='gain', ax=ax)

plt.title('Gain Curve（ゲインカーブ）', fontsize=14, fontweight='bold')
plt.xlabel('ターゲティング比率')
plt.ylabel('累積Gain')
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

# Qiniカーブ
fig, ax = plt.subplots(figsize=(10, 6))

plot_gain(y_test_cm, cate_pred_flat, t_test,
plot_type='qini', ax=ax)

plt.title('Qini Curve（キニカーブ）', fontsize=14, fontweight='bold')
plt.xlabel('ターゲティング比率')
plt.ylabel('Qini係数')
plt.grid(alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

# Qini係数とAUUCの計算
qini = qini_score(y_test_cm, cate_pred_flat, t_test)
auuc = auuc_score(y_test_cm, cate_pred_flat, t_test)

print(f"\n=== 評価指標 ===")
print(f"Qini係数: {qini:.3f}")
print(f"AUUC: {auuc:.3f}")

# 実務的解釈
print("\n=== 実務的解釈 ===")
print("上位20%の顧客にターゲティングした場合の期待効果:")

# CATEで顧客をソート
test_df_eval = pd.DataFrame({
'cate': cate_pred_flat,
'treatment': t_test,
'outcome': y_test_cm
})

test_df_eval = test_df_eval.sort_values('cate', ascending=False)

# 上位20%
top_20_pct = int(len(test_df_eval) * 0.2)
top_20_df = test_df_eval.head(top_20_pct)

# 効果の計算
treated_top20 = top_20_df[top_20_df['treatment'] == 1]
control_top20 = top_20_df[top_20_df['treatment'] == 0]

if len(treated_top20) > 0 and len(control_top20) > 0:
effect_top20 = treated_top20['outcome'].mean() - control_top20['outcome'].mean()
print(f" 上位20%のCATE: {effect_top20:.2%}")
print(f" 全体平均との比較: {effect_top20 / (y_test_cm[t_test==1].mean() - y_test_cm[t_test==0].mean()):.1f}倍")

# ランダムターゲティングとの比較
print("\nランダムターゲティング vs CATE順ターゲティング:")
print(" → CATEを使うことで、同じ予算でより高い効果が期待できる")

print("=== Doubly Robust Learner ===")

# データ準備（Week 1-2と同じデータを使用）
X_dr = df[['age', 'income', 'is_urban']].values
treatment_dr = df['treatment'].values
y_dr = df['outcome'].values

# 訓練/テストデータ分割
X_train_dr, X_test_dr, t_train_dr, t_test_dr, y_train_dr, y_test_dr = train_test_split(
X_dr, treatment_dr, y_dr, test_size=0.2, random_state=42
)

# ベースモデル
base_model_dr = XGBRegressor(n_estimators=100, max_depth=3, random_state=42)

# DRLearner
learner_dr = BaseDRLearner(learner=base_model_dr)
learner_dr.fit(X=X_train_dr, treatment=t_train_dr, y=y_train_dr)

# CATE推定
cate_pred_dr = learner_dr.predict(X=X_test_dr)

# 真のCATEと比較
true_cate_dr = df.loc[df.index[len(X_train_dr):], 'true_cate'].values

mae_dr = mean_absolute_error(true_cate_dr, cate_pred_dr
.flatten())
rmse_dr = np.sqrt(mean_squared_error(true_cate_dr, cate_pred_dr.flatten()))

print(f"\n=== DRLearner CATE推定精度 ===")
print(f"MAE: {mae_dr:.3f}")
print(f"RMSE: {rmse_dr:.3f}")

# T-Learnerと比較
print(f"\n=== T-Learner vs DRLearner ===")
print(f"T-Learner MAE: {mae_t:.3f}")
print(f"DRLearner MAE: {mae_dr:.3f}")
print(f"改善率: {(mae_t - mae_dr) / mae_t * 100:.1f}%")

# 可視化
plt.figure(figsize=(10, 6))
plt.scatter(true_cate_dr, cate_pred_dr.flatten(), alpha=0.5, color='#e91e63')
plt.plot([true_cate_dr.min(), true_cate_dr.max()],
[true_cate_dr.min(), true_cate_dr.max()],
'r--', linewidth=2, label='理想的な予測')
plt.xlabel('真のCATE')
plt.ylabel('推定CATE（DRLearner）')
plt.title('DRLearner: CATE推定精度', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

# 傾向スコアの確認（内部的に計算されている）
print("\n=== Doubly Robust推定の仕組み ===")
print("1. 傾向スコアモデル: 処置を受ける確率 P(T=1|X) を推定")
print("2. 結果モデル: アウトカム E[Y|T,X] を推定")
print("3. DR推定量: 両方を組み合わせてCATEを計算")
print("4. どちらか一方が正しければ一貫性が保証される")

# バイアス削減の実証（シミュレーション）
print("\n=== バイアス削減の実証 ===")

# 意図的にバイアスのあるデータを生成
# 高所得者ほど処置を受けやすい（選択バイアス）
biased_treatment_prob = 1 / (1 + np.exp(-(0.002 * income - 1)))
biased_treatment = np.random.binomial(1, biased_treatment_prob)

df_biased = df.copy()
df_biased['treatment'] = biased_treatment

# バイアスのあるデータでナイーブ推定（単純な平均の差）
ate_naive = (df_biased[df_biased['treatment']==1]['outcome'].mean() -
df_biased[df_biased['treatment']==0]['outcome'].mean())

print(f"ナイーブ推定（選択バイアスあり）: {ate_naive:.3f}")
print(f"真のATE（全体平均）: {df['true_cate'].mean():.3f}")
print(f"バイアス: {abs(ate_naive - df['true_cate'].mean()):.3f}")

# DRLearnerで推定（バイアス削減）
X_biased = df_biased[['age', 'income', 'is_urban']].values
t_biased = df_biased['treatment'].values
y_biased = df_biased['outcome'].values

learner_dr_biased = BaseDRLearner(learner=XGBRegressor(n_estimators=100, max_depth=3, random_state=42))
learner_dr_biased.fit(X=X_biased, treatment=t_biased, y=y_biased)
cate_dr_biased = learner_dr_biased.predict(X=X_biased)

ate_dr = cate_dr_biased.mean()

print(f"DRLearner推定（バイアス削減後）: {ate_dr:.3f}")
print(f"バイアス削減: {abs(ate_naive - df['true_cate'].mean()) - abs(ate_dr - df['true_cate'].mean()):.3f}")
print("\n✅ DRLearnerは選択バイアスに強い！")

print("=== EconML: Double Machine Learning ===")

# データ準備
X_dml = df[['age', 'income', 'is_urban']].values
T_dml = df['treatment'].values.reshape(-1, 1) # EconMLは2次元を期待
Y_dml = df['outcome'].values

# LinearDML
# model_y: 結果モデル（Y を X で予測）
# model_t: 処置モデル（T を X で予測）
dml = LinearDML(
model_y=RandomForestRegressor(n_estimators=100, max_depth=5, random_state=42),
model_t=RandomForestClassifier(n_estimators=100, max_depth=5, random_state=42),
discrete_treatment=True,
cv=3 # クロスフィッティングの分割数
)

# 訓練
dml.fit(Y=Y_dml, T=T_dml, X=X_dml)

# CATE推定
cate_dml = dml.effect(X_dml)

# 信頼区間の計算
cate_interval = dml.effect_interval(X_dml, alpha=0.05) # 95%信頼区間

print(f"\n=== DML CATE推定 ===")
print(f"平均CATE: {cate_dml.mean():.3f}")
print(f"CATE範囲: [{cate_dml.min():.3f}, {cate_dml.max():.3f}]")

# 真のCATEと比較
mae_dml = mean_absolute_error(df['true_cate'], cate_dml)
rmse_dml = np.sqrt(mean_squared_error(df['true_cate'], cate_dml))

print(f"\n=== DML精度 ===")
print(f"MAE: {mae_dml:.3f}")
print(f"RMSE: {rmse_dml:.3f}")

# セグメント別のCATE（95%信頼区間付き）
print(f"\n=== セグメント別CATE（95%信頼区間付き）===")

# 年齢層別
df_dml = df.copy()
df_dml['cate_dml'] = cate_dml
df_dml['cate_lower'] = cate_interval[0]
df_dml['cate_upper'] = cate_interval[1]
df_dml['age_group'] = pd.cut(df_dml['age'], bins=[0, 30, 50, 100],
labels=['若年層', '中年層', 'シニア層'])

segment_cate_dml = df_dml.groupby('age_group').agg({
'cate_dml': 'mean',
'cate_lower': 'mean',
'cate_upper': 'mean'
}).round(2)

print(segment_cate_dml)

# 可視化: 信頼区間付きCATE
plt.figure(figsize=(12, 6))

# サンプルを100件に絞って可視化（見やすくするため）
sample_indices = np.random.choice(len(df_dml), 100, replace=False)
df_sample = df_dml.iloc[sample_indices].sort_values('cate_dml')

plt.errorbar(range(len(df_sample)),
df_sample['cate_dml'],
yerr=[df_sample['cate_dml'] - df_sample['cate_lower'],
df_sample['cate_upper'] - df_sample['cate_dml']],
fmt='o', color='#e91e63', alpha=0.6, capsize=3)

plt.axhline(y=0, color='red', linestyle='--', alpha=0.5)
plt.xlabel('サンプル（CATE順）')
plt.ylabel('推定CATE（95%信頼区間付き）')
plt.title('DML: 個別CATEの推定（信頼区間付き）', fontsize=14, fontweight='bold')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

print("\n✅ DMLの利点: 機械学習の予測力 + 統計的推論（信頼区間）")
print("✅ 信頼区間が0を含まないセグメント → 統計的に有意な効果あり")

# 高次元データの生成（50個の特徴量）
np.random.seed(42)
n = 1000
p = 50 # 特徴量の数

# ランダムな特徴量
X_high = np.random.randn(n, p)

# 真の因果効果に関連する特徴量は最初の5つだけ
true_beta = np.zeros(p)
true_beta[:5] = [2, -1.5, 1, 0.5, -0.8]

# 処置割り当て（最初の3つの特徴量に依存）
treatment_logit = X_high[:, 0] * 0.5 + X_high[:, 1] * 0.3 - X_high[:, 2] * 0.2
treatment_prob = 1 / (1 + np.exp(-treatment_logit))
T_high = np.random.binomial(1, treatment_prob).reshape(-1, 1)

# 結果変数（処置効果は特徴量に依存）
true_cate_high = X_high @ true_beta
Y_high = (
X_high[:, 0] * 0.5 + # ベースライン効果
X_high[:, 1] * 0.3 +
T_high.flatten() * true_cate_high + # 処置効果
np.random.randn(n) * 0.5 # ノイズ
)

# DMLで推定
dml_high = LinearDML(
model_y=RandomForestRegressor(n_estimators=100, max_depth=5, random_state=42),
model_t=RandomForestClassifier(n_estimators=100, max_depth=5, random_state=42),
discrete_treatment=True,
cv=3
)

dml_high.fit(Y=Y_high, T=T_high, X=X_high)

# CATE推定
cate_high_pred = dml_high.effect(X_high)

# 精度評価
mae_high = mean_absolute_error(true_cate_high, cate_high_pred)
rmse_high = np.sqrt(mean_squared_error(true_cate_high, cate_high_pred))

print(f"\n=== 高次元データ（50特徴量）でのDML精度 ===")
print(f"MAE: {mae_high:.3f}")
print(f"RMSE: {rmse_high:.3f}")

# 可視化
plt.figure(figsize=(10, 6))
plt.scatter(true_cate_high, cate_high_pred, alpha=0.5, color='#e91e63')
plt.plot([true_cate_high.min(), true_cate_high.max()],
[true_cate_high.min(), true_cate_high.max()],
'r--', linewidth=2, label='理想的な予測')
plt.xlabel('真のCATE')
plt.ylabel('推定CATE（DML）')
plt.title('高次元データでのDML推定', fontsize=14, fontweight='bold')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()

print("\n=== DMLが高次元データで強い理由 ===")
print("1. 機械学習モデルが多数の特徴量を扱える")
print("2. クロスフィッティングで過学習を防ぐ")
print("3. ニューサンスパラメータの推定を分離")
print("4. 従来の線形回帰では不可能な規模のデータを扱える")

print("=== SHAP値の計算と可視化 ===")

# Random Forestモデルで問い合わせ予測（Week 1-2のモデルを使用）
X_shap = df[['age', 'income', 'is_urban', 'treatment']]
y_shap = df['outcome']

# 訓練
rf_shap = RandomForestClassifier(n_estimators=100, max_depth=5, random_state=42)
rf_shap.fit(X_shap, y_shap)

# TreeExplainer（Random Forest用）
explainer = shap.TreeExplainer(rf_shap)

# SHAP値の計算（100サンプルで計算、全データだと時間がかかる）
shap_values = explainer.shap_values(X_shap.iloc[:100])

# クラス1（問い合わせあり）のSHAP値を使用
if isinstance(shap_values, list):
shap_values_class1 = shap_values[1]
else:
shap_values_class1 = shap_values

# Summary Plot（全体的な特徴量の重要度）
plt.figure(figsize=(10, 6))
shap.summary_plot(shap_values_class1, X_shap.iloc[:100],
plot_type="bar", show=False)
plt.title('SHAP Summary Plot（平均絶対SHAP値）', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# Beeswarm Plot（値の分布も含む）
plt.figure(figsize=(10, 6))
shap.summary_plot(shap_values_class1, X_shap.iloc[:100], show=False)
plt.title('SHAP Summary Plot（分布付き）', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.show()

# Dependence Plot（特定の特徴量の効果）
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

for i, feature in enumerate(X_shap.columns):
ax = axes[i // 2, i % 2]
shap.dependence_plot(feature, shap_values_class1, X_shap.iloc[:100],
ax=ax, show=False)
ax.set_title(f'{feature} のSHAP Dependence Plot', fontsize=12, fontweight='bold')

plt.tight_layout()
plt.show()

print("\n=== SHAP値の解釈 ===")
print("1. Summary Plot: どの特徴量が予測に重要か")
print("2. Beeswarm Plot: 特徴量の値が高い/低いで効果が変わるか")
print("3. Dependence Plot: 特徴量と予測の関係（非線形も捉える）")

# Feature Importance
feature_importance_rf = pd.DataFrame({
'特徴量': X_shap.columns,
'Feature Importance': rf_shap.feature_importances_
}).sort_values('Feature Importance', ascending=False)

# SHAP値の平均絶対値
shap_importance = pd.DataFrame({
'特徴量': X_shap.columns,
'SHAP値（平均絶対値）': np.abs(shap_values_class1).mean(axis=0)
}).sort_values('SHAP値（平均絶対値）', ascending=False)

print("\n=== Feature Importance ===")
print(feature_importance_rf)

print("\n=== SHAP値（平均絶対値）===")
print(shap_importance)

# 可視化: 比較
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

axes[0].barh(feature_importance_rf['特徴量'],
feature_importance_rf['Feature Importance'],
color='#e91e63')
axes[0].set_xlabel('Feature Importance')
axes[0].set_title('Feature Importance（不純度減少）', fontsize=12, fontweight='bold')
axes[0].invert_yaxis()
axes[0].grid(axis='x', alpha=0.3)

axes[1].barh(shap_importance['特徴量'],
shap_importance['SHAP値（平均絶対値）'],
color='#2196f3')
axes[1].set_xlabel('SHAP値（平均絶対値）')
axes[1].set_title('SHAP値（Shapley値
）', fontsize=12, fontweight='bold')
axes[1].invert_yaxis()
axes[1].grid(axis='x', alpha=0.3)

plt.tight_layout()
plt.show()

# 交絡のあるデータ生成
np.random.seed(42)
n = 1000

# 交絡変数: 年齢
age_confound = np.random.randint(20, 70, n)

# 年齢が高いほど広告を受けやすい（選択バイアス）
treatment_prob_confound = 1 / (1 + np.exp(-(0.05 * age_confound - 2)))
treatment_confound = np.random.binomial(1, treatment_prob_confound)

# 年齢が高いほど問い合わせしやすい（年齢の直接効果）
# 広告の真の因果効果は小さい
outcome_prob = (
0.01 * age_confound + # 年齢の直接効果（大）
0.05 * treatment_confound + # 広告の因果効果（小）
np.random.randn(n) * 0.1
)
outcome_confound = (outcome_prob > np.median(outcome_prob)).astype(int)

df_confound = pd.DataFrame({
'age': age_confound,
'treatment': treatment_confound,
'outcome': outcome_confound
})

# 予測モデル（SHAP値用）
X_confound = df_confound[['age', 'treatment']]
y_confound = df_confound['outcome']

rf_confound = RandomForestClassifier(n_estimators=100, random_state=42)
rf_confound.fit(X_confound, y_confound)

# SHAP値の計算
explainer_confound = shap.TreeExplainer(rf_confound)
shap_values_confound = explainer_confound.shap_values(X_confound)

if isinstance(shap_values_confound, list):
shap_values_confound = shap_values_confound[1]

# 平均SHAP値
shap_mean = pd.DataFrame({
'特徴量': X_confound.columns,
'平均SHAP値': np.abs(shap_values_confound).mean(axis=0)
}).sort_values('平均SHAP値', ascending=False)

print("\n=== 予測モデルのSHAP値 ===")
print(shap_mean)
print("\n→ 「年齢」が最も予測に貢献")

# 因果効果（ATE）の推定
ate_naive = (df_confound[df_confound['treatment']==1]['outcome'].mean() -
df_confound[df_confound['treatment']==0]['outcome'].mean())

print(f"\n=== ナイーブなATE推定 ===")
print(f"広告ありの問い合わせ率: {df_confound[df_confound['treatment']==1]['outcome'].mean():.2%}")
print(f"広告なしの問い合わせ率: {df_confound[df_confound['treatment']==0]['outcome'].mean():.2%}")
print(f"見かけ上のATE: {ate_naive:.2%}")
print("\n→ 広告の効果が大きく見える（交絡のため）")

# 真の因果効果（設計時のパラメータ）
print(f"\n=== 真の因果効果 ===")
print(f"真のATE: 5% （データ生成時に設定）")
print(f"年齢の直接効果: 大（0.01 × 年齢）")

print("\n=== 結論 ===")
print("✅ SHAP値: 「年齢」が予測に最重要 → 正しい（予測の観点）")
print("⚠️ しかし: 「年齢を変えれば問い合わせが増える」は因果関係ではない")
print("✅ 因果推論: 広告の真の効果を測定するには、RCTやDMLが必要")
print("\n💡 実務では: SHAP値（予測の説明）と因果推論（施策の効果測定）を使い分ける")

print("=== ハイパーパラメータチューニング ===")

# データ準備（Week 1-2のデータを使用）
X_tune = df[['age', 'income', 'is_urban', 'treatment']]
y_tune = df['outcome']

# ベースラインモデル（デフォルトパラメータ）
rf_baseline = RandomForestClassifier(random_state=42)
baseline_scores = cross_val_score(rf_baseline, X_tune, y_tune, cv=5, scoring='roc_auc')

print(f"\n=== ベースラインモデル（デフォルトパラメータ）===")
print(f"AUC（平均）: {baseline_scores.mean():.3f} ± {baseline_scores.std():.3f}")

# GridSearchCV: パラメータの全組み合わせを試す
param_grid = {
'n_estimators': [50, 100, 200],
'max_depth': [3, 5, 7, None],
'min_samples_split': [2, 5, 10],
'min_samples_leaf': [1, 2, 4]
}

print(f"\n=== GridSearchCV（全組み合わせ: {3*4*3*3}通り）===")
print("探索中...")

grid_search = GridSearchCV(
RandomForestClassifier(random_state=42),
param_grid,
cv=5,
scoring='roc_auc',
n_jobs=-1, # 並列処理
verbose=0
)

grid_search.fit(X_tune, y_tune)

print(f"\n最良パラメータ: {grid_search.best_params_}")
print(f"最良AUC: {grid_search.best_score_:.3f}")
print(f"改善: {(grid_search.best_score_ - baseline_scores.mean()) / baseline_scores.mean() * 100:.1f}%")

# RandomizedSearchCV: ランダムにサンプリング（高速）
from scipy.stats import randint

param_dist = {
'n_estimators': randint(50, 300),
'max_depth': [3, 5, 7, 10, None],
'min_samples_split': randint(2, 20),
'min_samples_leaf': randint(1, 10)
}

print(f"\n=== RandomizedSearchCV（ランダムサンプリング: 50回）===")
print("探索中...")

random_search = RandomizedSearchCV(
RandomForestClassifier(random_state=42),
param_dist,
n_iter=50, # 試行回数
cv=5,
scoring='roc_auc',
n_jobs=-1,
random_state=42,
verbose=0
)

random_search.fit(X_tune, y_tune)

print(f"\n最良パラメータ: {random_search.best_params_}")
print(f"最良AUC: {random_search.best_score_:.3f}")

# 最良モデルで予測
best_model = grid_search.best_estimator_
y_pred_tuned = best_model.predict(X_tune)
y_pred_proba_tuned = best_model.predict_proba(X_tune)[:, 1]

# 評価
auc_tuned = roc_auc_score(y_tune, y_pred_proba_tuned)

print(f"\n=== チューニング後のモデル性能 ===")
print(f"訓練データAUC: {auc_tuned:.3f}")

# Feature Importanceの確認
feature_importance_tuned = pd.DataFrame({
'特徴量': X_tune.columns,
'重要度': best_model.feature_importances_
}).sort_values('重要度', ascending=False)

print(f"\n=== チューニング後のFeature Importance ===")
print(feature_importance_tuned)

# 可視化: チューニング前後の比較
plt.figure(figsize=(10, 6))
models = ['ベースライン', 'GridSearch', 'RandomizedSearch']
scores = [baseline_scores.mean(), grid_search.best_score_, random_search.best_score_]

plt.bar(models, scores, color=['#9c27b0', '#e91e63', '#2196f3'])
plt.ylabel('AUC')
plt.title('ハイパーパラメータチューニングの効果', fontsize=14, fontweight='bold')
plt.ylim([min(scores) - 0.02, max(scores) + 0.02])
plt.grid(axis='y', alpha=0.3)

for i, score in enumerate(scores):
plt.text(i, score + 0.005, f'{score:.3f}', ha='center', fontweight='bold')

plt.tight_layout()
plt.show()

print("\n=== チューニングのベストプラクティス ===")
print("1. まずはデフォルトパラメータで試す（ベースライン）")
print("2. 重要なパラメータを絞ってGridSearchCV")
print("3. 時間がない場合はRandomizedSearchCV")
print("4. 過学習に注意（訓練AUCとテストAUCの差を確認）")
print("5. 実務では計算時間とのトレードオフを考慮")

# ステップ1: データ準備
print("\n=== ステップ1: データ準備 ===")
# （実際は自社データを読み込む）
# df_project = pd.read_csv('marketing_data.csv')

# ステップ2: 探索的データ分析
print("\n=== ステップ2: EDA ===")
print("- 顧客属性の分布確認")
print("- 処置群/対照群のバランス確認")
print("- アウトカムの基本統計量")

# ステップ3: アップリフトモデリング
print("\n=== ステップ3: アップリフトモデリング ===")
print("Meta-Learners（S/T/X-Learner）で顧客セグメント別CATE推定")
# （Week 1-2のコードを活用）

# ステップ4: Doubly Robust推定
print("\n=== ステップ4: Doubly Robust推定 ===")
print("選択バイアスに対処し、頑健な因果効果を推定")
# （Week 3のコードを活用）

# ステップ5: SHAP値による解釈
print("\n=== ステップ5: SHAP値による予測モデル解釈 ===")
print("「どの属性が問い合わせ予測に効くか」を可視化")
# （Week 4のコードを活用）

# ステップ6: ビジネス示唆
print("\n=== ステップ6: ビジネス示唆 ===")
print("\n【発見事項】")
print("1. CATE分析の結果:")
print(" - 30代・都心在住・高所得者: CATE = +12% → 高効果セグメント")
print(" - 60代・郊外在住・中所得者: CATE = +2% → 低効果セグメント")
print(" - 一部セグメントで逆効果（CATE < 0）を確認") print("\n2. 予測モデル（SHAP値）の結果:") print(" - 問い合わせ予測に最も効く変数: 年齢、所得、都心在住") print(" - ただし因果効果（広告を出す効果）とは異なる") print("\n【アクションプラン】") print("1. 高効果セグメント（30代・都心・高所得）に予算の60%を集中") print("2. 逆効果セグメントには広告を停止（コスト削減）") print("3. 中間セグメントは小規模テストで効果を検証") print("4. 予測モデルで「問い合わせしやすい顧客」を特定し優先フォロー") print("\n【期待効果】") print("- CV数: +25%（現状100件/月 → 125件/月）") print("- 広告費: -15%（無駄な配信を削減）") print("- ROI: +47%（効率的なターゲティング）") print("\n【次のステップ】") print("- Phase 5でLLM×因果推論を学習") print("- 自動レポート生成、説明文の自動作成") print("- より高度な因果探索（VAR-LiNGAMなど）")