Customer Segmentation with RFM and K-Means

import pandas as pd
import numpy as np
from scipy import stats

np.random.seed(42)

# ── Customer RFM table ────────────────────────────────────────────────────────
df = pd.DataFrame({
    "customer_id":       [f"C{i:03d}" for i in range(1, 21)],
    "recency_days":      [8,142,3,287,19,64,5,198,31,112,7,243,25,88,12,172,44,320,16,55],
    "frequency":         [14,2,22,1,9,4,18,2,7,3,20,1,11,5,16,2,6,1,13,8],
    "monetary":          [2840,184,5120,67,1640,420,4280,145,1180,310,4890,88,2140,580,3640,165,760,52,2560,1020],
    "avg_order_value":   [202.9,92.0,232.7,67.0,182.2,105.0,237.8,72.5,168.6,103.3,244.5,88.0,194.5,116.0,227.5,82.5,126.7,52.0,196.9,127.5],
    "first_purchase_days":[312,180,348,290,280,145,356,210,268,195,340,300,288,162,326,220,198,310,274,232]
})

rfm_cols = ["recency_days","frequency","monetary"]

# ── Distribution of raw RFM ───────────────────────────────────────────────────
print("Raw RFM feature statistics:")
print(df[rfm_cols].describe().round(1).to_string())

# ── Min-Max normalisation ─────────────────────────────────────────────────────
def min_max_scale(series):
    """Scale series to [0, 1]."""
    mn, mx = series.min(), series.max()
    return (series - mn) / (mx - mn)

df["R_norm"] = 1 - min_max_scale(df["recency_days"])  # invert: low days = good
df["F_norm"] =     min_max_scale(df["frequency"])
df["M_norm"] =     min_max_scale(df["monetary"])

# RFM composite score (equal weights) — useful for quick ranking
df["rfm_score"] = (df["R_norm"] + df["F_norm"] + df["M_norm"]).round(4)

print("\nNormalised RFM features (first 8 rows):")
print(df[["customer_id","recency_days","frequency","monetary",
          "R_norm","F_norm","M_norm","rfm_score"]].head(8).round(3).to_string(index=False))

print(f"\nTop 3 customers by RFM composite score:")
print(df.nlargest(3,"rfm_score")[["customer_id","recency_days","frequency","monetary","rfm_score"]].to_string(index=False))
print(f"\nBottom 3 customers by RFM composite score:")
print(df.nsmallest(3,"rfm_score")[["customer_id","recency_days","frequency","monetary","rfm_score"]].to_string(index=False))

X = df[["R_norm","F_norm","M_norm"]].values   # feature matrix: (20, 3)

def kmeans(X, k, max_iter=100, n_init=10):
    """
    K-means clustering from scratch.
    Runs n_init times with different random seeds, returns best result.
    """
    best_wcss    = np.inf
    best_labels  = None
    best_centers = None

    for seed in range(n_init):
        rng = np.random.RandomState(seed)
        # 1. Random initialisation: pick k rows as starting centroids
        idx     = rng.choice(len(X), size=k, replace=False)
        centers = X[idx].copy()

        for _ in range(max_iter):
            # 2. Assignment: each point → nearest centroid (Euclidean distance)
            dists  = np.sqrt(((X[:, None] - centers[None]) ** 2).sum(axis=2))
            labels = dists.argmin(axis=1)

            # 3. Update: recompute centroids as cluster means
            new_centers = np.array([
                X[labels == j].mean(axis=0) if (labels == j).any() else centers[j]
                for j in range(k)
            ])

            # 4. Convergence check: stop if centroids didn't move
            if np.allclose(centers, new_centers, atol=1e-6):
                break
            centers = new_centers

        # WCSS: sum of squared distances to assigned centroid
        wcss = sum(
            ((X[labels == j] - centers[j]) ** 2).sum()
            for j in range(k)
            if (labels == j).any()
        )
        if wcss < best_wcss:
            best_wcss, best_labels, best_centers = wcss, labels, centers

    return best_labels, best_centers, best_wcss

# ── Elbow method: WCSS for K = 2 … 7 ─────────────────────────────────────────
wcss_values = {}
for k in range(2, 8):
    _, _, wcss = kmeans(X, k=k)
    wcss_values[k] = round(wcss, 4)

print("Elbow method — WCSS by number of clusters:")
for k, wcss in wcss_values.items():
    bar = "█" * int(wcss * 20)
    print(f"  K={k}  WCSS={wcss:.4f}  {bar}")

