Segmentación RFM y CLV predictivo

"""#8 Segmentación RFM y valor futuro del cliente (CLV predictivo).

Datos: Online Retail (UCI), 541,909 transacciones de un minorista en línea (dic-2010 a dic-2011).
1) Segmentación RFM con K-Means.  2) CLV: con 9 meses de historia se predice el gasto de los
3 meses siguientes y se valida contra lo que realmente ocurrió.
"""
import io

import numpy as np
import pandas as pd
from lightgbm import LGBMRegressor
from lifetimes import BetaGeoFitter, GammaGammaFitter
from lifetimes.utils import ConvergenceError, summary_data_from_transaction_data
from scipy.stats import spearmanr
from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score, davies_bouldin_score, silhouette_score
from sklearn.mixture import GaussianMixture
from sklearn.model_selection import KFold
from sklearn.preprocessing import StandardScaler

from common import SEED, bars, download, export, line, read_from_zip, table
from validacion import model_card

URL = "https://archive.ics.uci.edu/static/public/352/online+retail.zip"
CORTE = pd.Timestamp("2011-09-01")   # calibración: antes del corte · validación: 3 meses después


def load() -> pd.DataFrame:
    cache = download(URL, "online_retail.zip").with_suffix(".parquet")
    if cache.exists():
        return pd.read_parquet(cache)
    df = pd.read_excel(io.BytesIO(read_from_zip(cache.with_suffix(".zip"), ".xlsx")))
    df = df.dropna(subset=["CustomerID"])
    df = df[(df.Quantity > 0) & (df.UnitPrice > 0)]
    df["monto"] = df.Quantity * df.UnitPrice
    df = df[["InvoiceNo", "InvoiceDate", "CustomerID", "Country", "monto"]]
    df.to_parquet(cache)
    return df


def rfm(tx: pd.DataFrame, fecha_ref: pd.Timestamp) -> pd.DataFrame:
    f = tx.groupby("CustomerID").agg(ultima=("InvoiceDate", "max"), primera=("InvoiceDate", "min"),
                                     frecuencia=("InvoiceNo", "nunique"), monetario=("monto", "sum"))
    f["recencia"] = (fecha_ref - f.ultima).dt.days
    f["antiguedad"] = (fecha_ref - f.primera).dt.days
    f["ticket"] = f.monetario / f.frecuencia
    return f.drop(columns=["ultima", "primera"])


def estabilidad(Z: np.ndarray, k: int, n=20) -> tuple[float, float]:
    """ARI medio entre la segmentación completa y la de submuestras del 80 % (sobre las filas comunes)."""
    rng = np.random.default_rng(SEED)
    base = KMeans(k, n_init=10, random_state=SEED).fit_predict(Z)
    aris = []
    for i in range(n):
        idx = rng.choice(len(Z), int(len(Z) * 0.8), replace=False)
        lab = KMeans(k, n_init=5, random_state=i).fit_predict(Z[idx])
        aris.append(adjusted_rand_score(base[idx], lab))
    return float(np.mean(aris)), float(np.std(aris))


def clv_bgnbd(cal: pd.DataFrame, dias: int) -> pd.Series:
    """BG/NBD (frecuencia de compra) × Gamma-Gamma (valor por compra): el estándar probabilístico de CLV."""
    d = cal.assign(dia=cal.InvoiceDate.dt.normalize()).groupby(["CustomerID", "dia"]).monto.sum().reset_index()
    # Unidad semanal: con días el optimizador no converge (muchas compras el mismo día de la semana)
    s = summary_data_from_transaction_data(d, "CustomerID", "dia", monetary_value_col="monto",
                                           observation_period_end=CORTE - pd.Timedelta(days=1), freq="W")
    try:
        bgf = BetaGeoFitter(penalizer_coef=0.0).fit(s.frequency, s.recency, s["T"])
    except ConvergenceError:
        bgf = BetaGeoFitter(penalizer_coef=0.1).fit(s.frequency, s.recency, s["T"])
    compras = bgf.conditional_expected_number_of_purchases_up_to_time(dias / 7, s.frequency, s.recency, s["T"])
    rep = s[s.frequency > 0]
    ggf = GammaGammaFitter(penalizer_coef=0.01).fit(rep.frequency, rep.monetary_value)
    valor = ggf.conditional_expected_average_profit(s.frequency, s.monetary_value)
    valor[s.frequency == 0] = rep.monetary_value.mean()
    return (compras * valor).rename("bgnbd")


