2 Dimensional Kolmogorov-Smirnov test in astropy.stats (and 1D while we are at it)

### Close this if implemented upstream

* https://github.com/scipy/scipy/issues/17514

### Description

**Kolmogorov-Smirnov test (KS test)** can be used as a statistical test of whether an observed set of data is drawn from an expected distribution (or the distribution defined by another set of data), and is widely used in Astronomy. The 1D KS test is provided in various forms by `scipy.stats`. However, currently, there are no maintained or released packages that offer the 2D or multi-D version, (although code for it can be found in various forms on GitHub).

There are two seminal Astronomy papers—with nearly 800 citations among them—that developed the 2D and Multi-D approach for the KS test. Majority of their citations come form the past 10 years so this is increasingly useful for data-driven Astronomy.
- [Peacock J. A. 1983](https://ui.adsabs.harvard.edu/abs/1983MNRAS.202..615P/abstract): Initial 2D idea
- [Fasano & Franceschini 1987](https://ui.adsabs.harvard.edu/abs/1987MNRAS.225..155F/abstract): Multi-D extension

**I propose** that, at a minimum, the 2D KS test from the paper be made available via `astropy.stats`. Additionally, if the code is developed in a nicely abstracted manner, then it should be easy to add multi-D KS tests. The 1D one is in `scipy.stats`, and I know `astropy.stats` does not supplant it, but if people see a need, the scipy.stats version can be wrapped to resemble the 2D and multi-D version to provide a common interface, or it can be used in our unit tests.

The basic idea of the algorithm is to calculate the 2-D or N-D average maximum distance between two input sets of random variables (or between 1 and a distribution).

My initial idea of the API is that there are 1 to 3 `KSTest` classes representing the logically distinct algorithms of the 2D, and multi-D cases (and perhaps 1D). Additionally, there probably should be a functional interface for using them, something like `kstest_2d`, `kstest_nd`, (and `kstest_1d)`, to match the functional style of `astropy.stats`. Or there can be a single KSTest class, which is N-dimensional by default, but can fall back to simpler algorithms for the 2D and 1D versions.

Regardless of what's chosen, the intended usecase is to pass in datasets, and be able to get the two-sided KS statistic and probability using some chosen method. A skeleton:

```Python
from dataclasses import dataclass, field # just for nicer formatting.
from enum import Enum
from typing import Any, Callable

@dataclass # just for nicer formatting.
class KSTest:
dataset1_Xi: array_like # 1, 2 or N dimensional Array
dataset2_Xi: array_like
method: str | callable = 'asymptotic' # alternatively provide a callable.
method_kwargs: dict[str, Any] = field(default_factory=dict) # things like iteration limit, precision.
transform: callable = field(init=False)
dimensions: int = field(init=False)

def __post_init__(self):
if isinstance(self.method, str):
# will raise meaningful looking error if method not in Methods:
self.transform = getattr(self, Methods(self.method))
elif isinstance(self.method, Callable):
self.transform = self.method
else:
raise ValueError("method must be correct, etc.")
self.dimensions = self.dataset1_Xi.shape[1]

def calc_ks(self) -> KSStatistic:
# return the KS statistics D+/- using some algorithm.
if self.dimensions == 2:
return KSStatistic(dplus = calc_dplus(), dminus = calc_dminus())

def calc_prob(self, d: KSStatistic) -> float:
"""calculate a p-value for a given d, using the properties of the input datasets, and via given method"""
# different methods have different transforms on D, so may need other things from self such as len.
return self.transform(d, **self.method_args)
...

# @staticmethod # not a static method so we can use properties of input X1, X2, if relevant.
def asymptotic(self, d: KSStatistic, maxiter: int, precision: float) -> float:
... # calculates the probability with the asymptotic kuiper method for a given statistic. See Numerical Methods in C.

@dataclass
KSStatistic:
"""Some methods may require dplus, dminus, or other measures than just D to calculate probability."""
d: float = field(init=False)
dplus: float
dminus: float

__post_init__(self) ->:
self.d = max(self.dplus, self.dminus)

# This part below is just trying to structure the abstraction, not needed.
class Methods(StrEnum): # only in py3.11 but easy to create from Enum
asymptotic = 'asymptotic'
exact = 'exact'
other = ''

_Methods[

# Just one way of making a str enumeration
class StrEnum(str, enum.Enum):
"""Enum with string values."""

def __str__(self) -> str:
return str(self.value)
```

A KS statistic and probability can be derived. I would propose using Kuiper's statistic, an asymptotic variant of the KS test that can be robustly calculated and is itself robust to differences not just in the median but among the whole distribution. (The KS test basically compares maximum distances between 2 CDFs, and is invariant to simple coordinate transformations, so by wrapping the CDF into a circle, we can choose any part the be near the middle, and we get a result that is robust to differences throughout the CDF. Anyway, standard numerical approaches to this can be found in Numerical Recipes textbook or in `scipy` itself, as well as in papers listed by `scipy`.

Whether this is the publicly exposed class or not, the basic workflow is to find the "average maximum distance" between the CDFs, the so-called KS statistic, then do a test on that statistic to give a probability which can be converted into one to reject the null hypothesis that the CDFs come from the same distribution. Some sort of p-value. Like I said, `scipy.stats` already does this for the 1D case. The class would need to be able to at a minimum give the KS statistic and its p-value, via some method of calculating it, such as the asymptotic variant.

if we want to follow the `scipy.stats` API, then the input would also need to be possibly a callable the provides a CDF, and args to pass into a distribution function to provide said CDF. `scipy` already has helper methods to do this. I'll leave that up to future extension though. It would use the above class behind the scenes.

```Python
@dataclass # just for nicer formatting.
class KSTest:
dataset1_Xi: array_like | callable # 1, 2 or N dimensional Array, or callable to produce it.
dataset2_Xi: array_like | callable # Scipy also allows using a string to represent a distribution from scipy.stats, but this may not generalize to 2D/ND.
args: # To pass into the callables. I would push towards using `functools.partial` but this is `scipy` API.
N: int # sample size to draw from the first callable
```

### Additional context

I see there are existing approaches with open source licensed code. I've collected a few of them, and am willing to provide the basic algorithm/class contribution for the 2D case. I would appreciate comments from `astropy.stats` people and any usecases that I might not be thinking of.

Start by reviewing the proposed work in astropy.stats and the linked SciPy issue. Compare the possible 1D, 2D, and N-dimensional APIs, methods, and callable inputs before narrowing the scope. Done should be defined as an agreed interface that returns the requested KS statistic and probability, with tests and documentation identified.

Written by the indexing model from the issue text.