機械学習の全体地図（8）通しチュートリアル：解約予測

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
n = 3000
df = pd.DataFrame({
    "契約月数": rng.integers(1, 73, n),
    "月額料金": rng.normal(6000, 2000, n).clip(1500, 15000).round(-1),
    "プラン": rng.choice(["ライト", "スタンダード", "プレミアム"], n, p=[.4, .4, .2]),
    "支払方法": rng.choice(["カード", "口座振替", "コンビニ"], n, p=[.5, .35, .15]),
    "問い合わせ回数": rng.poisson(1.5, n),
})
# 解約のしやすさ（現実には分からない「本当の仕組み」）
z = (-1.6 - 0.045 * df["契約月数"] + 0.00022 * df["月額料金"]
     + 0.45 * df["問い合わせ回数"] + 0.9 * (df["支払方法"] == "コンビニ")
     - 0.6 * (df["プラン"] == "プレミアム"))
df["解約"] = (rng.random(n) < 1 / (1 + np.exp(-z))).astype(int)
# 現実のデータらしく欠損を混ぜる
df.loc[rng.random(n) < 0.05, "月額料金"] = np.nan
df.loc[rng.random(n) < 0.03, "支払方法"] = np.nan

print(df.shape)
print(df.head(3))

from sklearn.model_selection import train_test_split

X = df.drop(columns="解約")
y = df["解約"]
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=0)
print("学習用", X_train.shape, "テスト用", X_test.shape)
print("解約率 学習用", round(y_train.mean(), 3), "テスト用", round(y_test.mean(), 3))

train = X_train.assign(解約=y_train)
print("欠損の数:\n", train.isna().sum()[lambda s: s > 0])
print("\n支払方法ごとの解約率:\n", train.groupby("支払方法")["解約"].mean().round(3))
print("\n契約月数（解約した人／しなかった人の中央値）:\n",
      train.groupby("解約")["契約月数"].median())

from sklearn.preprocessing import FunctionTransformer, StandardScaler, OneHotEncoder
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import make_pipeline
from sklearn.impute import SimpleImputer

def add_features(d):
    d = d.copy()
    d["累計支払額"] = d["契約月数"] * d["月額料金"]
    d["月あたり問い合わせ"] = d["問い合わせ回数"] / d["契約月数"]
    return d

num_cols = ["契約月数", "月額料金", "問い合わせ回数", "累計支払額", "月あたり問い合わせ"]
cat_cols = ["プラン", "支払方法"]

def make_preprocess():
    return make_pipeline(
        FunctionTransformer(add_features),
        ColumnTransformer([
            ("num", make_pipeline(SimpleImputer(strategy="median", add_indicator=True),
                                  StandardScaler()), num_cols),
            ("cat", make_pipeline(SimpleImputer(strategy="most_frequent"),
                                  OneHotEncoder(handle_unknown="ignore")), cat_cols),
        ]))

print(make_preprocess().fit(X_train).transform(X_train.head(2)).round(2))

from sklearn.dummy import DummyClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import HistGradientBoostingClassifier
from sklearn.model_selection import StratifiedKFold, cross_val_score

cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=0)
candidates = {
    "ベースライン": DummyClassifier(strategy="prior"),
    "ロジスティック回帰": LogisticRegression(max_iter=1000),
    "勾配ブースティング": HistGradientBoostingClassifier(random_state=0),
}
for name, clf in candidates.items():
    model = make_pipeline(make_preprocess(), clf)
    s = cross_val_score(model, X_train, y_train, cv=cv, scoring="average_precision")
    print(f"{name:<10} PR-AUC {s.mean():.3f} (±{s.std():.3f})")

from sklearn.model_selection import cross_val_predict
from sklearn.metrics import classification_report
from sklearn.inspection import permutation_importance

model = make_pipeline(make_preprocess(), LogisticRegression(max_iter=1000))
proba = cross_val_predict(model, X_train, y_train, cv=cv, method="predict_proba")[:, 1]
print(classification_report(y_train, proba >= 0.5, digits=3))

# 学習用の中をさらに分けて、どの項目が効いているかを見る
X_tr, X_va, y_tr, y_va = train_test_split(X_train, y_train, test_size=0.25,
                                          stratify=y_train, random_state=0)
r = permutation_importance(model.fit(X_tr, y_tr), X_va, y_va,
                           scoring="average_precision", n_repeats=10, random_state=0)
print(pd.Series(r.importances_mean, index=X_va.columns).sort_values(ascending=False).round(3))

from sklearn.model_selection import GridSearchCV, TunedThresholdClassifierCV
from sklearn.metrics import average_precision_score, precision_score, recall_score

search = GridSearchCV(
    make_pipeline(make_preprocess(), LogisticRegression(max_iter=1000)),
    {"logisticregression__C": [0.01, 0.1, 1, 10]},
    cv=cv, scoring="average_precision").fit(X_train, y_train)
print("最良のC:", search.best_params_, "PR-AUC:", round(search.best_score_, 3))

# 見逃しと誤報のバランス（F1）が最良になる閾値を、学習用データの交差検証で選ぶ
final = TunedThresholdClassifierCV(search.best_estimator_, scoring="f1", cv=cv).fit(X_train, y_train)
print("選ばれた閾値:", round(final.best_threshold_, 3))

# テスト用データで、最後に1回だけ評価する
p_test = final.estimator_.predict_proba(X_test)[:, 1]
pred = final.predict(X_test)
print(f"テスト PR-AUC {average_precision_score(y_test, p_test):.3f} "
      f"適合率 {precision_score(y_test, pred):.3f} 再現率 {recall_score(y_test, pred):.3f}")

import joblib

joblib.dump(final, "churn_model.joblib")          # 前処理ごと1ファイルに保存
loaded = joblib.load("churn_model.joblib")         # 本番環境で読み込む

new_customers = pd.DataFrame({
    "契約月数": [3, 60], "月額料金": [9800, np.nan], "プラン": ["ライト", "プレミアム"],
    "支払方法": ["コンビニ", "未登録"], "問い合わせ回数": [4, 0]})  # 欠損・未知の値も処理される
print(loaded.estimator_.predict_proba(new_customers)[:, 1].round(3))
print(loaded.predict(new_customers))

from scipy.stats import ks_2samp

# 翌月のデータ（値上げで月額料金が上がった想定）
next_month = X_test.copy()
next_month["月額料金"] = next_month["月額料金"] * 1.15

for col in ["契約月数", "月額料金"]:
    stat, p = ks_2samp(X_train[col].dropna(), next_month[col].dropna())
    print(f"{col}: p値={p:.4f}", "← 分布が変わった可能性" if p < 0.01 else "")
print("解約と判定した割合 学習時:", round(final.predict(X_train).mean(), 3),
      "翌月:", round(loaded.predict(next_month).mean(), 3))