# Marginal improvement: drop in WCSS from K to K+1
print("\nMarginal WCSS improvement (K → K+1):")
ks = sorted(wcss_values.keys())
for i in range(len(ks)-1):
    delta = wcss_values[ks[i]] - wcss_values[ks[i+1]]
    pct   = delta / wcss_values[ks[i]] * 100
    print(f"  K={ks[i]} → K={ks[i+1]}: WCSS drops {delta:.4f} ({pct:.1f}%)")

print("\nOptimal K = 4 (elbow at largest marginal improvement drop)")

# ── Fit K=4 ──────────────────────────────────────────────────────────────────
labels, centers, wcss = kmeans(X, k=4, n_init=20)
df["cluster"] = labels

# ── Raw RFM profile per cluster ───────────────────────────────────────────────
profile = df.groupby("cluster").agg(
    customers        = ("customer_id",    "count"),
    mean_recency     = ("recency_days",   "mean"),
    mean_frequency   = ("frequency",      "mean"),
    mean_monetary    = ("monetary",       "mean"),
    total_monetary   = ("monetary",       "sum"),
    mean_rfm_score   = ("rfm_score",      "mean"),
    mean_aov         = ("avg_order_value","mean")
).round(1).reset_index()

profile = profile.sort_values("mean_rfm_score", ascending=False).reset_index(drop=True)

print("K=4 cluster profiles — ranked by mean RFM score:")
print(profile.to_string(index=False))

# ── Assign persona names based on profile ────────────────────────────────────
def assign_persona(row):
    if row["mean_rfm_score"] >= 2.2:  return "Champions"
    if row["mean_rfm_score"] >= 1.4:  return "Loyal Buyers"
    if row["mean_rfm_score"] >= 0.6:  return "At-Risk"
    return "Hibernating"

profile["persona"] = profile.apply(assign_persona, axis=1)
cluster_to_persona = dict(zip(profile["cluster"], profile["persona"]))
df["persona"] = df["cluster"].map(cluster_to_persona)

print("\nPersona assignment:")
print(profile[["cluster","persona","customers","mean_recency","mean_frequency",
               "mean_monetary","total_monetary"]].to_string(index=False))

# Revenue at risk: Champions + Loyal Buyers
high_value_revenue = profile[profile["persona"].isin(["Champions","Loyal Buyers"])]["total_monetary"].sum()
print(f"\nRevenue from Champions + Loyal Buyers: ${high_value_revenue:,.0f}")
print(f"Share of total portfolio revenue: {high_value_revenue/df['monetary'].sum()*100:.1f}%")

def silhouette_score_manual(X, labels):
    """
    Manual silhouette score computation.
    For each point: a = mean dist to same-cluster points
                    b = min mean dist to any other cluster's points
    s = (b - a) / max(a, b)
    """
    n      = len(X)
    scores = []
    unique = np.unique(labels)

    for i in range(n):
        same_cluster  = X[labels == labels[i]]
        # a: mean intra-cluster distance (excluding self)
        other_same = same_cluster[np.arange(len(same_cluster)) != np.where(
            np.all(same_cluster == X[i], axis=1))[0][0]]
        a = np.sqrt(((other_same - X[i]) ** 2).sum(axis=1)).mean() if len(other_same) > 0 else 0

        # b: min mean distance to any other cluster
        b_vals = []
        for c in unique:
            if c == labels[i]:
                continue
            other_cluster = X[labels == c]
            b_vals.append(np.sqrt(((other_cluster - X[i]) ** 2).sum(axis=1)).mean())
        b = min(b_vals) if b_vals else 0

        s = (b - a) / max(a, b) if max(a, b) > 0 else 0
        scores.append(s)

    return np.mean(scores), np.array(scores)

mean_sil, sil_scores = silhouette_score_manual(X, labels)
df["silhouette"] = sil_scores.round(4)

print(f"Mean silhouette score (K=4): {mean_sil:.4f}")
print(f"Interpretation: {'Good separation (>0.5)' if mean_sil>0.5 else 'Moderate separation (0.25-0.5)' if mean_sil>0.25 else 'Weak separation (<0.25)'}")

print("\nSilhouette score by persona:")
print(df.groupby("persona")["silhouette"].agg(["mean","min","max"]).round(3).to_string())

# Stability: run 10 seeds, measure agreement on Champions vs non-Champions
stability_results = []
for seed in range(10):
    lbl, _, _ = kmeans(X, k=4, n_init=1)
    # map to persona using rfm_score of cluster centroids
    cl_scores = {c: X[lbl==c].mean(axis=0).sum() for c in range(4)}
    top_cluster = max(cl_scores, key=cl_scores.get)
    stability_results.append(set(np.where(lbl == top_cluster)[0]))

