Statistical Detection

PSI, KL divergence, KS test, and Jensen-Shannon divergence: the four statistical methods that catch drift before users do, with worked Python examples and decision guidance.

The Core Problem These Methods Solve

You have a baseline: the distribution of inputs (or outputs) at deployment time. You have a current window: the distribution of the same features over the last day, week, or batch. The question is: has the distribution changed significantly, or is the difference just normal sampling variation?

Statistical drift detection methods answer that question with a number. Each method produces a score; you compare that score to a threshold and make a binary decision: drift detected or not detected. The methods differ in what they measure, what they assume, and how they handle edge cases.

Method 1: Population Stability Index (PSI)

PSI is the workhorse of drift detection in credit risk, insurance, and enterprise ML. It compares two distributions by bucketing values and measuring how much each bucket's proportion changed.

How it works

Divide the feature's value range into N buckets (typically 10). For each bucket, calculate the expected proportion (from the baseline) and the actual proportion (from the current window). PSI is the sum over all buckets of:

Worked example

A customer age feature at training had this distribution (simplified to 5 buckets):

PSI = 0.0811 + 0.0077 + 0.0269 + 0.0409 + 0.0067 = 0.163

Interpreting PSI

Our example PSI of 0.163 falls in the "investigate" range: a 13-point swing in the youngest cohort is significant and warrants investigation.

import numpy as np

def calculate_psi(baseline, current, buckets=10, eps=1e-4):
    """Population Stability Index between two arrays."""
    breakpoints = np.percentile(baseline, np.linspace(0, 100, buckets + 1))
    breakpoints[0] = -np.inf
    breakpoints[-1] = np.inf

    base_counts = np.histogram(baseline, bins=breakpoints)[0]
    curr_counts = np.histogram(current, bins=breakpoints)[0]

    base_pct = base_counts / len(baseline) + eps
    curr_pct = curr_counts / len(current) + eps

    psi = np.sum((curr_pct - base_pct) * np.log(curr_pct / base_pct))
    return psi

# Usage
baseline_ages = ...   # array from training data
current_ages  = ...   # array from last 7 days of production traffic
psi = calculate_psi(baseline_ages, current_ages)
print(f"PSI: {psi:.4f} - {'Drift detected' if psi > 0.2 else 'Monitor' if psi > 0.1 else 'Stable'}")

Method 2: Kolmogorov-Smirnov (KS) Test

The KS test is a non-parametric test that measures the maximum absolute difference between the cumulative distribution functions (CDFs) of two samples. Unlike PSI, it makes no assumption about binning and uses the data's actual distribution.

The KS statistic D is the largest vertical gap between the two CDFs. When D is large (and the p-value is small), the distributions are significantly different.

from scipy import stats

# Two-sample KS test
ks_stat, p_value = stats.ks_2samp(baseline_feature, current_feature)

print(f"KS Statistic: {ks_stat:.4f}")
print(f"p-value: {p_value:.4f}")

if p_value < 0.05:
    print("Drift detected: distributions are significantly different")
else:
    print("No significant drift detected at alpha=0.05")

Method 3: KL Divergence

Kullback-Leibler (KL) divergence measures the information loss when one distribution is used to approximate another. Unlike PSI and KS, it is asymmetric: KL(P||Q) is not equal to KL(Q||P). This makes it useful when you have a clear "true" distribution (the baseline) and want to measure how far the current distribution departs from it.

KL divergence has a critical limitation: if any bucket in the current distribution has zero probability while the baseline has non-zero probability, the calculation is undefined. Always add a small epsilon (smoothing constant) or use Jensen-Shannon instead for production monitoring.

Method 4: Jensen-Shannon (JS) Divergence

JS divergence is the symmetric, bounded version of KL divergence. It is bounded between 0 (identical distributions) and 1 (completely disjoint distributions, using base-2 logarithm). Because it is symmetric and handles zero probabilities gracefully, it is the preferred method for monitoring output distributions and embedding space similarity.

from scipy.spatial.distance import jensenshannon
import numpy as np

def distribution_to_histogram(data, bins=50):
    """Convert raw data to a normalized probability distribution."""
    counts, _ = np.histogram(data, bins=bins, density=True)
    return counts / counts.sum()

baseline_dist = distribution_to_histogram(baseline_outputs)
current_dist  = distribution_to_histogram(current_outputs)

js_dist = jensenshannon(baseline_dist, current_dist)  # range: 0 to 1
print(f"JS Distance: {js_dist:.4f}")
# Threshold guidance: > 0.1 is worth investigating; > 0.2 is significant drift

Authoritative References

The implementations above use widely available open-source libraries. The following official documentation pages are stable references for the underlying functions:

• SciPy: scipy.stats.ks_2samp - Two-sample Kolmogorov-Smirnov test
• SciPy: scipy.spatial.distance.jensenshannon - Jensen-Shannon distance
• NumPy: numpy.histogram - Histogram computation used in PSI and JS

Choosing the Right Method for Your Data