Customer Segmentation with k-means and PCA: Scaling, Choosing k, Stability and Deployment

"""Customer segments for Tallow & Thyme's meal kits. You write the functions marked TODO, one step at a time."""
import numpy as np
import pandas as pd
from sklearn.metrics import adjusted_rand_score, silhouette_score

FEATURES = ["orders_90d", "avg_basket_eur", "veg_share", "weekend_share", "discount_share", "tenure_days",
            "support_tickets"]
LOG_FEATURES = ["avg_basket_eur", "tenure_days"]   # long right tails: a few huge values would dominate


def load(path="customers.csv"):
    return pd.read_csv(path)


# ---------- Step 1: make distances mean something ----------

class Scaler:
    """Learns each column's mean and standard deviation from the data it is fitted on, and applies those."""

    def fit(self, X):
        # TODO (Step 1): store the column means and standard deviations in self.mean_ and self.std_
        # (a std of 0 becomes 1); return self.
        raise NotImplementedError("Step 1: write Scaler.fit()")

    def transform(self, X):
        # TODO (Step 1): (X - self.mean_) / self.std_.
        raise NotImplementedError("Step 1: write Scaler.transform()")


def features(df):
    """The FEATURES columns as a float array, with np.log applied to the LOG_FEATURES."""
    # TODO (Step 1): df[FEATURES] as floats, np.log on the LOG_FEATURES columns, .to_numpy().
    raise NotImplementedError("Step 1: write features()")


def prepare(df):
    """(X, scaler): features(df) standardised by a Scaler fitted on them."""
    # TODO (Step 1): fit a Scaler on features(df) and return (its transform, the scaler).
    raise NotImplementedError("Step 1: write prepare()")


# ---------- Step 2: k-means ----------

def init_centroids(X, k, rng):
    """k-means++: the first centroid is a random row; each next one is a row drawn with probability proportional to
    its squared distance from the nearest centroid chosen so far."""
    raise NotImplementedError("init_centroids() arrives in Step 2")


def assign(X, centroids):
    """Index of the nearest centroid (squared Euclidean distance) for every row."""
    raise NotImplementedError("assign() arrives in Step 2")


def kmeans(X, k, seed=0, max_iter=100):
    """(labels, centroids, inertia). Alternate assign() and moving each centroid to the mean of its rows until no
    label changes; a centroid that loses all its rows keeps its position. inertia = sum of squared distances of rows
    to their centroid."""
    raise NotImplementedError("kmeans() arrives in Step 2")


def best_of(X, k, n_init=5):
    """The kmeans() result with the lowest inertia over seeds 0 .. n_init - 1."""
    raise NotImplementedError("best_of() arrives in Step 2")


# ---------- Step 3: how many segments ----------

def stability(X, k, n_boot=5):
    """Mean adjusted Rand index between best_of(X, k) labels and the labels a model fitted on a bootstrap sample
    (np.random.default_rng(b).choice(len(X), len(X)), b = 0 .. n_boot - 1) gives every row of X."""
    raise NotImplementedError("stability() arrives in Step 3")


def scan_k(X, ks):
    """{k: {"inertia", "silhouette", "stability"}} for each k; silhouette_score(X, labels, sample_size=1500,
    random_state=0) of best_of(X, k)."""
    raise NotImplementedError("scan_k() arrives in Step 3")


def choose_k(scan, min_stability=0.9):
    """The k with the highest silhouette among those whose stability is at least min_stability."""
    raise NotImplementedError("choose_k() arrives in Step 3")


K = None


# ---------- Step 4: a map of the customers ----------

def pca(X, n=2):
    """(scores, components, explained_ratio) from the SVD of the centred X: components are the top n right
    singular vectors (one per row), scores = centred X @ components.T, explained_ratio = each component's share of
    the total variance (all components, so the n values sum to less than 1)."""
    raise NotImplementedError("pca() arrives in Step 4")


def top_loadings(component, names=FEATURES, n=3):
    """The n (name, weight) pairs with the largest absolute weight in one component, largest first."""
    raise NotImplementedError("top_loadings() arrives in Step 4")


# ---------- Step 5: put the segments to work ----------

def profile(df, labels):
    """One row per cluster (index = cluster number): the median of every FEATURES column in its original units,
    plus "share", the fraction of customers in the cluster."""
    raise NotImplementedError("profile() arrives in Step 5")


def segment_new(new_df, scaler, centroids):
    """Cluster numbers for new customers, measured with the scaler and centroids learned from the original ones."""
    raise NotImplementedError("segment_new() arrives in Step 5")