consumer_behavior_mining.py

from __future__ import annotations

import json
import os
import sys
import warnings
import importlib.util
from pathlib import Path

import numpy as np
import pandas as pd
import joblib

os.environ.setdefault("LOKY_MAX_CPU_COUNT", "1")

from sklearn.base import BaseEstimator, ClassifierMixin, clone
from sklearn.decomposition import PCA
from sklearn.ensemble import HistGradientBoostingClassifier, RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (
    accuracy_score,
    calinski_harabasz_score,
    confusion_matrix,
    davies_bouldin_score,
    f1_score,
    precision_recall_fscore_support,
    precision_score,
    recall_score,
    roc_auc_score,
    silhouette_score,
)
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler


warnings.filterwarnings(
    "ignore",
    message="Could not find the number of physical cores.*",
    category=UserWarning,
)

RANDOM_STATE = 42
DATA_PATH = Path("marketing_campaign.csv")
OUTPUT_DIR = Path("outputs")
MIN_TOUCH_RESPONSE_PROBABILITY = 0.15
DISTURBANCE_PENALTY_WEIGHT = 3.0
CLUSTER_SIMILARITY = "cosine"
CLUSTER_LOG_AMOUNTS = True

PRIVACY_FIELDS = [
    "Income",
    "Education",
    "Marital_Status",
    "Kidhome",
    "Teenhome",
    "Year_Birth",
    "Dt_Customer",
]

DROP_FIELDS = ["ID", "Z_CostContact", "Z_Revenue"]

SPEND_COLS = [
    "MntWines",
    "MntFruits",
    "MntMeatProducts",
    "MntFishProducts",
    "MntSweetProducts",
    "MntGoldProds",
]

CATEGORY_NAMES = {
    "MntWines": "wine",
    "MntFruits": "fruit",
    "MntMeatProducts": "meat",
    "MntFishProducts": "fish",
    "MntSweetProducts": "sweet",
    "MntGoldProds": "gold",
}

PURCHASE_COLS = [
    "NumDealsPurchases",
    "NumWebPurchases",
    "NumCatalogPurchases",
    "NumStorePurchases",
]

CAMPAIGN_COLS = [
    "AcceptedCmp1",
    "AcceptedCmp2",
    "AcceptedCmp3",
    "AcceptedCmp4",
    "AcceptedCmp5",
]

BEHAVIOR_BASE_COLS = [
    "Recency",
    "NumDealsPurchases",
    "NumWebPurchases",
    "NumCatalogPurchases",
    "NumStorePurchases",
    "NumWebVisitsMonth",
    "AcceptedCmp1",
    "AcceptedCmp2",
    "AcceptedCmp3",
    "AcceptedCmp4",
    "AcceptedCmp5",
    "Complain",
]


