1 Why Models Decay: The World Doesn't Hold Still

Back in Lessons 19–20 you learned to evaluate a model with metrics like accuracy, AUC, and F1 on a held-out test set. That test set was a snapshot — a sample of reality at one moment in time. The implicit assumption behind every supervised model is that the data it sees in production is drawn from the same distribution as its training data. Statisticians call this the i.i.d. assumption: training and future production examples are independent and identically distributed.

That assumption is false almost everywhere, almost immediately. Consider a few concrete examples:

• Credit risk model: trained on applicants from 2022. In 2026, interest rates have changed, a recession has hit certain industries, and the income distribution of applicants has shifted upward as inflation pushed nominal wages higher.
• E-commerce recommender: trained before a new product category (say, refurbished electronics) launched. The model has never seen interaction patterns for this category, so its score distribution for affected users starts behaving strangely.
• Retail demand forecasting: trained on pre-pandemic shopping patterns. Lockdowns push grocery delivery demand up 400% and restaurant-supply demand down 80% — almost overnight, the "normal" the model learned no longer exists.
• Industrial sensor model: a temperature or vibration sensor on a production line gets physically recalibrated or replaced. The raw signal range shifts even though the underlying machine behavior hasn't changed at all.

In every case, the model's weights are bit-for-bit identical to the day it was deployed. Nothing "broke" in the engineering sense. Yet the relationship the model is being asked to apply no longer matches the relationship it learned. This phenomenon is called model decay or model staleness, and it is one of the most common — and most quietly dangerous — failure modes in production ML, because there's no stack trace, no 500 error, no alert. The model happily returns a confident, wrong prediction.

The Monitoring Mental Model

Think of monitoring as comparing two snapshots: a frozen reference distribution (typically your training set, or a known-good recent production window) against a current distribution (a rolling window of recent production traffic — say, the last day or week). If the two snapshots look statistically similar, the model is still operating in familiar territory. If they diverge significantly, something about the world has changed, and the model's learned relationship may no longer apply.

2 Three Kinds of Drift: Precise Definitions

"Drift" is used loosely in casual conversation, but it actually refers to three mathematically distinct phenomena. Getting this right matters because each type calls for a different response. Recall from probability (Lesson 7) that a supervised model is, at its core, an approximation of the conditional distribution P(Y|X) — the probability of the label given the features. We can decompose any change in the world into changes to P(X), P(Y|X), or P(Y).

1. Data Drift / Covariate Shift — P(X) changes, P(Y|X) stays fixed

The distribution of the input features shifts, but the underlying relationship between features and target is unchanged. If you saw the same feature values again, the correct label would still be the same. Example: a bank starts marketing to a younger demographic, so the age distribution of applicants shifts younger — but for any given age and income, the actual probability of default hasn't changed. This is the type of drift PSI and CSI are built to detect, because it shows up directly in the feature/score distributions without needing labels.

2. Concept Drift — P(Y|X) itself changes

The relationship between inputs and the correct output changes. The exact same feature vector now implies a different label than it used to. Example: pre-2020, "books a hotel for 14 consecutive days" was a weak signal of a business traveler; during the pandemic, the same feature pattern became a strong signal of remote-work relocation. The mapping from X to Y itself moved. Concept drift is more dangerous and harder to detect from features alone — it typically requires monitoring prediction performance against (delayed) ground truth, or proxy/business metrics, since the feature distribution can look perfectly normal while the correct answer has quietly flipped.

3. Label / Prior Drift — P(Y) changes

The overall base rate of the outcome shifts, independent of any single feature. Example: a fraud model trained when fraud occurred in 0.3% of transactions now operates in an environment where a new fraud ring has pushed the base rate to 1.1%. Feature distributions might look identical for individual legitimate vs fraudulent transactions, but the mix has changed. This affects calibration (a 5% predicted-fraud score now corresponds to a different real-world probability) and can also be partly visible in the prediction-score distribution, which is exactly what PSI on model outputs is designed to catch.

3 The Population Stability Index (PSI)

PSI originated in the credit-risk industry decades ago as a way to monitor whether a credit scorecard's output distribution had shifted from the population it was built on. It has since become the de facto standard metric for monitoring any model's prediction distribution (and, as you'll see in Section 4, any individual feature's distribution) in production ML.

The Core Idea: Bucket, Compare, Weight by Divergence