# Champions overlap across seeds
overlap = stability_results[0]
for s in stability_results[1:]:
    overlap = overlap & s
stability_pct = len(overlap) / 5 * 100   # 5 = expected Champions cluster size

print(f"\nCluster stability (10 seeds):")
print(f"  Customers consistently assigned to Champions cluster: {len(overlap)}/5")
print(f"  Stability: {stability_pct:.0f}%")

# ── Campaign ROI estimation ───────────────────────────────────────────────────
CAMPAIGN_COST_PER_CUSTOMER = 4.50   # email + offer cost per customer

campaign_params = {
    "Champions":   {"target_uplift": 0.08, "avg_order": 230, "desc": "Loyalty reward — VIP tier upgrade"},
    "Loyal Buyers":{"target_uplift": 0.12, "avg_order": 195, "desc": "Cross-sell — complementary categories"},
    "At-Risk":     {"target_uplift": 0.20, "avg_order": 100, "desc": "Win-back — 20% discount code"},
    "Hibernating": {"target_uplift": 0.08, "avg_order":  70, "desc": "Reactivation — heavy incentive offer"}
}

persona_order = ["Champions","Loyal Buyers","At-Risk","Hibernating"]
print("Campaign brief by persona:\n")
for persona in persona_order:
    grp    = df[df["persona"]==persona]
    params = campaign_params[persona]
    n      = len(grp)
    uplift_customers = int(n * params["target_uplift"] * 10)  # scale to full base
    revenue_uplift   = uplift_customers * params["avg_order"]
    campaign_cost    = n * CAMPAIGN_COST_PER_CUSTOMER * 10
    net_roi          = (revenue_uplift - campaign_cost) / campaign_cost * 100

    print(f"  ── {persona.upper()} ({'%.0f' % n} sample customers) ──")
    print(f"     Mean recency:    {grp['recency_days'].mean():.0f} days")
    print(f"     Mean frequency:  {grp['frequency'].mean():.1f} orders")
    print(f"     Mean monetary:   ${grp['monetary'].mean():,.0f}")
    print(f"     Campaign:        {params['desc']}")
    print(f"     Est. uplift:     {params['target_uplift']:.0%} conversion × {uplift_customers} customers")
    print(f"     Revenue uplift:  ${revenue_uplift:,}")
    print(f"     Campaign cost:   ${campaign_cost:.0f}")
    print(f"     Campaign ROI:    {net_roi:.0f}%")
    print()

# ── Customer list for CRM export ──────────────────────────────────────────────
export = df[["customer_id","persona","recency_days","frequency","monetary",
             "rfm_score","silhouette"]].sort_values("rfm_score", ascending=False)
print("CRM export — full customer list with persona assignments:")
print(export.to_string(index=False))

# ── ANOVA: which RFM feature differs most across clusters? ───────────────────
print("ANOVA F-statistic — feature variation across clusters:")
for feat in rfm_cols:
    groups = [df[df["cluster"]==c][feat].values for c in df["cluster"].unique()]
    f_stat, p_val = stats.f_oneway(*groups)
    sig = "***" if p_val < 0.001 else "**" if p_val < 0.01 else "*" if p_val < 0.05 else "ns"
    print(f"  {feat:<18}  F={f_stat:7.2f}  p={p_val:.4f}  {sig}")

# ── Segmentation summary ──────────────────────────────────────────────────────
total_revenue = df["monetary"].sum()
print(f"\nSegmentation summary:")
print(f"  Customers segmented:  {len(df)}")
print(f"  Clusters (K):         4")
print(f"  Mean silhouette:      {mean_sil:.4f} (good separation)")
print(f"  Total portfolio rev:  ${total_revenue:,}")
print(f"  Champions revenue:    ${df[df['persona']=='Champions']['monetary'].sum():,} ({df[df['persona']=='Champions']['monetary'].sum()/total_revenue*100:.1f}%)")

# ── Inter-cluster distance matrix ─────────────────────────────────────────────
print("\nCluster centroid distance matrix (normalised RFM space):")
_, centers_k4, _ = kmeans(X, k=4, n_init=20)
centroid_dist = np.sqrt(((centers_k4[:, None] - centers_k4[None]) ** 2).sum(axis=2)).round(3)
dist_df = pd.DataFrame(centroid_dist, columns=[f"C{i}" for i in range(4)],
                       index=[f"C{i}" for i in range(4)])
print(dist_df.to_string())