class SimilarityKMeans:
    """K-Means variant using cosine or normalized Pearson similarity for assignment."""

    def __init__(
        self,
        n_clusters: int,
        similarity: str = "cosine",
        max_iter: int = 300,
        n_init: int = 20,
        random_state: int = RANDOM_STATE,
    ) -> None:
        if similarity not in {"cosine", "pearson"}:
            raise ValueError("similarity must be 'cosine' or 'pearson'")
        self.n_clusters = n_clusters
        self.similarity = similarity
        self.max_iter = max_iter
        self.n_init = n_init
        self.random_state = random_state
        self.scaler = StandardScaler()
        self.centers_: np.ndarray | None = None
        self.labels_: np.ndarray | None = None

    @staticmethod
    def _l2_normalize(x: np.ndarray) -> np.ndarray:
        norms = np.linalg.norm(x, axis=1, keepdims=True)
        return x / np.where(norms == 0, 1.0, norms)

    def _prepare_fit(self, x: pd.DataFrame) -> np.ndarray:
        scaled = self.scaler.fit_transform(x)
        if self.similarity == "pearson":
            scaled = scaled - scaled.mean(axis=1, keepdims=True)
        return self._l2_normalize(scaled)

    def _prepare_predict(self, x: pd.DataFrame) -> np.ndarray:
        scaled = self.scaler.transform(x)
        if self.similarity == "pearson":
            scaled = scaled - scaled.mean(axis=1, keepdims=True)
        return self._l2_normalize(scaled)

    def fit(self, x: pd.DataFrame) -> "SimilarityKMeans":
        prepared = self._prepare_fit(x)
        rng = np.random.default_rng(self.random_state)
        best_labels = None
        best_centers = None
        best_objective = -np.inf

        for _ in range(self.n_init):
            init_idx = rng.choice(len(prepared), size=self.n_clusters, replace=False)
            centers = prepared[init_idx].copy()
            labels = np.zeros(len(prepared), dtype=int)

            for _iteration in range(self.max_iter):
                similarities = prepared @ centers.T
                new_labels = similarities.argmax(axis=1)
                if np.array_equal(new_labels, labels):
                    break
                labels = new_labels
                for cluster_id in range(self.n_clusters):
                    members = prepared[labels == cluster_id]
                    if len(members) == 0:
                        centers[cluster_id] = prepared[rng.integers(0, len(prepared))]
                    else:
                        centers[cluster_id] = members.mean(axis=0)
                centers = self._l2_normalize(centers)

            objective = float(np.sum(np.max(prepared @ centers.T, axis=1)))
            if objective > best_objective:
                best_objective = objective
                best_labels = labels.copy()
                best_centers = centers.copy()

        assert best_labels is not None and best_centers is not None
        self.labels_ = best_labels
        self.centers_ = best_centers
        return self

    def fit_predict(self, x: pd.DataFrame) -> np.ndarray:
        return self.fit(x).labels_

    def predict(self, x: pd.DataFrame) -> np.ndarray:
        if self.centers_ is None:
            raise ValueError("SimilarityKMeans is not fitted.")
        prepared = self._prepare_predict(x)
        return (prepared @ self.centers_.T).argmax(axis=1)

    def transform_for_metrics(self, x: pd.DataFrame) -> np.ndarray:
        if self.centers_ is None:
            raise ValueError("SimilarityKMeans is not fitted.")
        return self._prepare_predict(x)


SimilarityKMeans.__module__ = "consumer_behavior_mining"
sys.modules.setdefault("consumer_behavior_mining", sys.modules[__name__])


class XGBoostBinaryWrapper(BaseEstimator, ClassifierMixin):
    def __init__(self, scale_pos_weight: float = 1.0) -> None:
        self.scale_pos_weight = scale_pos_weight
        self.model = None
        self.classes_ = np.array([0, 1])

    def fit(self, x: pd.DataFrame, y: pd.Series) -> "XGBoostBinaryWrapper":
        from xgboost import XGBClassifier

        self.model = XGBClassifier(
            objective="binary:logistic",
            n_estimators=500,
            max_depth=3,
            learning_rate=0.03,
            subsample=0.9,
            colsample_bytree=0.9,
            reg_lambda=1.0,
            eval_metric="logloss",
            scale_pos_weight=self.scale_pos_weight,
            random_state=RANDOM_STATE,
            n_jobs=1,
        )
        self.model.fit(x, y)
        return self

    def predict_proba(self, x: pd.DataFrame) -> np.ndarray:
        return self.model.predict_proba(x)


class LightGBMBinaryWrapper(BaseEstimator, ClassifierMixin):
    def __init__(self, scale_pos_weight: float = 1.0) -> None:
        self.scale_pos_weight = scale_pos_weight
        self.model = None
        self.classes_ = np.array([0, 1])

    def fit(self, x: pd.DataFrame, y: pd.Series) -> "LightGBMBinaryWrapper":
        from lightgbm import LGBMClassifier

        self.model = LGBMClassifier(
            objective="binary",
            n_estimators=500,
            max_depth=3,
            learning_rate=0.03,
            subsample=0.9,
            colsample_bytree=0.9,
            scale_pos_weight=self.scale_pos_weight,
            random_state=RANDOM_STATE,
            n_jobs=1,
            verbose=-1,
        )
        self.model.fit(x, y)
        return self

    def predict_proba(self, x: pd.DataFrame) -> np.ndarray:
        return self.model.predict_proba(x)


class CatBoostBinaryWrapper(BaseEstimator, ClassifierMixin):
    def __init__(self, scale_pos_weight: float = 1.0) -> None:
        self.scale_pos_weight = scale_pos_weight
        self.model = None
        self.classes_ = np.array([0, 1])

    def fit(self, x: pd.DataFrame, y: pd.Series) -> "CatBoostBinaryWrapper":
        from catboost import CatBoostClassifier

        self.model = CatBoostClassifier(
            iterations=500,
            depth=3,
            learning_rate=0.03,
            loss_function="Logloss",
            eval_metric="F1",
            class_weights=[1.0, self.scale_pos_weight],
            random_seed=RANDOM_STATE,
            verbose=False,
        )
        self.model.fit(x, y)
        return self

    def predict_proba(self, x: pd.DataFrame) -> np.ndarray:
        return self.model.predict_proba(x)


