cd/** permutation p-value omits the observed statistic at 10 sites

[shivamlalakiya]

Summary

Ten sites across the cd/ drift detectors compute a permutation p-value as
k / n_permutations, where k is the count of permuted statistics at least
as extreme as the observed one, but the observed statistic itself is
computed separately (with permute=False) and is not one of the k
permuted draws counted:

cd/pytorch/mmd.py:142               p_val = (mmd2 <= mmd2_permuted).float().mean()
cd/tensorflow/mmd.py:127            p_val = (mmd2 <= mmd2_permuted).mean()
cd/keops/mmd.py:181                 p_val = (mmd2 <= mmd2_permuted).float().mean()
cd/pytorch/lsdd.py:175              p_val = (lsdd <= lsdd_permuted).float().mean()
cd/tensorflow/lsdd.py:158           p_val = tf.reduce_mean(tf.cast(lsdd <= lsdd_permuted, float))
cd/pytorch/learned_kernel.py:220    p_val = (mmd2 <= mmd2_permuted).float().mean()
cd/tensorflow/learned_kernel.py:197 p_val = (mmd2 <= mmd2_permuted).mean()
cd/keops/learned_kernel.py:240      p_val = (mmd2 <= mmd2_permuted).float().mean()
cd/pytorch/context_aware.py:174     p_val = (stat <= permuted_stats).float().mean()
cd/tensorflow/context_aware.py:163  p_val = tf.reduce_mean(tf.cast(stat <= permuted_stats, float))

All ten confirmed against the current alibi-detect 0.13.0 wheel from PyPI
(same line numbers).

The usual construction for a permutation p-value (Phipson & Smyth,
Permutation P-values Should Never Be Zero, SAGMB 9(1) Art. 39, 2010,
doi:10.2202/1544-6115.1585) includes the observed statistic among the B+1
values compared, (1 + k) / (B + 1), precisely so the p-value can't land on
zero. The shipped form is k / B.

Reproduction

cd/mmd.py:35 defaults n_permutations=100. 400 pairs of i.i.d. N(0,1)
draws (no actual drift), through the public MMDDrift API, pytorch
backend:

>>> import numpy as np
>>> from alibi_detect.cd import MMDDrift
>>> ps = []
>>> for _ in range(400):
...     x_ref = np.random.randn(50, 4).astype(np.float32)
...     x_test = np.random.randn(50, 4).astype(np.float32)
...     d = MMDDrift(x_ref, backend='pytorch', p_val=.05, n_permutations=100)
...     ps.append(float(d.predict(x_test)['data']['p_val']))
>>> ps = np.array(ps)
>>> ps.min()
0.0
>>> (ps == 0).mean(), (ps == 0).sum()
(0.015, 6)
>>> (ps <= 0.05).mean()
0.065
>>> sorted(set(np.round(ps[ps < 0.05], 6).tolist()))
[0.0, 0.01, 0.02, 0.03, 0.04]

A p-value of exactly zero on exchangeable, no-drift data, at a rate (6/400)
in the same range as the corrected floor 1/(n_permutations+1) = 1/101 ≈ 0.0099.

Second, cd/base.py:642 (and the same pattern at :282, :483, :804, :1202,
:1002/:1005) decides with drift_pred = int(p_val < self.p_val), a strict
inequality on the k/B lattice. Combined, the realised rejection rate can
exceed the nominal one whenever alpha * n_permutations is not an integer.
Worked directly: at B=101, alpha=0.05, the strict decision rejects for
k in {0,...,5}, i.e. 6 of the 102 equally spaced lattice points, giving
a realised level of 6/102 = 0.05882 against the nominal 0.05:

>>> import math
>>> B, alpha = 101, 0.05
>>> n_reject_k = math.floor(alpha * B - 1e-12) + 1
>>> n_reject_k, n_reject_k / (B + 1)
(6, 0.058823529411764705)

Third, the same missing addend appears in a standalone public helper with
no callers inside the package, so no existing test exercises it:
utils/statstest.py:37-38:

k += dist <= dist_permutation
...
return k / n_permutations, dist, dist_permutations

Scope

This isn't a one-line fix. It's the same reduction, spelled slightly
differently by each backend/detector-family pair, at all ten sites listed
above (three backends: pytorch, tensorflow, keops; four detector families:
MMD, LSDD, learned-kernel, context-aware MMD), plus the standalone helper.
A general shape would be (1 + k) / (n_permutations + 1) per site, with
the observed statistic folded into the comparison instead of being
computed and excluded separately. Filing as a single issue rather than ten
separate ones since it's one construction repeated, but happy to split
into one PR per backend if that's easier to review.

One related, much smaller item while in this code: utils/statstest.py:63
does below_threshold = p_sorted < q_threshold where the standard BH
procedure rejects at <=. That one is boundary-only and conservative in
the other direction (e.g. fdr([0.05], 0.05) returns no rejections where
the textbook procedure rejects), so it's a much smaller and different
question from the ten sites above; mentioning it here rather than opening
a separate report since it's the same file.