Metadata-Version: 2.5
Name: pted
Version: 2.0.0
Summary: Implementation of a Permutation Test using the Energy Distance for two sample tests and posterior coverage tests
Project-URL: Homepage, https://github.com/ConnorStoneAstro/pted
Project-URL: Documentation, https://github.com/ConnorStoneAstro/pted
Project-URL: Repository, https://github.com/ConnorStoneAstro/pted
Project-URL: Issues, https://github.com/ConnorStoneAstro/pted/issues
Author-email: Connor Stone <connorstone628@gmail.com>
License: MIT License
        
        Copyright (c) 2025 Connor Stone, PhD
        
        Permission is hereby granted, free of charge, to any person obtaining a copy
        of this software and associated documentation files (the "Software"), to deal
        in the Software without restriction, including without limitation the rights
        to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the Software is
        furnished to do so, subject to the following conditions:
        
        The above copyright notice and this permission notice shall be included in all
        copies or substantial portions of the Software.
        
        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
        LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
        OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
        SOFTWARE.
License-File: LICENSE
Keywords: Bayesian,machine learning,pytorch,statistics
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Requires-Python: >=3.11
Requires-Dist: numpy
Requires-Dist: scipy
Requires-Dist: tqdm
Provides-Extra: dev
Requires-Dist: jax<1,>=0.7; extra == 'dev'
Requires-Dist: matplotlib; extra == 'dev'
Requires-Dist: pytest-cov<5,>=4.1; extra == 'dev'
Requires-Dist: pytest-mock<4,>=3.12; extra == 'dev'
Requires-Dist: pytest<9,>=8.0; extra == 'dev'
Provides-Extra: jax
Requires-Dist: jax<1,>=0.7; extra == 'jax'
Provides-Extra: torch
Requires-Dist: torch<3,>=2.0; extra == 'torch'
Description-Content-Type: text/markdown

# PTED: Permutation Test using the Energy Distance