XGBoostBinaryWrapper.__module__ = "consumer_behavior_mining"
LightGBMBinaryWrapper.__module__ = "consumer_behavior_mining"
CatBoostBinaryWrapper.__module__ = "consumer_behavior_mining"


def safe_divide(numerator: pd.Series, denominator: pd.Series) -> pd.Series:
    return numerator / denominator.replace(0, np.nan)


def classify_operating_segment(
    total_spend: pd.Series | np.ndarray,
    response_probability: pd.Series | np.ndarray,
    spend_threshold: float,
    response_threshold: float,
) -> np.ndarray:
    high_spend = np.asarray(total_spend, dtype=float) >= float(spend_threshold)
    high_response = np.asarray(response_probability, dtype=float) >= float(response_threshold)
    return np.select(
        [
            high_spend & high_response,
            high_spend & ~high_response,
            ~high_spend & high_response,
            ~high_spend & ~high_response,
        ],
        [
            "高消费高响应型",
            "高消费低响应型",
            "低消费高响应型",
            "低消费低响应型",
        ],
        default="未分类",
    )


def load_data(path: Path) -> pd.DataFrame:
    if not path.exists():
        raise FileNotFoundError(f"Missing data file: {path.resolve()}")
    df = pd.read_csv(path, sep="\t")
    missing = set(SPEND_COLS + BEHAVIOR_BASE_COLS + ["Response"]) - set(df.columns)
    if missing:
        raise ValueError(f"Required columns are missing: {sorted(missing)}")
    return df


def build_behavior_features(
    df: pd.DataFrame,
    spend_frame: pd.DataFrame | None = None,
    log_amounts: bool = False,
) -> pd.DataFrame:
    raw_spend = df[SPEND_COLS].copy() if spend_frame is None else spend_frame[SPEND_COLS].copy()
    spend = np.log1p(raw_spend) if log_amounts else raw_spend
    base = df[BEHAVIOR_BASE_COLS].copy()

    features = pd.concat([base, spend], axis=1)
    total_spend = spend.sum(axis=1)
    raw_total_spend = raw_spend.sum(axis=1)
    total_purchases = df[PURCHASE_COLS].sum(axis=1)
    web_visits = df["NumWebVisitsMonth"]

    features["TotalSpend"] = np.log1p(raw_total_spend) if log_amounts else raw_total_spend
    features["TotalPurchases"] = total_purchases
    features["CampaignHistory"] = df[CAMPAIGN_COLS].sum(axis=1)
    features["DealSensitivity"] = safe_divide(df["NumDealsPurchases"], total_purchases).fillna(0)
    features["WebPurchaseRatio"] = safe_divide(df["NumWebPurchases"], total_purchases).fillna(0)
    features["CatalogPurchaseRatio"] = safe_divide(df["NumCatalogPurchases"], total_purchases).fillna(0)
    features["StorePurchaseRatio"] = safe_divide(df["NumStorePurchases"], total_purchases).fillna(0)
    features["VisitPurchaseRatio"] = safe_divide(df["NumWebPurchases"], web_visits).fillna(0)

    for col in SPEND_COLS:
        features[f"{CATEGORY_NAMES[col]}_share"] = safe_divide(spend[col], total_spend).fillna(0)

    features = features.replace([np.inf, -np.inf], 0).fillna(0)
    return features


def make_preference_dataset(df: pd.DataFrame) -> tuple[pd.DataFrame, pd.Series]:
    preferred_col = df[SPEND_COLS].idxmax(axis=1)
    masked_spend = df[SPEND_COLS].copy()
    for idx, col in preferred_col.items():
        masked_spend.at[idx, col] = 0
    x_pref = build_behavior_features(df, masked_spend)
    y_pref = preferred_col.map(CATEGORY_NAMES)
    return x_pref, y_pref


