Classifier Two-Sample Tests (C2STs)

A classifier two-sample test turns the question are these two samples from the same distribution into a supervised problem. Label the baseline rows 0 and the current rows 1, train a classifier to tell them apart, and read its accuracy. If both samples come from the same distribution and the classes are balanced, the classifier has no population-level signal to find; held-out accuracy still varies in a finite sample.

That is the whole idea, and it scales to images, text and embeddings where KS shift has no multivariate form. Two details decide whether it is a test or a self-portrait.

Which accuracy

Lopez-Paz and Oquab define the statistic on a held-out split, not on the data the classifier was fitted to. The reason shows up immediately when there is nothing to find.

import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split

rng = np.random.default_rng(1234)
n, d = 500, 10
A = rng.normal(size=(n, d))
B = rng.normal(size=(n, d))
X = np.vstack([A, B])
y = np.concatenate([np.zeros(n), np.ones(n)])
Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=0, stratify=y)

print("%-30s %14s %14s" % ("classifier", "training acc", "held-out acc"))
for name, clf in (("logistic regression", LogisticRegression(max_iter=5000)),
                  ("random forest, 300 trees", RandomForestClassifier(
                      n_estimators=300, random_state=0)),
                  ("random forest, no pruning", RandomForestClassifier(
                      n_estimators=300, min_samples_leaf=1, max_features=None,
                      random_state=0))):
    clf.fit(Xtr, ytr)
    print("%-30s %14.4f %14.4f" % (name, clf.score(Xtr, ytr), clf.score(Xte, yte)))

classifier                       training acc   held-out acc
logistic regression                    0.5580         0.5120
random forest, 300 trees               1.0000         0.4700
random forest, no pruning              1.0000         0.4620

Both samples are ten-dimensional standard normals. There is no difference. The forests reach 1.0000 on the rows they were fitted to and 0.4700 and 0.4620 on rows they have not seen; even logistic regression finds 0.5580 in sample. A C2ST computed in sample does not overstate the drift a little. A sufficiently flexible classifier can memorize the training labels and produce high training accuracy even when both samples come from the same distribution. That accuracy is not valid evidence of drift.

The null is already known

For iid fair random test labels independent of features and a fixed trained classifier, the correct count has a Binomial\((n_{te},1/2)\) null. The code instead fixes class counts by stratification: those correctness indicators need not be independent Bernoulli(1/2). An always-zero classifier then has exactly 50% accuracy on a balanced test set and zero variance. Use a null matching the design.

from scipy.stats import hypergeom


def conditional_p(y, pred):
    # Exact fixed-count label-permutation tail, conditional on held-out predictions.
    successes = int(np.sum((y == 1) & (pred == 1)))
    return hypergeom.sf(successes - 1, len(y), int(np.sum(y)), int(np.sum(pred)))


rng = np.random.default_rng(88)
d, reps = 8, 600
print("%7s %7s %11s %11s %11s" % ("n / side", "n_te", "mean acc", "sd of acc", "reject 5%"))
for n in (100, 400, 1600):
    accs, pvals = [], []
    for _ in range(reps):
        X = np.vstack([rng.normal(size=(n, d)), rng.normal(size=(n, d))])
        y = np.concatenate([np.zeros(n), np.ones(n)])
        Xtr, Xte, ytr, yte = train_test_split(X, y, test_size=0.5, random_state=0,
                                              stratify=y)
        pred = LogisticRegression(max_iter=5000).fit(Xtr, ytr).predict(Xte)
        accs.append(np.mean(pred == yte))
        pvals.append(conditional_p(yte, pred))
    accs = np.array(accs)
    n_te = len(yte)
    print("%7d %7d %11.4f %11.4f %10.1f%%"
          % (n, n_te, accs.mean(), accs.std(ddof=1),
             (np.array(pvals) <= .05).mean() * 100))

n / side    n_te    mean acc   sd of acc   reject 5%
    100     100      0.5033      0.0490        3.8%
    400     400      0.5001      0.0247        4.5%
   1600    1600      0.4997      0.0125        4.8%

Across 600 replications per size, the mean held-out accuracy stays at 0.50 and its standard deviation shrinks from 0.0490 to 0.0125 as the test set grows from 100 to 1,600 rows, roughly \( 1/(2\sqrt{n_{te}}) \). The exact test rejects at 3.8%, 4.5% and 4.8%, at or below the nominal 5% each time.

For the stratified examples, condition on the trained classifier and its held-out predictions, then permute test labels while preserving class counts. The number of positive labels among predicted-positive rows has a hypergeometric distribution, so its upper-tail probability gives the exact conditional permutation p-value without simulation. In conditional_p, successes counts the test rows that are labeled 1 and predicted 1, and the test-set size, the number of 1 labels and the number of positive predictions fix the hypergeometric distribution. The test set must remain untouched by training, tuning or classifier selection, and test rows must be exchangeable under the null.