PSI compares two distributions — a fixed reference distribution (almost always your training data, or a "golden" early-production window) and a current distribution (a recent slice of live traffic) — by binning both into the same buckets and measuring how much the proportion of observations in each bucket has shifted.

Step 1: Bin the Reference Distribution

You first decide on bin edges using only the reference distribution — this is critical, the bins are frozen once defined and reused for every future comparison. The standard approach (recall histograms/binning from Lesson 6) is to use deciles of the reference distribution: 10 bins, each containing exactly 10% of the reference population, with edges at the 10th, 20th, ..., 90th percentiles. Equal-frequency (quantile) binning is preferred over equal-width binning because it guarantees every bin starts with a non-trivial reference count, which avoids division-by-zero and unstable ratios later in the formula.

Step 2: Compute Bucket Percentages for Both Populations

Using the exact same bin edges, compute what percentage of the reference population falls in each bin, and separately what percentage of the current population falls in each bin.

Step 3: Apply the PSI Formula

For each bin i, with reference_pct_i and current_pct_i as the proportions (not counts) of each population in that bin:

# PSI formula, applied bin by bin and summed:
#
#   PSI = Σ_i  (current_pct_i − reference_pct_i) × ln(current_pct_i / reference_pct_i)
#
# Intuition for a single bin:
#   - If current_pct_i == reference_pct_i  → that bin contributes 0 (no shift)
#   - If current_pct_i > reference_pct_i   → positive contribution (bin over-represented now)
#   - If current_pct_i < reference_pct_i   → STILL a positive contribution, because both the
#     difference term and the log term flip sign together — every bin contributes >= 0
#   - Bins with a LARGER relative change contribute disproportionately more (the log term
#     amplifies large ratios), so PSI is most sensitive to bins that have nearly emptied out
#     or suddenly filled up, not just bins with small absolute movement.
print("PSI is a sum of always-non-negative terms: 0 means identical distributions")