def choose_best_kmeans(x: pd.DataFrame) -> tuple[SimilarityKMeans, pd.DataFrame, dict]:
    best = None
    rows = []
    for k in range(2, 9):
        model = SimilarityKMeans(
            n_clusters=k,
            similarity=CLUSTER_SIMILARITY,
            n_init=30,
            random_state=RANDOM_STATE,
        )
        labels = model.fit_predict(x)
        metric_features = model.transform_for_metrics(x)
        metrics = {
            "k": k,
            "similarity": CLUSTER_SIMILARITY,
            "silhouette": silhouette_score(metric_features, labels, metric="cosine"),
            "davies_bouldin": davies_bouldin_score(metric_features, labels),
            "calinski_harabasz": calinski_harabasz_score(metric_features, labels),
        }
        rows.append(metrics)
        if best is None or metrics["silhouette"] > best[0]["silhouette"]:
            best = (metrics, model, labels)

    assert best is not None
    metrics, model, labels = best
    return model, pd.DataFrame(rows), {"labels": labels, **metrics}


def fit_preference_models(x: pd.DataFrame, y: pd.Series) -> tuple[dict, pd.DataFrame, pd.DataFrame]:
    x_train, x_test, y_train, y_test = train_test_split(
        x, y, test_size=0.2, random_state=RANDOM_STATE, stratify=y
    )
    models = {
        "logistic_regression": Pipeline(
            [
                ("scale", StandardScaler()),
                (
                    "model",
                    LogisticRegression(
                        max_iter=5000,
                        class_weight="balanced",
                        random_state=RANDOM_STATE,
                    ),
                ),
            ]
        ),
        "random_forest": RandomForestClassifier(
            n_estimators=500,
            min_samples_leaf=3,
            class_weight="balanced_subsample",
            random_state=RANDOM_STATE,
            n_jobs=1,
        ),
        "hist_gradient_boosting": HistGradientBoostingClassifier(
            learning_rate=0.06,
            max_iter=350,
            l2_regularization=0.1,
            random_state=RANDOM_STATE,
        ),
    }
    rows = []
    best = None
    predictions = pd.DataFrame(index=x_test.index)
    predictions["true_preference"] = y_test
    for name, model in models.items():
        model.fit(x_train, y_train)
        pred = model.predict(x_test)
        proba = model.predict_proba(x_test)
        confidence = proba.max(axis=1)
        row = {
            "model": name,
            "macro_f1": f1_score(y_test, pred, average="macro"),
            "weighted_f1": f1_score(y_test, pred, average="weighted"),
            "accuracy": accuracy_score(y_test, pred),
            "mean_confidence": float(np.mean(confidence)),
        }
        rows.append(row)
        if best is None or row["macro_f1"] > best["metrics"]["macro_f1"]:
            best = {"name": name, "model": model, "metrics": row, "x_test": x_test, "y_test": y_test}

    assert best is not None
    best_pred = best["model"].predict(x_test)
    best_proba = best["model"].predict_proba(x_test)
    predictions["predicted_preference"] = best_pred
    predictions["preference_confidence"] = best_proba.max(axis=1)
    return best, pd.DataFrame(rows).sort_values("macro_f1", ascending=False), predictions


def best_threshold_by_f1(y_true: pd.Series, proba: np.ndarray) -> tuple[float, float]:
    best_threshold = 0.5
    best_f1 = -1.0
    for threshold in np.linspace(0.05, 0.95, 181):
        pred = (proba >= threshold).astype(int)
        score = f1_score(y_true, pred, zero_division=0)
        if score > best_f1:
            best_threshold = float(threshold)
            best_f1 = float(score)
    return best_threshold, best_f1