Chance level is a statement about the classifier

The test inherits whatever the classifier can express. Two samples with identical means and identical spreads that differ only in how the features move together are a clean case, because no straight line separates them. The first block tests the same correlated pair three ways: a logistic regression on the two raw features, the same model given the product of the two features as a third column, and a random forest on the raw features. Each test uses a stratified 50/50 split and the conditional_p defined above, and this block and the next continue the same session.

def correlated_pair(rng, n, rho):
    base = rng.normal(size=(n, 2))
    a, b = rng.normal(size=n), rng.normal(size=n)
    cur = np.column_stack([a, rho * a + np.sqrt(1 - rho ** 2) * b])
    return (np.vstack([base, cur]),
            np.concatenate([np.zeros(n), np.ones(n)]))


rng = np.random.default_rng(101)
X, y = correlated_pair(rng, 2000, 0.8)
prod = np.column_stack([X[:, 0], X[:, 1], X[:, 0] * X[:, 1]])

print("%-42s %10s %11s" % ("test", "accuracy", "p"))
for name, feats, clf in (
        ("logistic regression, two raw features", X, LogisticRegression(max_iter=5000)),
        ("logistic regression, plus the product", prod, LogisticRegression(max_iter=5000)),
        ("random forest, two raw features", X, RandomForestClassifier(
            n_estimators=300, random_state=0))):
    Ftr, Fte, ytr, yte = train_test_split(feats, y, test_size=0.5, random_state=0,
                                          stratify=y)
    pred = clf.fit(Ftr, ytr).predict(Fte)
    acc = np.mean(pred == yte)
    print("%-42s %10.4f %11.2e" % (name, acc, conditional_p(yte, pred)))

test                                         accuracy           p
logistic regression, two raw features          0.4950    6.89e-01
logistic regression, plus the product          0.6065    5.22e-23
random forest, two raw features                0.6555    3.07e-45

On the same held-out features, a linear classifier has little power to detect this correlation change (accuracy 0.4950, p = 0.689); adding the product feature (0.6065) or using a forest (0.6555) exposes it, with p-values near 5e-23 and 3e-45. Failure to reject says this test found insufficient evidence, not that the distributions match.

The second experiment asks a different question: with a fixed total number of rows, how should they be divided between fitting the classifier and testing it? It draws 150 rows per side from a weaker alternative (correlation 0.45), gives the classifier the product feature, and varies the share of the 300 rows held out for testing. Each setting is repeated 400 times; the last column is the fraction of repetitions in which conditional_p is at most 0.05.

rng2 = np.random.default_rng(55)
print("%-14s %8s %10s %12s" % ("train / test", "n_te", "mean acc", "power at 5%"))
for frac in (0.95, 0.8, 0.5, 0.2, 0.05):
    accs, pvals = [], []
    for _ in range(400):
        Xw, yw = correlated_pair(rng2, 150, 0.45)
        Fw = np.column_stack([Xw[:, 0], Xw[:, 1], Xw[:, 0] * Xw[:, 1]])
        Ftr, Fte, ytr, yte = train_test_split(Fw, yw, test_size=frac, random_state=0,
                                              stratify=yw)
        pred = LogisticRegression(max_iter=5000).fit(Ftr, ytr).predict(Fte)
        accs.append(np.mean(pred == yte))
        pvals.append(conditional_p(yte, pred))
    accs, n_te = np.array(accs), len(yte)
    print("%-14s %8d %10.4f %11.4f"
          % ("%d / %d" % (100 - frac * 100, frac * 100), n_te, accs.mean(),
             (np.array(pvals) <= .05).mean()))

train / test       n_te   mean acc  power at 5%
5 / 95              285     0.5276      0.3450
20 / 80             240     0.5465      0.4175
50 / 50             150     0.5588      0.3775
80 / 20              60     0.5600      0.1850
95 / 5               15     0.5752      0.0800

Held-out accuracy rises from 0.5276 to 0.5752 as the classifier gets more training rows, yet power peaks at 0.4175 for a 20/80 split and falls to 0.0800 when only 15 test rows remain. The split experiment reports held-out accuracy and rejection frequency under the stated shift. A larger training set may improve the classifier while leaving fewer test observations. The best split depends on classifier, alternative and sample size; no split universally maximizes power.

Choose the classifier and split before examining test results, and calibrate significance to the actual sampling design. MMD and energy distance offer complementary tests, with their own finite-sample power limitations. All of them measure distribution differences, not business harm.

References

• Lopez-Paz, D., & Oquab, M. (2017). Revisiting classifier two-sample tests. International Conference on Learning Representations (ICLR).
• Gretton, A., Borgwardt, K. M., Rasch, M. J., Schölkopf, B., & Smola, A. (2012). A kernel two-sample test. Journal of Machine Learning Research, 13, 723–773.