Each term in the sum is non-negative (it's a form of symmetrized KL divergence applied bin-by-bin), so PSI itself is always ≥ 0, with 0 meaning the two distributions are identical across every bin.

From-Scratch PSI Implementation

Here is a complete, dependency-light implementation using only NumPy and Pandas. Notice that the bin edges are computed once from the reference distribution and then reused — this is the detail most homemade implementations get wrong.

import numpy as np
import pandas as pd

def compute_psi(reference: np.ndarray, current: np.ndarray, n_bins: int = 10,
                 epsilon: float = 1e-4) -> dict:
    """
    Compute the Population Stability Index between a reference distribution
    (e.g. training data or model scores at deploy time) and a current
    distribution (e.g. last 7 days of production traffic).

    Parameters
    ----------
    reference : 1D array of numeric values (the BASELINE — bins are derived from this)
    current   : 1D array of numeric values (the CURRENT window being checked)
    n_bins    : number of quantile bins to use (10 deciles is the industry default)
    epsilon   : small constant added to avoid log(0) / division-by-zero when a
                bin's percentage is exactly zero in one of the two populations

    Returns
    -------
    dict with the overall psi value and a per-bin breakdown DataFrame.
    """
    reference = np.asarray(reference, dtype=float)
    current   = np.asarray(current, dtype=float)

    # Step 1: derive bin edges from the REFERENCE distribution only (deciles).
    # np.unique handles the edge case where many reference values are identical
    # (e.g. a feature with a spike at 0), which would otherwise create empty bins.
    quantiles  = np.linspace(0, 1, n_bins + 1)
    bin_edges  = np.unique(np.quantile(reference, quantiles))
    bin_edges[0]  = -np.inf   # extend outer edges so new data is never "out of range"
    bin_edges[-1] =  np.inf

    # Step 2: assign every observation in both populations to a bin
    ref_bins = pd.cut(reference, bins=bin_edges, include_lowest=True)
    cur_bins = pd.cut(current,   bins=bin_edges, include_lowest=True)

    ref_counts = ref_bins.value_counts(sort=False)
    cur_counts = cur_bins.value_counts(sort=False)

    ref_pct = (ref_counts / ref_counts.sum()).clip(lower=epsilon)
    cur_pct = (cur_counts / cur_counts.sum()).clip(lower=epsilon)

    # Step 3: PSI formula, bin by bin
    psi_per_bin = (cur_pct - ref_pct) * np.log(cur_pct / ref_pct)

    breakdown = pd.DataFrame({
        "bin":            ref_pct.index.astype(str),
        "reference_pct":  ref_pct.values,
        "current_pct":    cur_pct.values,
        "psi_contribution": psi_per_bin.values,
    })

    return {
        "psi": float(psi_per_bin.sum()),
        "breakdown": breakdown,
    }


# ── Example: a credit score feature drifting upward in production ──
np.random.seed(42)
reference_scores = np.random.normal(loc=650, scale=50, size=10_000)   # training-time scores
current_scores   = np.random.normal(loc=680, scale=55, size=2_000)    # last week's scores (shifted up)

result = compute_psi(reference_scores, current_scores, n_bins=10)
print(f"Overall PSI: {result['psi']:.4f}\n")
print(result["breakdown"].to_string(index=False))

Interpreting PSI: The Threshold Table

PSI values are interpreted using thresholds that originated in retail credit scoring and have become the near-universal industry convention across ML monitoring tools:

In the example above, PSI ≈ 0.119 lands in the "moderate shift, investigate" band — consistent with the fact that we deliberately shifted the mean from 650 to 680 (about 0.5 standard deviations), a meaningful but not extreme change.

4 The Characteristic Stability Index (CSI): Finding the Culprit Feature

CSI uses exactly the same mathematical machinery as PSI — same binning strategy, same formula — but applied independently to each input feature rather than to the final model score. Some practitioners simply say "we compute PSI per-feature"; "CSI" is the traditional credit-risk name for that specific application. If PSI on the model's output score is the smoke alarm for the whole house, CSI is walking room to room with a thermal camera to find which room is actually on fire.

Why You Need Both

A model can have a perfectly stable output-score PSI even while one input feature has drifted substantially, if other features happen to compensate, or if the drifted feature has low importance in the model. Conversely, a single important feature drifting can move the output-score PSI a lot. Computing CSI for every feature gives you a ranked list of "which inputs look the most different from training," which is invaluable for root-causing an alert — and for distinguishing genuine data drift from a pipeline bug (e.g., a unit conversion bug that suddenly reports income in cents instead of dollars will show up as an enormous CSI on exactly one feature).

From-Scratch CSI Implementation

Because CSI is "PSI applied per column," the cleanest implementation reuses the same compute_psi function from Section 3 and simply loops over the feature columns of a reference and current DataFrame:

import numpy as np
import pandas as pd

def compute_csi_report(reference_df: pd.DataFrame, current_df: pd.DataFrame,
                        feature_cols: list, n_bins: int = 10) -> pd.DataFrame:
    """
    Compute the Characteristic Stability Index for every feature in feature_cols,
    by applying the PSI formula independently to each column.

    Returns a DataFrame sorted by descending CSI (worst-drifting feature first),
    with a verdict column applying the standard PSI thresholds.
    """
    rows = []
    for col in feature_cols:
        ref_values = reference_df[col].dropna().values
        cur_values = current_df[col].dropna().values

        result = compute_psi(ref_values, cur_values, n_bins=n_bins)
        csi = result["psi"]

        if csi < 0.10:
            verdict = "stable"
        elif csi < 0.25:
            verdict = "moderate shift — investigate"
        else:
            verdict = "significant shift — retrain/recalibrate"

        rows.append({"feature": col, "csi": round(csi, 4), "verdict": verdict})

    report = pd.DataFrame(rows).sort_values("csi", ascending=False).reset_index(drop=True)
    return report


# ── Example: a credit risk model with 5 features, one of which has drifted ──
np.random.seed(7)
n_ref, n_cur = 8_000, 1_500

reference_df = pd.DataFrame({
    "annual_income":      np.random.lognormal(mean=10.8, sigma=0.4, size=n_ref),
    "credit_utilization": np.random.beta(2, 5, size=n_ref),
    "age":                np.random.normal(40, 12, size=n_ref).clip(18, 85),
    "num_open_accounts":  np.random.poisson(4, size=n_ref),
    "months_since_default": np.random.exponential(36, size=n_ref),
})

current_df = pd.DataFrame({
    # income shifted up ~18% (inflation + applicant mix change) — the "culprit"
    "annual_income":      np.random.lognormal(mean=10.97, sigma=0.42, size=n_cur),
    "credit_utilization": np.random.beta(2, 5, size=n_cur),          # unchanged
    "age":                np.random.normal(41, 12, size=n_cur).clip(18, 85),  # tiny, harmless shift
    "num_open_accounts":  np.random.poisson(4, size=n_cur),          # unchanged
    "months_since_default": np.random.exponential(36, size=n_cur),   # unchanged
})

feature_cols = list(reference_df.columns)
csi_report = compute_csi_report(reference_df, current_df, feature_cols)
print(csi_report.to_string(index=False))

This is exactly the workflow that catches drift before it shows up as a measurable accuracy drop: the output-score PSI alert tells you something moved; the CSI report tells you it was annual_income specifically, with every other feature stable. That immediately focuses the investigation — is this a real macroeconomic shift, an upstream data pipeline bug, or a change in applicant mix from a new marketing channel?

5 Choosing Bins Well: Practical Pitfalls

The binning strategy is where most homegrown PSI implementations go subtly wrong. A few practical rules, learned the hard way by teams running this in production:

Freeze the Bins at Reference Time

Bin edges must be computed once, from the reference distribution, and then reused unchanged for every future comparison window. If you recompute quantile edges from the current window each time, you are comparing the current distribution's shape against a constantly moving target rather than against the original baseline — this silently disables drift detection, since equal-frequency bins drawn from the current data will always show roughly 10% per bin by construction.

Decile Binning Is a Default, Not a Law

10 equal-frequency bins is the standard default and works well for most continuous features with reasonably smooth distributions. For features with a small number of distinct values (e.g. an integer count feature like num_open_accounts ranging 0–8), fewer bins (or exact-value bins) often make more sense — 10 quantile bins on a feature with only 9 distinct values will produce duplicate edges and uneven bins. For very large reference samples and smooth continuous features, 20 bins gives finer resolution at the cost of more noise in low-traffic-window comparisons.

Small Current-Window Sample Sizes Inflate PSI

PSI is a sample statistic, and like any statistic computed on percentages, it gets noisier as the current window's sample size shrinks. Monitoring against a current window of only 50 predictions can produce a misleadingly high PSI purely from sampling noise, not real drift. A common rule of thumb is to require at least a few hundred observations (ideally 1,000+) in the current window before trusting a PSI reading, and to widen the time window (e.g. from 1 day to 7 days) for low-traffic models or rare feature values.

The Epsilon Floor Matters More Than It Looks

Whenever a bin's percentage is exactly 0 in either population, ln(0) or division by 0 will crash or return inf/nan. Clipping percentages to a small floor (commonly 0.0001) before taking the log avoids this, but be aware that this floor itself becomes the dominant driver of PSI for the affected bin if a category has fully appeared or disappeared — which is usually exactly the signal you want to surface, not suppress.

6 Building an Automated Drift Monitoring Pipeline

PSI and CSI are only useful if they run continuously and alert someone. Recall from Lessons 65–67 how models are served behind an endpoint that logs every request and prediction; drift monitoring sits downstream of that logging, as a scheduled batch job.

Pipeline Architecture

A typical production drift-monitoring pipeline has four stages, usually run as a daily or hourly scheduled job (cron, Airflow DAG, or a managed scheduler like SageMaker Model Monitor or a cloud function on a timer):

1. Baseline capture (one-time, at deploy time): save the reference distribution's raw values (or pre-computed bin edges + reference percentages) for the model score and every input feature, as an artifact alongside the model.
2. Window extraction (scheduled): pull the last N hours/days of logged predictions and their input features from the serving logs or a feature store.
3. Compute PSI/CSI: run the score-level PSI and the per-feature CSI report against the frozen reference baseline.
4. Alert and record: compare against thresholds, write the results to a monitoring dashboard/time series, and fire an alert (Slack, PagerDuty, email) if any threshold is breached.

import json
import datetime as dt
import numpy as np
import pandas as pd
from pathlib import Path

class DriftMonitor:
    """
    A minimal scheduled drift-monitoring job. In production this would be
    triggered by a scheduler (Airflow, cron, a cloud Lambda on a timer) and
    would read from a real feature store / prediction log instead of an
    in-memory DataFrame.
    """

    def __init__(self, reference_path: str, score_col: str, feature_cols: list,
                 psi_warn: float = 0.10, psi_alert: float = 0.25, n_bins: int = 10):
        self.reference = pd.read_parquet(reference_path)   # frozen baseline snapshot
        self.score_col = score_col
        self.feature_cols = feature_cols
        self.psi_warn = psi_warn
        self.psi_alert = psi_alert
        self.n_bins = n_bins

    def run(self, current_window: pd.DataFrame) -> dict:
        """Run one monitoring cycle against a window of recent production data."""
        timestamp = dt.datetime.utcnow().isoformat()

        # 1. Score-level PSI — the cheap, first-line tripwire
        score_psi = compute_psi(
            self.reference[self.score_col].values,
            current_window[self.score_col].values,
            n_bins=self.n_bins,
        )["psi"]

        # 2. Per-feature CSI — the diagnostic breakdown
        csi_report = compute_csi_report(
            self.reference, current_window, self.feature_cols, n_bins=self.n_bins
        )

        # 3. Decide overall status
        breached_features = csi_report[csi_report["csi"] >= self.psi_alert]
        if score_psi >= self.psi_alert or not breached_features.empty:
            status = "ALERT"
        elif score_psi >= self.psi_warn or (csi_report["csi"] >= self.psi_warn).any():
            status = "WARN"
        else:
            status = "OK"

        result = {
            "timestamp": timestamp,
            "status": status,
            "score_psi": round(float(score_psi), 4),
            "n_current_rows": len(current_window),
            "top_drifting_features": csi_report.head(3).to_dict(orient="records"),
        }

        self._notify_if_needed(result)
        return result

    def _notify_if_needed(self, result: dict):
        if result["status"] in ("WARN", "ALERT"):
            # In production: post to Slack/PagerDuty/email here.
            print(f"[{result['status']}] Drift monitor — score PSI = {result['score_psi']}")
            for feat in result["top_drifting_features"]:
                print(f"    {feat['feature']:<20s} csi={feat['csi']:.4f}  ({feat['verdict']})")
        else:
            print(f"[OK] Drift monitor — score PSI = {result['score_psi']} — no action needed")


# Example of a single scheduled run (illustrative — paths/data are placeholders)
# monitor = DriftMonitor(
#     reference_path="s3://ml-artifacts/credit-model/v3/reference_baseline.parquet",
#     score_col="risk_score",
#     feature_cols=["annual_income", "credit_utilization", "age",
#                   "num_open_accounts", "months_since_default"],
# )
# today_window = load_predictions_from_feature_store(days=1)
# report = monitor.run(today_window)
print("DriftMonitor.run() is designed to be invoked by a daily scheduled job")

7 Beyond PSI: Other Drift Detection Tools (Context)

PSI and CSI are popular because they are simple, interpretable as a single number, and have decades of industry-accepted thresholds — but they are not the only tools in this space. It's worth knowing the alternatives so you can recognize them in other teams' tooling or pick the right tool for a specific situation.

KL Divergence

Kullback-Leibler divergence measures how one probability distribution diverges from a reference distribution — PSI is in fact mathematically a symmetrized variant of KL divergence applied to binned data (PSI is sometimes written as KL(current‖reference) + KL(reference‖current) computed bin-wise). Raw KL divergence is asymmetric (KL(P‖Q) ≠ KL(Q‖P)) and has no universally agreed interpretation thresholds, which is part of why the symmetric, threshold-calibrated PSI became the practical industry standard instead.

Kolmogorov-Smirnov (K-S) Test

The K-S test compares two continuous distributions using their empirical cumulative distribution functions (CDFs), producing both a distance statistic and a p-value. Unlike PSI, it doesn't require choosing bins at all, and it gives you a formal hypothesis-test framework ("reject the null hypothesis that these come from the same distribution at p < 0.05"). It is sensitive primarily to differences in the overall shape/location of two distributions and works well for continuous numeric features, but doesn't extend as naturally to categorical features the way a bucketed PSI/CSI does, and a statistically significant K-S result on a huge production sample doesn't always indicate a practically meaningful shift.

Tooling: Evidently AI and whylogs

Rather than hand-rolling every check, most production teams reach for a dedicated data/ML monitoring library. Evidently AI is an open-source Python library that generates drift reports and dashboards, computing PSI, K-S tests, Wasserstein distance, and more, per column, with built-in HTML reports and threshold configuration. whylogs (WhyLabs) takes a "logging" approach — it computes lightweight statistical profiles of data at the point of ingestion (so you don't need to store raw data for monitoring) and compares profiles over time, integrating with most ML pipeline orchestrators. Both are worth knowing by name; in a real production setting you would typically configure one of these rather than maintaining a bespoke PSI script — but understanding the PSI/CSI math underneath, as you now do, is what lets you correctly interpret, debug, and tune whatever the dashboard shows you.