def fit_response_models(x: pd.DataFrame, y: pd.Series) -> tuple[dict, pd.DataFrame, pd.DataFrame]:
    x_train_val, x_test, y_train_val, y_test = train_test_split(
        x, y, test_size=0.2, random_state=RANDOM_STATE, stratify=y
    )
    x_train, x_val, y_train, y_val = train_test_split(
        x_train_val, y_train_val, test_size=0.25, random_state=RANDOM_STATE, stratify=y_train_val
    )

    pos_weight = (len(y_train) - int(y_train.sum())) / max(int(y_train.sum()), 1)
    sample_weight = np.where(y_train == 1, pos_weight, 1.0)

    models = {
        "logistic_regression": Pipeline(
            [
                ("scale", StandardScaler()),
                (
                    "model",
                    LogisticRegression(
                        max_iter=5000,
                        class_weight="balanced",
                        random_state=RANDOM_STATE,
                    ),
                ),
            ]
        ),
        "random_forest": RandomForestClassifier(
            n_estimators=600,
            min_samples_leaf=3,
            class_weight="balanced_subsample",
            random_state=RANDOM_STATE,
            n_jobs=1,
        ),
        "hist_gradient_boosting": HistGradientBoostingClassifier(
            learning_rate=0.05,
            max_iter=450,
            l2_regularization=0.1,
            random_state=RANDOM_STATE,
        ),
    }
    if importlib.util.find_spec("xgboost") is not None:
        models["xgboost"] = XGBoostBinaryWrapper(scale_pos_weight=pos_weight)
    if importlib.util.find_spec("lightgbm") is not None:
        models["lightgbm"] = LightGBMBinaryWrapper(scale_pos_weight=pos_weight)
    if importlib.util.find_spec("catboost") is not None:
        models["catboost"] = CatBoostBinaryWrapper(scale_pos_weight=pos_weight)

    rows = []
    best = None
    for name, model in models.items():
        if name == "hist_gradient_boosting":
            model.fit(x_train, y_train, sample_weight=sample_weight)
        else:
            model.fit(x_train, y_train)
        val_proba = model.predict_proba(x_val)[:, 1]
        threshold, val_f1 = best_threshold_by_f1(y_val, val_proba)
        test_proba = model.predict_proba(x_test)[:, 1]
        test_pred = (test_proba >= threshold).astype(int)
        precision, recall, f1, _ = precision_recall_fscore_support(
            y_test, test_pred, average="binary", zero_division=0
        )
        row = {
            "model": name,
            "threshold": threshold,
            "validation_f1": val_f1,
            "test_f1": f1,
            "precision": precision,
            "recall": recall,
            "macro_f1": f1_score(y_test, test_pred, average="macro", zero_division=0),
            "roc_auc": roc_auc_score(y_test, test_proba),
        }
        rows.append(row)
        if best is None or row["test_f1"] > best["metrics"]["test_f1"]:
            best = {
                "name": name,
                "model": model,
                "metrics": row,
                "threshold": threshold,
                "x_test": x_test,
                "y_test": y_test,
                "proba": test_proba,
                "pred": test_pred,
            }

    assert best is not None
    response_predictions = pd.DataFrame(index=x_test.index)
    response_predictions["true_response"] = y_test
    response_predictions["response_probability"] = best["proba"]
    response_predictions["predicted_response"] = best["pred"]
    return best, pd.DataFrame(rows).sort_values("test_f1", ascending=False), response_predictions


def make_recommendations(
    df: pd.DataFrame,
    x_pref: pd.DataFrame,
    best_pref: dict,
    best_response: dict,
    response_predictions: pd.DataFrame,
) -> tuple[pd.DataFrame, pd.DataFrame]:
    test_idx = response_predictions.index
    pref_model = best_pref["model"]
    pref_pred = pref_model.predict(x_pref.loc[test_idx])
    pref_proba = pref_model.predict_proba(x_pref.loc[test_idx])
    pref_conf = pref_proba.max(axis=1)

    result = pd.DataFrame(index=test_idx)
    result["true_response"] = response_predictions["true_response"]
    result["response_probability"] = response_predictions["response_probability"]
    result["recommended_category"] = pref_pred
    result["preference_confidence"] = pref_conf
    result["total_spend"] = df.loc[test_idx, SPEND_COLS].sum(axis=1)
    result["recency"] = df.loc[test_idx, "Recency"]
    result["complain"] = df.loc[test_idx, "Complain"]
    result["campaign_history"] = df.loc[test_idx, CAMPAIGN_COLS].sum(axis=1)
    spend_threshold = float(df[SPEND_COLS].sum(axis=1).median())
    result["operating_segment"] = classify_operating_segment(
        result["total_spend"],
        result["response_probability"],
        spend_threshold,
        best_response["threshold"],
    )

    expected_spend = np.log1p(result["total_spend"])
    recency_risk = 1 - (result["recency"] / result["recency"].max()).fillna(0)
    low_response_history = (result["campaign_history"] == 0).astype(float)
    low_promo_interest = ((result["campaign_history"] == 0) & (df.loc[test_idx, "NumDealsPurchases"] == 0)).astype(float)
    disturbance_cost = (
        0.35 * recency_risk
        + 2.0 * result["complain"]
        + 0.25 * low_response_history
        + 0.40 * low_promo_interest
    )
    raw_score = result["response_probability"] * result["preference_confidence"] * expected_spend
    result["disturbance_cost"] = disturbance_cost
    result["recommendation_score"] = raw_score - DISTURBANCE_PENALTY_WEIGHT * disturbance_cost
    result["touch_allowed"] = (
        (result["response_probability"] >= MIN_TOUCH_RESPONSE_PROBABILITY)
        & (result["recommendation_score"] > 0)
        & (result["complain"] == 0)
    ).astype(int)

    strategies = {
        "low_disturbance_model": result["recommendation_score"].where(result["touch_allowed"] == 1, -np.inf),
        "baseline_total_spend": result["total_spend"],
        "baseline_campaign_history": result["campaign_history"],
        "baseline_random": pd.Series(
            np.random.default_rng(RANDOM_STATE).random(len(result)),
            index=result.index,
        ),
    }

    rows = []
    true_positive_count = int(result["true_response"].sum())
    for strategy_name, score in strategies.items():
        for k_ratio in [0.05, 0.10, 0.15, 0.20, 0.30]:
            k = max(1, int(round(len(result) * k_ratio)))
            chosen_idx = score.nlargest(k).index
            chosen = result.loc[chosen_idx]
            y_true = result["true_response"].to_numpy()
            y_pred = np.zeros(len(result), dtype=int)
            chosen_positions = result.index.get_indexer(chosen.index)
            y_pred[chosen_positions] = 1
            rows.append(
                {
                    "strategy": strategy_name,
                    "top_k_ratio": k_ratio,
                    "k": k,
                    "f1_at_k": f1_score(y_true, y_pred, zero_division=0),
                    "precision_at_k": precision_score(y_true, y_pred, zero_division=0),
                    "recall_at_k": recall_score(y_true, y_pred, zero_division=0),
                    "response_rate_at_k": float(chosen["true_response"].mean()),
                    "captured_responders": int(chosen["true_response"].sum()),
                    "total_responders": true_positive_count,
                }
            )

    result = result.sort_values("recommendation_score", ascending=False)
    return result, pd.DataFrame(rows)