def run():
    tx = load()
    ref = tx.InvoiceDate.max().normalize() + pd.Timedelta(days=1)
    r = rfm(tx, ref)

    # ---------- 1) Segmentación ----------
    L = np.log1p(r[["recencia", "frecuencia", "monetario"]])
    Z = StandardScaler().fit_transform(L)
    rng = np.random.default_rng(SEED)
    muestra = rng.choice(len(Z), min(3000, len(Z)), replace=False)
    diag = []
    for k in range(2, 9):
        km_k = KMeans(k, n_init=10, random_state=SEED).fit(Z)
        diag.append({"k": k, "inercia": km_k.inertia_,
                     "silueta": silhouette_score(Z[muestra], km_k.labels_[muestra]),
                     "davies_bouldin": davies_bouldin_score(Z, km_k.labels_)})
    diag = pd.DataFrame(diag)
    K = 4   # elección de negocio (4 segmentos accionables) validada abajo con silueta y estabilidad
    km = KMeans(K, n_init=10, random_state=SEED).fit(Z)
    r["cluster"] = km.labels_
    ari_m, ari_s = estabilidad(Z, K)
    gmm = GaussianMixture(K, random_state=SEED, n_init=3).fit(Z)
    ari_gmm = adjusted_rand_score(km.labels_, gmm.predict(Z))
    # Regla de negocio clásica: quintiles RFM → 4 grupos por puntaje total
    score_rfm = (pd.qcut(-r.recencia, 5, labels=False, duplicates="drop") + pd.qcut(r.frecuencia.rank(method="first"), 5, labels=False)
                 + pd.qcut(r.monetario, 5, labels=False))
    ari_regla = adjusted_rand_score(km.labels_, pd.qcut(score_rfm.rank(method="first"), K, labels=False))

    perfil = r.groupby("cluster").agg(clientes=("recencia", "size"), recencia=("recencia", "median"),
                                      frecuencia=("frecuencia", "median"), monetario=("monetario", "median"),
                                      ingreso=("monetario", "sum"))
    nombres = {perfil.monetario.idxmax(): "Campeones", perfil.recencia.idxmax(): "Dormidos"}
    resto = perfil.drop(index=list(nombres)).sort_values("monetario", ascending=False).index
    nombres.update(dict(zip(resto, ["Leales", "Nuevos / prometedores"])))
    perfil.index = perfil.index.map(nombres)
    perfil["% ingreso"] = (perfil.ingreso / perfil.ingreso.sum() * 100).round(1)
    perfil["% clientes"] = (perfil.clientes / perfil.clientes.sum() * 100).round(1)
    perfil = perfil.sort_values("monetario", ascending=False)
    # Perfil estandarizado (z) de cada segmento para la gráfica de coordenadas paralelas
    zc = pd.DataFrame(Z, columns=["Recencia", "Frecuencia", "Monetario"]).groupby(r["cluster"].map(nombres).values).mean()

    # ---------- 2) CLV: 9 meses de historia → gasto de los 3 meses siguientes ----------
    cal, val = tx[tx.InvoiceDate < CORTE], tx[tx.InvoiceDate >= CORTE]
    dias = (tx.InvoiceDate.max().normalize() - CORTE).days + 1
    Xc = rfm(cal, CORTE)
    Xc["meses_activos"] = cal.groupby("CustomerID").InvoiceDate.apply(lambda s: s.dt.to_period("M").nunique())
    Xc["reino_unido"] = (cal.groupby("CustomerID").Country.first() == "United Kingdom").astype(int)
    real = val.groupby("CustomerID").monto.sum().reindex(Xc.index, fill_value=0).values
    # (a) regla simple, (b) BG/NBD + Gamma-Gamma, (c) LightGBM con predicción fuera de fold
    pred = {"Regla: promedio histórico": (Xc.monetario / np.maximum(Xc.antiguedad, 90) * dias).values,
            "BG/NBD + Gamma-Gamma": clv_bgnbd(cal, dias).reindex(Xc.index).fillna(0).values}
    oof = np.zeros(len(Xc))
    for tr, te in KFold(5, shuffle=True, random_state=SEED).split(Xc):
        m = LGBMRegressor(n_estimators=400, learning_rate=0.03, num_leaves=15, min_child_samples=30,
                          objective="regression_l1", random_state=SEED, verbose=-1)
        oof[te] = m.fit(Xc.iloc[tr], real[tr]).predict(Xc.iloc[te])
    pred["LightGBM (CV 5-fold)"] = oof

    def concentracion(score):
        top = np.argsort(-score)[: int(len(score) * 0.2)]
        return real[top].sum() / real.sum()
    def mae_ci(e):
        b = [np.mean(np.abs(e[i])) for i in (rng.integers(0, len(e), len(e)) for _ in range(300))]
        return np.mean(np.abs(e)), *np.percentile(b, [2.5, 97.5])
    evals = []
    for nombre, pr in pred.items():
        mae, lo, hi = mae_ci(real - pr)
        evals.append({"Modelo": nombre, "MAE £": f"{mae:,.0f} [{lo:,.0f}–{hi:,.0f}]", "mae": mae,
                      "Spearman": round(spearmanr(real, pr)[0], 3), "Top 20 % concentra": f"{concentracion(pr):.0%}",
                      "Gasto total predicho vs. real": f"{pr.sum() / real.sum():.0%}"})
    evals = pd.DataFrame(evals)
    mejor = evals.loc[evals.mae.idxmin(), "Modelo"]
    dec = pd.DataFrame({"pred": pred[mejor], "real": real})
    dec["decil"] = pd.qcut(-dec.pred.rank(method="first"), 10, labels=[f"D{i}" for i in range(1, 11)])
    dec = dec.groupby("decil", observed=True).agg(pred=("pred", "mean"), real=("real", "mean")).reset_index()

    perfil_tbl = perfil.reset_index().rename(columns={"cluster": "Segmento"})[
        ["Segmento", "clientes", "% clientes", "recencia", "frecuencia", "monetario", "% ingreso"]]
    perfil_tbl.columns = ["Segmento", "Clientes", "% clientes", "Recencia (días, mediana)",
                          "Compras (mediana)", "Gasto (mediana, £)", "% del ingreso"]
    perfil_tbl["Gasto (mediana, £)"] = perfil_tbl["Gasto (mediana, £)"].round(0)
    diag_tbl = diag.round(3).rename(columns={"inercia": "Inercia", "silueta": "Silueta (↑)",
                                             "davies_bouldin": "Davies-Bouldin (↓)"})
    camp = perfil.iloc[0]
    sil_k = float(diag.set_index("k").loc[K, "silueta"])
    k_sil = int(diag.loc[diag.silueta.idxmax(), "k"])
    e = evals.set_index("Modelo")
    result = {
        "title": "Segmentación RFM y CLV predictivo",
        "summary": [
            f"{len(r):,} clientes en {K} segmentos. Validación: silueta {sil_k:.2f} (la máxima es k = {k_sil}, pero k = {K} "
            f"es más accionable), estabilidad bootstrap ARI {ari_m:.2f} ± {ari_s:.2f} y acuerdo con un modelo de mezclas "
            f"gaussianas ARI {ari_gmm:.2f}.",
            f"Los \"{perfil.index[0]}\" son {camp['% clientes']:.0f} % de los clientes y generan {camp['% ingreso']:.0f} % del "
            f"ingreso. Acuerdo con la regla clásica de quintiles RFM: ARI {ari_regla:.2f}.",
            "CLV a 3 meses fuera de muestra — MAE: " + ", ".join(f"{k} £{v:,.0f}" for k, v in e.mae.items()) +
            f". Los IC bootstrap se traslapan (empate práctico) y todos subestiman el gasto total "
            f"(temporada navideña fuera de la historia): {', '.join(f'{k.split(' (')[0]} {v}' for k, v in e['Gasto total predicho vs. real'].items())}.",
            "Cambio de enfoque: además de ML, un modelo probabilístico (BG/NBD + Gamma-Gamma) que se explica con dos "
            "ideas de negocio —qué tan seguido compra y cuánto gasta por compra— y que no necesita etiquetas.",
        ],
        "metrics": [
            {"label": "Silueta k = 4", "value": f"{sil_k:.2f}"},
            {"label": "Estabilidad (ARI)", "value": f"{ari_m:.2f} ± {ari_s:.2f}"},
            {"label": "Mejor CLV", "value": mejor.split(" (")[0]},
            {"label": "MAE mejor CLV", "value": f"£{e.loc[mejor, 'mae']:,.0f}"},
        ],
        "charts": [
            line("perfil", "Perfil de cada segmento (z-score de log RFM)",
                 {seg: [[c, round(float(v), 3)] for c, v in fila.items()] for seg, fila in zc.iterrows()},
                 x_label="Dimensión", y_label="Desviaciones estándar vs. promedio", fmt="dec",
                 subtitle="Recencia alta = hace mucho que no compra"),
            bars("ingreso", "% del ingreso por segmento", perfil.index, perfil["% ingreso"], fmt="pct"),
            line("sil", "Elección de k: silueta (más alto es mejor)", {"Silueta": diag[["k", "silueta"]].round(3).values.tolist()},
                 x_label="Número de clusters", y_label="Silueta", fmt="dec", markers=[[K, "k elegido"]]),
            line("codo", "Método del codo", {"Inercia": diag[["k", "inercia"]].round(1).values.tolist()},
                 x_label="Número de clusters", y_label="Inercia", markers=[[K, "k elegido"]]),
            bars("clv", "CLV: error absoluto medio por modelo (£, menor es mejor)", evals.Modelo, evals.mae.round(0).values,
                 fmt="gbp"),
            line("dec", f"CLV ({mejor}): gasto predicho vs. real por decil",
                 {"Predicho": [[d, round(float(v), 1)] for d, v in zip(dec.decil.astype(str), dec.pred)],
                  "Real": [[d, round(float(v), 1)] for d, v in zip(dec.decil.astype(str), dec.real)]},
                 x_label="Decil de CLV predicho", y_label="£ en los 3 meses siguientes", fmt="gbp"),
        ],
        "tables": [table("Perfil de segmentos", perfil_tbl),
                   table("Comparación de modelos de CLV (validación temporal)", evals.drop(columns="mae")),
                   table("Diagnóstico de número de clusters", diag_tbl)],
        "card": model_card(
            "K-Means sobre log-RFM estandarizado; CLV con regla simple, BG/NBD + Gamma-Gamma y LightGBM.",
            ["Clustering: silueta, Davies-Bouldin, codo, estabilidad bootstrap (ARI) y acuerdo con GMM y regla RFM.",
             "CLV: validación temporal (historia hasta ago-2011, objetivo sep–dic 2011); LightGBM con predicción "
             "fuera de fold; MAE con IC bootstrap, Spearman y concentración del top 20 %."],
            ["Las devoluciones y facturas sin cliente se excluyen.",
             "BG/NBD supone que un cliente \"muere\" una sola vez y compra como proceso de Poisson mientras vive."],
            ["Un solo año de historia: el trimestre objetivo incluye la temporada navideña.",
             "Clientes mayoristas extremos dominan el error absoluto; para decisiones conviene ver el orden (Spearman)."]),
    }
    sample = r[["recencia", "frecuencia", "monetario"]].round(2).sample(1500, random_state=1)
    export("p08-segmentacion-clv", result, sample)
    return result


if __name__ == "__main__":
    run()