![PyPI - Version](https://img.shields.io/pypi/v/pted?style=flat-square)
[![CI](https://github.com/ConnorStoneAstro/pted/actions/workflows/ci.yml/badge.svg)](https://github.com/ConnorStoneAstro/pted/actions/workflows/ci.yml)
[![Code style: black](https://img.shields.io/badge/code%20style-black-000000.svg)](https://github.com/psf/black)
![PyPI - Downloads](https://img.shields.io/pypi/dm/pted)
[![codecov](https://codecov.io/gh/ConnorStoneAstro/pted/graph/badge.svg?token=5LISJ5BN17)](https://codecov.io/gh/ConnorStoneAstro/pted)
[![DOI](https://zenodo.org/badge/966938845.svg)](https://doi.org/10.5281/zenodo.15353928)

PTED (pronounced "ted") takes in `x` and `y` two datasets and determines if they
were sampled from the same underlying distribution. It produces a p-value under
the null hypothesis that they are sampled from the same distribution. The
samples may be multi-dimensional, and the p-value is "exact" meaning it has a
correctly calibrated type I error rate regardless of the data distribution.

![pted logo](media/pted_logo.png)

> **NOTE:** PTED is basically just the Maximum Mean Discrepancy (MMD) except using the Euclidean distance as the kernel. It took me a while to realize this, and in the meantime I made a nice python package. So I consider it a win for everyone. 

## Install

To install PTED, run the following:

```bash
pip install pted
```

If you want to run PTED on GPUs using PyTorch, then also install torch:

```bash
pip install torch
```

If you want to use JAX arrays as inputs, then also install jax:

```bash
pip install jax
```

The two functions are ``pted.pted`` and ``pted.pted_coverage_test``. For
information about each argument, just use ``help(pted.pted)`` or
``help(pted.pted_coverage_test)``. 

## What does PTED do?

You can think of it like a multi-dimensional KS-test! Although it works entirely
differently from the KS-test, this gives you some idea of how useful it is! It
is used for two sample testing and posterior coverage tests. In some cases it is
even more sensitive than the KS-test, but likely not all cases.

PTED is useful for:

* "were these two samples drawn from the same distribution?" this works even with noise, so long as the noise distribution is also the same for each sample
* Evaluate the coverage of a posterior sampling procedure, and check over/under-confidence
* Check for MCMC chain convergence. Split the chain in half or take two chains, that's two samples that PTED can work with (PTED assumes samples are independent, make sure to thin your chain accordingly!)
* Evaluate the performance of a generative ML model. PTED is powerful here as it can detect overfitting to the training sample (ensure `two_tailed = True` to check this).
* Evaluate if a simulator generates true "data-like" samples
* PTED (or just the energy distance) can be a distance metric for Approximate Bayesian Computing posteriors
* Check for drift in a time series, comparing samples before/after some cutoff time
* The energy distance is differentiable, so it can be used to train machine learning algorithms

And much more!

## Example: Two-Sample-Test

```python
from pted import pted
import numpy as np

x = np.random.normal(size = (500, 10)) # (n_samples_x, n_dimensions)
y = np.random.normal(size = (400, 10)) # (n_samples_y, n_dimensions)

p_value = pted(x, y)
print(f"p-value: {p_value:.3f}") # expect uniform random from 0-1
```

## Example: Coverage Test

```python
from pted import pted_coverage_test
import numpy as np

g = np.random.normal(size = (100, 10)) # ground truth (n_simulations, n_dimensions)
s = np.random.normal(size = (200, 100, 10)) # posterior samples (n_samples, n_simulations, n_dimensions)

p_value = pted_coverage_test(g, s)
print(f"p-value: {p_value:.3f}") # expect uniform random from 0-1
```

Note, you can also provide a filename via a parameter: `sbc_histogram = "sbc_hist.pdf"` and this will generate an SBC histogram from the test[^1].

You can also generate a Probability Integral Transform (PIT) plot via `pit_plot = "pit.pdf"` for `pted_coverage_test`. The PIT plot shows the empirical CDF of the p-values against the expected uniform CDF (1:1 diagonal), along with a shaded band determined by `pit_confidence` (95% by default). The band is simultaneous: under the null the whole curve stays inside it with probability `pit_confidence`, so any portion of the curve that pokes out is evidence of non-uniform p-values at significance `1 - pit_confidence`, indicating a potentially miscalibrated posterior.

The band is built from the exact null distribution of the p-values. The ECDF can only move where the p-values can actually land, so rather than constraining order statistics the band constrains the counts there.
Working with counts rather than order statistics is what makes this usable here. Permutation p-values are discrete — after enumeration they are exactly uniform on a grid of `L = len(permute) + 1` points — and ties are precisely what a binomial count expects.
`pted_coverage_test` works out the lattice size and passes it in for you.
Ultimately, this means that the fact we only use a finite number of simulations to compute the coverage test does not affect the validity of the bands in the PIT plot; any point that sticks out of the grey shaded area indicates possible miscalibration at the given significance level (default 95%).

## How it works

### Two sample test

PTED uses the energy distance of the two samples `x` and `y`, this is computed as:

$$d = \frac{2}{n_xn_y}\sum_{i,j}||x_i - y_j|| - \frac{1}{n_x(n_x-1)}\sum_{i\neq j}||x_i - x_j|| - \frac{1}{n_y(n_y-1)}\sum_{i\neq j}||y_i - y_j||$$

The within-group sums run over distinct pairs, skipping the `i = j` terms which
are zero by construction. That makes `d` an unbiased estimator of the population
energy distance, so under the null it scatters around zero and can come out
slightly negative (and strongly negative if `x` and `y` are *more* alike than
chance allows, e.g. if they share samples). Only its rank among the permuted
values matters. For equal sample sizes this ranks the permutations identically
to the `1/n^2` normalisation of Székely & Rizzo, so the p-value is unchanged; for
unequal sizes the two differ slightly.

The energy distance measures distances between pairs of points[^2]. It becomes more
positive if the `x` and `y` samples tend to be further from each other than from
themselves. We demonstrate this in the figure below, where the `x` samples are
drawn from a (thick) circle, while the `y` samples are drawn from a (thick)
line.

![pted demo test](media/test_PTED.png)

In the left figure, we show the two distributions, which by eye are clearly not
drawn from the same distribution (circle and line). In the center figure we show
the individual distance measurements as histograms. To compute the energy
distance, we would sum all the elements in these histograms rather than binning
them. You can also see a schematic of the distance matrix, which represents
every pair of samples and is colour coded the same as the histograms. In the
right figure we show the energy distance as a vertical line, the grey
distribution is explained below.

The next element of PTED is the permutation test. For this we combine the `x`
and `y` samples into a single collection `z`. We then randomly shuffle (permute)
the `z` collection and break it back into `x` and `y`, now with samples randomly
swapped between the two distributions (though they are the same size as before).
If we compute the energy distance again, we will get very different results.
This time we are sure that the null hypothesis is true, `x` and `y` have been
drawn from the same distribution (`z`), and so the energy distance will be quite
low. If we do this many times and track the permuted energy distances we get a
distribution, this is the grey distribution in the right figure. Below we show
an example of what this looks like.

![pted demo permute](media/permute_PTED.png)

Here we see the `x` and `y` samples have been scrambled in the left figure. In
the center figure we see the components of the energy distance matrix are much
more consistent because `x` and `y` now follow the same distribution (a mixture
of the original circle and line distribution). In the right figure we now see
that the vertical line is situated well within the grey distribution. Indeed the
grey distribution is just a histogram of many re-runs of this procedure. We
compute a p-value by taking the fraction of the energy distances that are
greater than the current one.

### Coverage test

In the coverage test we have some number of simulations ``nsim`` where there is
a true value ``g`` and some posterior samples ``s``. The procedure goes like
this, first you sample from your prior: ``g ~ Prior(G)``. Then you sample from
your likelihood: ``x ~ Likelihood(X | g)``. Then you sample from your posterior:
``s ~ Posterior(S | x)``, you will want many samples ``s``. You repeat this
procedure ``nsim`` times. The ``g`` and ``s`` samples are what you need for the
test.

Internally, for each simulation separately we use PTED to compute a p-value,
essentially asking the question "was ``g`` drawn from the distribution that
generated ``s``?". Individually, these tests are possibly not especially
informative (unless the sampler is really bad), however their p-values must have
been drawn from ``U(0,1)`` under the null-hypothesis[^3]. Thus we just need a
way to combine their statistical power. It turns out that for some ``p ~
U(0,1)``, we have that ``- 2 ln(p)`` is chi2 distributed with ``dof = 2``. This
means that we can sum the chi2 values for the PTED test on each simulation and
compare with a chi2 distribution with ``dof = 2 * nsim``. We use a simple
doubling procedure (``2 * min(p_right, p_left)``) to get the p-value meaning
that if your posterior is underconfident or overconfident, you will get a small
p-value that can be used to reject the null.

## Example: Containment Test

`pted` asks "are these the same distribution?". `pted_containment_test` asks the
weaker, directional question "**does y cover x?**" — a sample drawn from a
tighter distribution than `y` passes, while one that spreads beyond it, sits off
to one side, or throws a few points clear of it fails.

```python
from pted import pted_containment_test
import numpy as np

y = np.random.normal(size = (500, 10))
inside  = np.random.normal(size = (100, 10)) * 0.5
outside = np.random.normal(size = (100, 10)) * 1.6

print(pted_containment_test(inside,  y))  # large: contained
print(pted_containment_test(outside, y))  # small: not contained
```

Unlike every other test in the package this one is **not symmetric** — swapping
the arguments asks the other question and will usually give a different answer.
That asymmetry is the point, and it is why the energy distance alone cannot
answer it: the energy distance is just as large when `x` is *tighter* than `y`
as when it is broader, so a one-tailed energy test rejects 99.8% of the time on
a perfectly contained sample.

Instead, every `x` point gets a **depth** — its mean distance to the points
labelled `y` — so peripheral points score high. This is nearly identical to
`pted_coverage_test`, and we use the same `-2 sum log p` formula to combine the
p-values for each `x`. However, instead of interpreting the result as a $\chi^2$
distribution, we use a secondary permutation test to calibrate it. 

The test is exact where `x` and `y` share a distribution and conservative inside
the null, which is the goal for a containment test. It measures depth rather
than literal support: a tight cluster of `x` sitting in a low-density pocket
well inside `y` counts as contained. So use a heavy dose of caution when
interpreting the results.

### Read the plot as well as the p-value

`-2 sum log p` is a sum whose per-point floor is zero against a null mean of
two, so a bulk of `x` sitting deep inside `y` banks slack that can hide a
handful of points `y` cannot reach at all. In testing, 5% of `x` placed five
sigma outside the prior predictive went undetected, while every other failure
mode — data broader than `y`, data offset into its tail — was caught cleanly. So
pass `pit_plot`:

```python
p = pted_containment_test(data, prior_predictive, pit_plot = "containment.pdf")
```

The plot shows the empirical CDF of the per-point depth p-values — one step per
point of `x` — on a logarithmic p axis, because the whole diagnostic lives at
the left edge. It carries two reference marks and no diagonal.

The **upper bound** is the one-sided simultaneous ceiling for `n1` p-values that
really are uniform: the same as the ordinary PIT plot, with the lower edge
dropped. Only the upper edge means anything here. Any part of the CDF below this
line indicates the `x` values are likely contained at the threshold level (95%
by default).

The **vertical line** marks threshold p-value (0.05 by default). Parts of the
CDF to the right of this line are embedded in `y` at the threshold level and so
likely are contained. Points to its left are peripheral relative to `y` and so
this may suggest the values are not contained. However, when many values of `x`
are being tested some leakage to low p-values are expected, which is why the
uniform-upper-bound line is also plotted.

Read them together rather than as a decision rule. A curve with points to the
left of the vertical line and above the uniform-upper-bound is a fair indication
that `x` is not contained in `y`. Though results either way are not definitive,
as is the nature of null hypothesis testing, and even moreso here.

## Interpreting the results

### Two sample test

This is a null hypothesis test, thus we are specifically asking the question:
"if `x` and `y` were drawn from the same distribution, how likely am I to have
observed an energy distance as extreme as this?" This is fundamentally different
from the question "how likely is it that `x` and `y` were drawn from the same
distribution?" Which is really what we would like to ask, but I am unaware of
how we would do that in a meaningful way. It is also important to note that we
are specifically looking at extreme energy distances, so we are not even really
talking about the probability densities directly. If there was a transformation
between `x` and `y` that the energy distance was insensitive to, then the two
distributions could potentially be arbitrarily different without PTED
identifying it. For example, since the default energy distance is computed with
the Euclidean distance, a single dimension in which the values are orders of
magnitude larger than the others could make it so that all other dimensions are
ignored and could be very different. For this reason we suggest using the metric
`mahalanobis` if this is a potential issue in your data.

### Coverage Test

For the coverage test we apply the PTED two sample test to each simulation
separately. We then combine the resulting p-values using chi squared where the
resulting degrees of freedom is 2 times the number of simulations. Because of
this, we can detect underconfidence or overconfidence. Specifically we detect
the average over/under confidence, it is possible to be overconfident in some
parts of the posterior and underconfident in others. Underconfidence is when the
posterior distribution is too large, it covers the ground truth by spreading too
thin and not fully exploiting the information in the prior/likelihood of the
posterior sampling process. Sometimes this is acceptable/expected, for example
when using Approximate Bayesian Computation one expects the posterior to be at
least slightly underconfident. Overconfidence is when the posterior is too
narrow and so the ground truth appears as an outlier from its perspective. This
can occur in two main ways, one is by having a too narrow posterior. This could
occur if the measurement uncertainty estimates were too low or there were
sources of error not accounted for in the model. Another way is if your
posterior is biased, you may have an appropriately broad posterior, but it is in
the wrong part of your parameter space. PTED has no way to distinguish these or
other modes of overconfidence, however just knowing under/over-confidence can be
powerful. As such, by default the PTED coverage test will warn users as to which
kind of failure mode they are in if the `warn_confidence` parameter is not
`None` (default is 1e-3).

### Necessary but not Sufficient

PTED is a null hypothesis test. This means we assume the null hypothesis is true
and compute a probability for how likely we are to have a pair of datasets with
a certain energy distance. If PTED gives a very low p-value then it is probably
safe to reject that null hypothesis (at the significance given by the p-value).
However, if the p-value is high and you cannot reject the null, then that does
not mean the two samples were drawn from the same distribution! Merely that PTED
could not find any significant discrepancies. The samples could have been drawn
from the same distribution, or PTED could be insensitive to the deviation, or
maybe the test needs more samples. In some sense PTED (like all null hypothesis
tests) is "necessary but not sufficient" in that failing the test is bad news
for the null, but passing the test is possibly inconclusive[^5]. Use your judgement,
and contact me or some smarter stat-oriented person if you are unsure about the
results you are getting!

## Arguments

### Two Sample Test

```python
def pted(
    x: Union[np.ndarray, "Tensor", "jax.Array"],
    y: Union[np.ndarray, "Tensor", "jax.Array"],
    permutations: int = 1000,
    return_all: bool = False,
    n_landmarks: Optional[int] = None,
    two_tailed: bool = True,
    prog_bar: bool = False,
    batch_size: Optional[int] = None,
    rng=None,
) -> Union[float, tuple[float, np.ndarray, float]]:
```

* **x** *(Union[np.ndarray, Tensor, jax.Array])*: first set of samples. Shape (N, *D)
* **y** *(Union[np.ndarray, Tensor, jax.Array])*: second set of samples. Shape (M, *D)
* **permutations** *(int)*: number of permutations to run. This determines how accurately the p-value is computed.
* **return_all** *(bool)*: if True, return the test statistic and the permuted statistics with the p-value. If False, just return the p-value. bool (default: False)
* **n_landmarks** *(Optional[int])*: if not None, estimate the energy distance from a rectangular distance matrix instead of the full pairwise matrix. `n_landmarks` points are drawn from the pooled sample as "landmarks", and only the distance from every sample to each landmark is computed, so the cost drops from `O(n^2 d)` to `O(n m d)` for `m = n_landmarks`. A value covering the whole pooled sample, or None, runs the exact full-matrix computation.
* **two_tailed** *(bool)*: if True, compute a two-tailed p-value. This is useful if you want to reject the null hypothesis when x and y are either too similar or too different. If False, only checks for dissimilarity but is more sensitive. Default is True.
* **prog_bar** *(bool)*: if True, show a progress bar to track the progress of permutation tests. Default is False.
* **batch_size** *(Optional[int])*: number of permutations evaluated per matrix product. Larger values are faster (especially on GPU) at `O(batch_size * n)` extra memory. None picks the largest size that keeps the batch-scaling tensors within 1 GiB and bounds one batch's matrix product, so a run cannot allocate unboundedly whatever the sample shape.
* **rng**: seed, `np.random.Generator`, or None to draw from the global numpy state, so `np.random.seed` still controls reproducibility.

### Coverage test

```python
def pted_coverage_test(
    g: Union[np.ndarray, "Tensor", "jax.Array"],
    s: Union[np.ndarray, "Tensor", "jax.Array"],
    permutations: int = 1000,
    warn_confidence: Optional[float] = 1e-3,
    return_all: bool = False,
    n_landmarks: Optional[int] = None,
    sbc_histogram: Optional[str] = None,
    sbc_bins: Optional[int] = None,
    pit_plot: Optional[str] = None,
    pit_confidence: float = 0.95,
    prog_bar: bool = False,
    batch_size: Optional[int] = None,
    rng=None,
) -> Union[float, tuple[np.ndarray, np.ndarray, float]]:
```

* **g** *(Union[np.ndarray, Tensor, jax.Array])*: Ground truth samples. Shape (n_sims, *D)
* **s** *(Union[np.ndarray, Tensor, jax.Array])*: Posterior samples. Shape (n_samples, n_sims, *D)
* **permutations** *(int)*: number of permutations to run. This determines how accurately the p-value is computed.
* **return_all** *(bool)*: if True, return the test statistic and the permuted statistics with the p-value. If False, just return the p-value. bool (default: False)
* **n_landmarks** *(Optional[int])*: if not None, estimate the energy distance from a rectangular distance matrix instead of the full pairwise matrix. `n_landmarks` points are drawn from the pooled sample as "landmarks", and only the distance from every sample to each landmark is computed, so the cost drops from `O(n^2 d)` to `O(n m d)` for `m = n_landmarks`. A value covering the whole pooled sample, or None, runs the exact full-matrix computation.

  Because the ground truth is a single point, the per-simulation test runs in the `singleton` regime, where the permutation subgroup reaches only `n - m` distinct label assignments. Keep `n_landmarks` well below the number of posterior samples.
* **sbc_histogram** *(Optional[str])*: If given, the path/filename to save a Simulation-Based-Calibration histogram.
* **sbc_bins** *(Optional[int])*: If given, force the histogram to have the provided number of bins. Otherwise, select an appropriate size: ~sqrt(N).
* **pit_plot** *(Optional[str])*: If given, the path/filename to save a Probability Integral Transform (PIT) plot of the per-simulation p-values against the expected uniform distribution, with a shaded KS confidence band.
* **pit_confidence** *(float)*: Confidence level for the PIT plot's simultaneous band. Default is 0.95 (95%). Only used when `pit_plot` is not None.
* **prog_bar** *(bool)*: if True, show a progress bar to track the progress of simulations. Default is False.
* **batch_size** *(Optional[int])*: number of permutations evaluated per matrix product. Larger values are faster (especially on GPU) at `O(batch_size * n)` extra memory. None picks the largest size that keeps the batch-scaling tensors within 1 GiB and bounds one batch's matrix product, so a run cannot allocate unboundedly whatever the sample shape.
* **rng**: seed, `np.random.Generator`, or None to draw from the global numpy state, so `np.random.seed` still controls reproducibility.

## GPU Compatibility

PTED works on both CPU and GPU. All that is needed is to pass the `x` and `y` as
PyTorch Tensors or JAX Arrays on the appropriate device.

Example with PyTorch:
```python
from pted import pted
import numpy as np
import torch

x = np.random.normal(size = (500, 10)) # (n_samples_x, n_dimensions)
y = np.random.normal(size = (400, 10)) # (n_samples_y, n_dimensions)

p_value = pted(torch.tensor(x), torch.tensor(y))
print(f"p-value: {p_value:.3f}") # expect uniform random from 0-1
```

Example with JAX:
```python
from pted import pted
import numpy as np
import jax.numpy as jnp

x = np.random.normal(size = (500, 10)) # (n_samples_x, n_dimensions)
y = np.random.normal(size = (400, 10)) # (n_samples_y, n_dimensions)

p_value = pted(jnp.array(x), jnp.array(y))
print(f"p-value: {p_value:.3f}") # expect uniform random from 0-1
```

## Memory and Compute limitations

If a GPU isn't enough to get PTED running fast enough for you, or if you are
running into memory limitations, there are still options! We can use an
approximation of the energy distance; the test stays exact, it just becomes less
sensitive than it would otherwise be. Instead of building the full `n x n`
pairwise distance matrix (`n = n_samp_x + n_samp_y`), PTED builds a smaller
rectangular `n x m` matrix of distances from every sample to `m` "landmark"
points drawn from the pooled sample — the Nyström-style subsampling familiar
from kernel methods. Set `n_landmarks = m` to choose how many.

Building the full matrix costs `O(n^2 d)` and each permutation `O(n^2)`; the
rectangular matrix costs `O(n m d)` and each permutation `O(n m)`. So for a
fixed number of landmarks PTED grows linearly with dataset size, much like other
large scale (machine learning oriented) two sample tests. Permutations are
evaluated in batches as a single matrix product rather than one at a time, and
everything heavy stays on whichever backend your arrays live on.

Example:
```python
from pted import pted
import numpy as np

x = np.random.normal(size = (5000, 10)) # (n_samples_x, n_dimensions)
y = np.random.normal(size = (4000, 10)) # (n_samples_y, n_dimensions)

p_value = pted(x, y, n_landmarks = 200)
print(f"p-value: {p_value:.3f}") # expect uniform random from 0-1
```

### Enumerating small permutation groups

If the permutation subgroup turns out to have no more than `permutations`
members, PTED evaluates *all* of it rather than sampling. Sampling a small
group keeps redrawing the observed labelling, and each of those ties inflates
the p-value under the `>=` convention: with a group of 151 and 199 draws,
`P(p <= 0.02)` came out at 0.0169 instead of 0.0200. Walking the group instead
puts it at 0.0195, leaving only the unavoidable granularity of a 151-point
lattice — and it costs *less* compute, since 151 evaluations beat 199. This
mostly matters for `pted_coverage_test`, where the ground truth is a single
point and the group is only `n - m` large.

One consequence: when this kicks in, `return_all` gives back
`reference_size - 1` permuted statistics rather than `permutations` of them.
The p-value is then exact rather than sampled.

### Why using landmarks is still exact

The subtlety is that the landmarks cannot simply be re-split between the two
groups after each shuffle. The landmark set `L` is chosen using the group labels
(so that both groups are represented), which means the labels are no longer
uniformly distributed once you condition on `L`. PTED handles this by confining
permutations to the subgroup that fixes `L`: labels are shuffled within the
landmark positions and within their complement, never between. `L` then stays put
under every permutation and the observed labelling is exchangeable with the permuted
ones. That is what makes the p-value exact. (`L` is chosen from sample positions
and group labels only, never from the data values — picking landmarks by maximin,
k-means or leverage would break the argument. "Landmark" here means only "a point
everything is measured against", never "a point chosen for its importance".)

What you give up is sensitivity: the null spread grows like `sqrt(n / m)`, so the
smallest detectable energy distance scales as `(n m)^-0.5` rather than `n^-1`.
The detection threshold degrades as one over the square root of the compute.

You also give up some p-value resolution, and this is most significant when one
group holds a single point — which is exactly the per-simulation test inside
`pted_coverage_test`. The lone point is kept out of `L` — it has no within-group
pairs, so landmarking it would buy nothing — and its label then roams over the
`n - m` positions outside, which is the whole reachable set. The smallest
attainable p-value is about `1 / (n - m)` however many permutations you draw, so
PTED raises a `PermutationResolutionWarning` when that set is too small to
resolve the p-value you asked for.

This is the sharpest case of a general rule: **landmarks are for `m << n`.** At
`m = n / 2` the rectangular matrix saves only a factor of two over the exact
test, which does not pay for the sensitivity and resolution it costs. Reach for
them when the full matrix will not fit or will not finish, not to shave a
constant factor.

### Does restricting the permutations still measure the energy distance?

Yes — the restriction is on how the statistic is *calibrated*, not on what it
measures. Every block mean still averages genuine distances between genuine
members of the two groups: the cross term pairs **all** `n_x` samples against the
large-group landmarks and **all** `n_y` samples against the small-group
landmarks, so no sample is demoted to a mere marker. Each block mean is an unbiased estimate of
exactly the population quantity the full test estimates, which makes the whole
statistic an unbiased estimator of the population energy distance — an
*incomplete U-statistic* in the sense of Janson (1984), averaging over a subset
of the pairs rather than computing a different function of them.

Freezing the per-group landmark counts is a restricted-randomisation device, the
same idea as conditioning on the margins in Fisher's exact test: it removes a
nuisance source of variability from the null rather than changing the estimand.
Empirically it is power-neutral. Compare it against the obvious alternative —
draw `L` uniformly by position, so that `L` is independent of the labels and the
full permutation group is legal — and both reject at the same rate (`n = 200`,
`m = 40`, 500 trials: 0.16 vs 0.16 at a 0.25σ shift, 0.31 vs 0.35 at 0.40σ). The
power lost relative to the full test comes from subsampling, not from the
subgroup.

What the subgroup buys is the freedom to choose `L` *by design*. Drawing `L`
blind to the labels leaves the smaller group with no landmarks at all — and its
within-group term unestimable — distressingly often once the samples are
unbalanced: with `n_x = 3`, `n_y = 200` and `m = 40`, 52% of permutations have no
x-landmarks. Choosing `L` so that it covers the small group fixes that, and the
subgroup is what keeps the test exact when you do.

## Citation

If you use PTED in your work, please include a citation to the [zenodo
record](https://doi.org/10.5281/zenodo.15353928) and also see below for
references of the underlying method.

## Reference

I didn't invent this test, I just think its neat. Here is a paper on the subject:

```
@article{szekely2004testing,
      title = {Testing for equal distributions in high dimension},
     author = {Sz{\'e}kely, G{\'a}bor J and Rizzo, Maria L and others},
    journal = {InterStat},
     volume = {5},
     number = {16.10},
      pages = {1249--1272},
       year = {2004},
  publisher = {Citeseer}
}
```

Permutation tests are a whole class of tests, with much literature. Here are
some starting points:

```
@book{good2013permutation,
  title={Permutation tests: a practical guide to resampling methods for testing hypotheses},
  author={Good, Phillip},
  year={2013},
  publisher={Springer Science \& Business Media}
}
```

```
@book{rizzo2019statistical,
  title={Statistical computing with R},
  author={Rizzo, Maria L},
  year={2019},
  publisher={Chapman and Hall/CRC}
}
```

There is also [the wikipedia
page](https://en.wikipedia.org/wiki/Permutation_test), and the more general
[scipy
implementation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.permutation_test.html),
and other [python implementations](https://github.com/qbarthelemy/PyPermut)

As for the posterior coverage testing, this is also an established technique.
See the references below for the nitty gritty details and to search further look for "Simulation-Based Calibration".

```
@article{Cook2006,
     title = {Validation of Software for Bayesian Models Using Posterior Quantiles},
    author = {Samantha R. Cook and Andrew Gelman and Donald B. Rubin},
   journal = {Journal of Computational and Graphical Statistics},
      year = {2006}
 publisher = {[American Statistical Association, Taylor & Francis, Ltd., Institute of Mathematical Statistics, Interface Foundation of America]},
       URL = {http://www.jstor.org/stable/27594203},
   urldate = {2026-01-09},
    number = {3},
    volume = {15},
     pages = {675--692},
      ISSN = {10618600},
}
```

```
@ARTICLE{Talts2018,
       author = {{Talts}, Sean and {Betancourt}, Michael and {Simpson}, Daniel and {Vehtari}, Aki and {Gelman}, Andrew},
        title = "{Validating Bayesian Inference Algorithms with Simulation-Based Calibration}",
      journal = {arXiv e-prints},
     keywords = {Statistics - Methodology},
         year = 2018,
        month = apr,
          eid = {arXiv:1804.06788},
        pages = {arXiv:1804.06788},
          doi = {10.48550/arXiv.1804.06788},
archivePrefix = {arXiv},
       eprint = {1804.06788},
 primaryClass = {stat.ME},
}
```

If you think those are neat, then you'll probably also like this paper, which uses HDP regions and a KS-test. It has the same feel as PTED but works differently, so the two are complimentary.

```
@article{Harrison2015,
  author = {Harrison, Diana and Sutton, David and Carvalho, Pedro and Hobson, Michael},
   title = {Validation of Bayesian posterior distributions using a multidimensional Kolmogorov–Smirnov test},
 journal = {Monthly Notices of the Royal Astronomical Society},
  volume = {451},
  number = {3},
   pages = {2610-2624},
    year = {2015},
   month = {06},
    issn = {0035-8711},
     doi = {10.1093/mnras/stv1110},
     url = {https://doi.org/10.1093/mnras/stv1110},
  eprint = {https://academic.oup.com/mnras/article-pdf/451/3/2610/4011597/stv1110.pdf},
}
```

[^1]: See the Simulation-Based Calibration paper by Talts et al. 2018 for what "SBC" is.
[^2]: Yes, I see it now. That formula is the same as for MMD. Just think of it as a special case of MMD that uses the Euclidean distance kernel. In fact this way it is often more sensitive than the RBF kernel typically used for MMD.
[^3]: Since PTED works by a permutation test, we only get the p-value from a discrete uniform distribution. By default we use 1000 permutations, if you are running an especially sensitive test you may need more permutations, but for most purposes this is sufficient.
[^4]: Actually, we take (q + 1) / (Nsamp + 1) rather than q/Nsamp where q is the number of posterior samples with posterior density greater than the ground truth. This just turns out to be a better estimator for finite Nsamp.
[^5]: actual "necessary but not sufficient" conditions are a different thing than null hypothesis tests, but they have a similar intuitive meaning.