def write_report(
    df: pd.DataFrame,
    cluster_metrics: pd.DataFrame,
    best_cluster: dict,
    preference_metrics: pd.DataFrame,
    response_metrics: pd.DataFrame,
    recommendation_metrics: pd.DataFrame,
    cluster_profile: pd.DataFrame,
    confusion: np.ndarray,
) -> None:
    best_pref = preference_metrics.iloc[0]
    best_resp = response_metrics.iloc[0]
    model_recommendation_metrics = recommendation_metrics[
        recommendation_metrics["strategy"] == "low_disturbance_model"
    ]
    best_reco = model_recommendation_metrics.sort_values("f1_at_k", ascending=False).iloc[0]

    report = f"""# 基于历史购买行为的消费者偏好识别与低打扰促销推荐

## 数据与隐私约束

- 数据集：Kaggle Customer Personality Analysis，文件 `marketing_campaign.csv`
- 样本量：{len(df)}
- 建模时不使用隐私字段：`{", ".join(PRIVACY_FIELDS)}`
- 删除无意义字段：`{", ".join(DROP_FIELDS)}`
- 核心输入仅来自历史消费、渠道购买、网站访问、促销响应历史和投诉记录。

## 任务设计

1. 消费习惯分群：用行为特征识别客户类型。
2. 品类偏好预测：隐藏客户最高消费品类后，预测其偏好品类。
3. 促销响应预测：预测 `Response=1`，主指标为 F1。
4. 低打扰推荐：综合响应概率、偏好置信度、历史消费和打扰成本，输出 Top-K 触达名单。

## 分群结果

- 最优聚类数：{int(best_cluster["k"])}
- Silhouette：{best_cluster["silhouette"]:.4f}
- Davies-Bouldin：{best_cluster["davies_bouldin"]:.4f}
- Calinski-Harabasz：{best_cluster["calinski_harabasz"]:.2f}

## 品类偏好预测

主指标优先使用 Macro-F1。

| 模型 | Macro-F1 | Weighted-F1 | Accuracy |
|---|---:|---:|---:|
"""
    for _, row in preference_metrics.iterrows():
        report += (
            f"| {row['model']} | {row['macro_f1']:.4f} | "
            f"{row['weighted_f1']:.4f} | {row['accuracy']:.4f} |\n"
        )

    report += f"""
最佳品类偏好模型：`{best_pref['model']}`，Macro-F1 = **{best_pref['macro_f1']:.4f}**。

## 促销响应预测

主指标优先使用 positive-class F1，即对真实响应客户 `Response=1` 的识别能力。

| 模型 | 阈值 | F1 | Precision | Recall | Macro-F1 | ROC-AUC |
|---|---:|---:|---:|---:|---:|---:|
"""
    for _, row in response_metrics.iterrows():
        report += (
            f"| {row['model']} | {row['threshold']:.3f} | {row['test_f1']:.4f} | "
            f"{row['precision']:.4f} | {row['recall']:.4f} | "
            f"{row['macro_f1']:.4f} | {row['roc_auc']:.4f} |\n"
        )

    report += f"""
最佳响应模型：`{best_resp['model']}`，F1 = **{best_resp['test_f1']:.4f}**。

响应模型混淆矩阵：

```text
{confusion.tolist()}
```

## 低打扰推荐评价

| 策略 | Top-K比例 | K | F1@K | Precision@K | Recall@K | Top-K响应率 | 捕获响应人数/总响应人数 |
|---|---:|---:|---:|---:|---:|---:|---:|
"""
    for _, row in recommendation_metrics.iterrows():
        report += (
            f"| {row['strategy']} | {row['top_k_ratio']:.0%} | {int(row['k'])} | {row['f1_at_k']:.4f} | "
            f"{row['precision_at_k']:.4f} | {row['recall_at_k']:.4f} | "
            f"{row['response_rate_at_k']:.4f} | {int(row['captured_responders'])}/{int(row['total_responders'])} |\n"
        )

    report += f"""
低打扰模型的最佳推荐切点：Top {best_reco['top_k_ratio']:.0%}，F1@K = **{best_reco['f1_at_k']:.4f}**。

## 输出文件

- `outputs/cluster_metrics.csv`：不同聚类数的分群指标。
- `outputs/cluster_profile.csv`：每个消费群体的行为画像。
- `outputs/preference_model_metrics.csv`：品类偏好预测指标。
- `outputs/response_model_metrics.csv`：促销响应预测指标。
- `outputs/recommendation_metrics.csv`：Top-K 推荐评价。
- `outputs/top_recommendations.csv`：推荐触达名单。
- `outputs/summary.json`：关键结果摘要。

## 结论

该方案避免使用收入、家庭构成、教育、婚姻、出生年份等隐私字段，仅用商家现实中更可能获得的历史购买行为进行建模。实验链路串行包含无监督分群、多分类偏好预测、二分类响应预测和约束推荐决策，并以 F1、Macro-F1、F1@K 作为优先评价指标。
"""

    OUTPUT_DIR.joinpath("experiment_report.md").write_text(report, encoding="utf-8")
    cluster_profile.to_csv(OUTPUT_DIR / "cluster_profile.csv", index=True, encoding="utf-8-sig")
    cluster_metrics.to_csv(OUTPUT_DIR / "cluster_metrics.csv", index=False, encoding="utf-8-sig")


def main() -> None:
    OUTPUT_DIR.mkdir(exist_ok=True)
    df = load_data(DATA_PATH)

    feature_df = build_behavior_features(df, log_amounts=CLUSTER_LOG_AMOUNTS)
    cluster_model, cluster_metrics, best_cluster = choose_best_kmeans(feature_df)
    labels = best_cluster["labels"]
    df_with_cluster = df.copy()
    df_with_cluster["cluster"] = labels

    cluster_profile = (
        df_with_cluster.groupby("cluster")[SPEND_COLS + PURCHASE_COLS + ["Recency", "NumWebVisitsMonth", "Complain", "Response"]]
        .mean()
        .round(3)
    )
    cluster_profile.insert(0, "sample_count", df_with_cluster.groupby("cluster").size())
    feature_df["BehaviorCluster"] = labels

    pca = PCA(n_components=2, random_state=RANDOM_STATE)
    pca_coords = pca.fit_transform(StandardScaler().fit_transform(feature_df.drop(columns=["BehaviorCluster"])))
    pd.DataFrame(
        {
            "pc1": pca_coords[:, 0],
            "pc2": pca_coords[:, 1],
            "cluster": labels,
            "response": df["Response"],
        }
    ).to_csv(OUTPUT_DIR / "cluster_pca_points.csv", index=False, encoding="utf-8-sig")

    x_pref, y_pref = make_preference_dataset(df)
    x_pref["BehaviorCluster"] = labels
    best_pref, preference_metrics, preference_predictions = fit_preference_models(x_pref, y_pref)

    x_response = feature_df.copy()
    y_response = df["Response"].astype(int)
    best_response, response_metrics, response_predictions = fit_response_models(x_response, y_response)

    final_preference_model = clone(best_pref["model"])
    final_preference_model.fit(x_pref, y_pref)

    full_pos_weight = (len(y_response) - int(y_response.sum())) / max(int(y_response.sum()), 1)
    if best_response["name"] == "xgboost":
        final_response_model = XGBoostBinaryWrapper(scale_pos_weight=full_pos_weight)
    elif best_response["name"] == "lightgbm":
        final_response_model = LightGBMBinaryWrapper(scale_pos_weight=full_pos_weight)
    elif best_response["name"] == "catboost":
        final_response_model = CatBoostBinaryWrapper(scale_pos_weight=full_pos_weight)
    else:
        final_response_model = clone(best_response["model"])
    if best_response["name"] == "hist_gradient_boosting":
        response_sample_weight = np.where(y_response == 1, full_pos_weight, 1.0)
        final_response_model.fit(x_response, y_response, sample_weight=response_sample_weight)
    else:
        final_response_model.fit(x_response, y_response)

    model_bundle = {
        "cluster_model": cluster_model,
        "preference_model": final_preference_model,
        "response_model": final_response_model,
        "response_threshold": float(best_response["threshold"]),
        "behavior_feature_columns": list(feature_df.drop(columns=["BehaviorCluster"]).columns),
        "response_feature_columns": list(x_response.columns),
        "preference_feature_columns": list(x_pref.columns),
        "spend_columns": SPEND_COLS,
        "purchase_columns": PURCHASE_COLS,
        "campaign_columns": CAMPAIGN_COLS,
        "total_spend_threshold": float(df[SPEND_COLS].sum(axis=1).median()),
        "privacy_fields_excluded": PRIVACY_FIELDS,
    }
    joblib.dump(model_bundle, OUTPUT_DIR / "model_bundle.joblib")

    recommendations, recommendation_metrics = make_recommendations(
        df, x_pref, best_pref, best_response, response_predictions
    )

    preference_metrics.to_csv(OUTPUT_DIR / "preference_model_metrics.csv", index=False, encoding="utf-8-sig")
    response_metrics.to_csv(OUTPUT_DIR / "response_model_metrics.csv", index=False, encoding="utf-8-sig")
    recommendation_metrics = recommendation_metrics.sort_values(
        ["strategy", "top_k_ratio"], ascending=[True, True]
    )
    recommendation_metrics.to_csv(OUTPUT_DIR / "recommendation_metrics.csv", index=False, encoding="utf-8-sig")
    preference_predictions.to_csv(OUTPUT_DIR / "preference_predictions.csv", index=True, encoding="utf-8-sig")
    response_predictions.to_csv(OUTPUT_DIR / "response_predictions.csv", index=True, encoding="utf-8-sig")
    recommendations.head(100).to_csv(OUTPUT_DIR / "top_recommendations.csv", index=True, encoding="utf-8-sig")

    confusion = confusion_matrix(best_response["y_test"], best_response["pred"])
    summary = {
        "sample_count": int(len(df)),
        "privacy_fields_excluded": PRIVACY_FIELDS,
        "best_cluster_k": int(best_cluster["k"]),
        "cluster_similarity": str(best_cluster["similarity"]),
        "cluster_log_amounts": CLUSTER_LOG_AMOUNTS,
        "best_cluster_silhouette": float(best_cluster["silhouette"]),
        "best_preference_model": str(preference_metrics.iloc[0]["model"]),
        "best_preference_macro_f1": float(preference_metrics.iloc[0]["macro_f1"]),
        "best_response_model": str(response_metrics.iloc[0]["model"]),
        "best_response_f1": float(response_metrics.iloc[0]["test_f1"]),
        "best_response_threshold": float(response_metrics.iloc[0]["threshold"]),
        "best_recommendation_f1_at_k": float(
            recommendation_metrics.loc[
                recommendation_metrics["strategy"] == "low_disturbance_model", "f1_at_k"
            ].max()
        ),
    }
    (OUTPUT_DIR / "summary.json").write_text(json.dumps(summary, ensure_ascii=False, indent=2), encoding="utf-8")

    write_report(
        df=df,
        cluster_metrics=cluster_metrics,
        best_cluster=best_cluster,
        preference_metrics=preference_metrics,
        response_metrics=response_metrics,
        recommendation_metrics=recommendation_metrics,
        cluster_profile=cluster_profile,
        confusion=confusion,
    )

    print(json.dumps(summary, ensure_ascii=False, indent=2))


if __name__ == "__main__":
    main()