# Discovering mesoscopic descriptions of collective movement
with neural stochastic modelling

Utkarsh Pratiush  Arshed Nabeel  Vishwesha Guttal  Prathosh AP

###### Abstract

Collective motion is an ubiquitous phenomenon in nature, inspiring engineers,
physicists and mathematicians to develop mathematical models and bio-inspired
designs. Collective motion at small to medium group sizes ($\sim$10-1000
individuals, also called the _‘mesoscale’_), can show nontrivial features due
to stochasticity. Therefore, characterizing both the deterministic and
stochastic aspects of the dynamics is crucial in the study of mesoscale
collective phenomena. Here, we use a physics-inspired, neural-network based
approach to characterize the stochastic group dynamics of interacting
individuals, through a _stochastic differential equation (SDE)_ that governs
the collective dynamics of the group. We apply this technique on both
synthetic and real-world datasets, and identify the deterministic and
stochastic aspects of the dynamics using _drift_ and _diffusion_ fields,
enabling us to make novel inferences about the nature of order in these
systems.

Machine Learning, ICML

##  1 Introduction

Collective motion is a phenomenon that observed in natural and synthetic
systems, and has fascinated physicists and biologists alike (Vicsek &
Zafeiris, 2012; Sumpter, 2010; Camazine et al., 2020). Many systems across
scales—such as microscopic organisms (Dinet et al., 2021; Be’er & Ariel,
2019), cells (Rørth, 2009; Alert & Trepat, 2020), human crowds (Chen et al.,
2018b), and synthetic active matter (Ramaswamy, 2010, 2017)—exhibit self-
organized collective movement. How the seemingly simple behaviour and
interactions of the individuals give rise to the complex self-organised
emergent dynamics of the group, is one of the central questions in the study
of collective dynamics.

At the level of individual organisms, animal behaviour is complex, and
modelling every aspect of individual stochastic animal behaviour seems an
unattainable goal. While the stochastic effects typically average out in the
limit of infinite (or sufficiently large) group sizes, real-world animal
groups are finite, and often small to medium sized (10 to 1000 individuals).
At these _‘meoscopic’_ scales, the individual-level stochasticity can affect
on the group dynamics in non-trivial ways. Therefore, a correct description of
the dynamics at the mesoscale should incorporate stochasticity (McKane &
Newman, 2004; Yates et al., 2009; Biancalani et al., 2014; Brückner et al.,
2019; Jhawar et al., 2020).

Owing to the complexity at the individual behavioural level, one often uses a
coarse-grained dynamical description of the collective system. A physicist’s
way of approaching this this is to quantify the state of the group, the so
called _order parameter_ , whose dynamics is then modelled to obtain a
parsimonious description of the system. In the context of collective motion,
order parameter can be the _group polarisation_ —a quantity analogous to the
magnetisation of spin systems (Vicsek et al., 1995). To study the dynamics of
order parameter at the mesoscale, the inherent stochasticity becomes crucial;
and the order parameter dynamics must be modelled as a stochastic process. A
commonly used framework for continuous stochastic processes such as this is
the _diffusion process_ , which can be described using a _stochastic
differential equation (SDE)_. An SDE describes a stochastic process by
disentangling it into a deterministic term, called the _drift_ , and a
stochastic term, called the _diffusion_ (Van Kampen, 1992; Gardiner, 2009).

Analytically deriving coarse-grained SDE descriptions starting from
individual-level interactions is an arduous task, even for highly simplified
toy models of collective behaviour (McKane & Newman, 2004; Biancalani et al.,
2014; Van Kampen, 1992; Jhawar et al., 2019). Indeed, for real-world
collective systems, one seldom has a detailed understanding of the individual
behaviour; therefore, obtaining a group-level description seems hopeless.

We address this challenge by tackling the inverse problem: we use a data-
driven approach to estimate mesoscopic SDEs directly from observed group
trajectory data. From the individual movement trajectories, we compute the
_polarization order parameter_ , quantifying the level of alignment in the
group. We then use a physics-inspired, neural-network based approach (Dietrich
et al., 2021; Evangelou et al., 2022) to fit an SDE model to describe the
temporal dynamics of polarization. To facilitate the interpretability of the
discovered neural SDE—something which is notoriously hard in neural-network-
based approaches—we propose a method to visualize the discovered drift and
diffusion as fields. The drift and diffusion fields enable us to readily
identify the deterministic equilibria, structure of stochasticity, etc.

We use this approach to study both simulated as well as real-world fish
(species _Etroplus suratensis_) schools. For a well-studied simulation model
of collective behaviour (Buhl et al., 2006; Biancalani et al., 2014; Dyson et
al., 2015; Jhawar et al., 2020; Jhawar & Guttal, 2020), the neural SDE models
can accurately recover the mean field SDEs—and hence identify underlying
interactions. We also discover SDEs from a real world dataset of schooling
fish, where we observe the signatures of _noise-induced order_ , a counter-
intuitive phenomenon where intrinsic noise serves to increase group-level
order rather than destroy it (Biancalani et al., 2014; Jhawar & Guttal, 2020).
This was first observed by (Jhawar et al., 2020), who used the conventional
SDE discovery methods (Jhawar & Guttal, 2020; Tabar, 2019). Our neural-
network-based approach recapitulates the same result; we emphasize that this
result is non-obvious, since neural SDEs capture a much wider hypothesis class
than the conventional SDEs.

The main contributions of this study are as follows:

  * •

We demonstrate the efficacy neural-network-based approaches (Dietrich et al.,
2021; Evangelou et al., 2022) for discovering mesoscale stochastic dynamics
for collective behaviour data, using both simulated datasets where the
mesoscale SDEs are analytically established.

  * •

To make the results of the _grey-box_ neural network techniques more
interpretable, we propose novel visualization techniques to visualize the
drift and diffusion fields in an easily interpretable form.

  * •

We apply these techniques to study the mesoscale dynamics in real-world
experimental datasets of fish schools (Jhawar et al., 2020), thereby
demonstrating the significance of these techniques in animal behaviour and
ecology. To the best of our knowledge, ours is the first study that applies
these techniques to a real world dataset of animal collective behaviour, or
even the broader field of ecology.

#### Related work:

A class of techniques for estimating SDEs from data involve computing the
_jump moments_ from the observed time-series, and estimating the drift and
diffusion coefficients of the SDE using certain conditional averages of the
jump-moments (Tabar, 2019; Gardiner, 2009). However, these conditional
averages are often noisy and inaccurate, especially when the available data is
limited. More recently, approaches have been developed to use jump moments to
obtain smooth estimates of the drift and diffusion coefficients, using kernel
density estimation, sparse regression etc. (Gorjão & Meirinhos, 2019; Nabeel
et al., 2022; Boninsegna et al., 2018; Wang et al., 2022). Another approach is
to use parametric estimation techniques, by assuming an appropriate parametric
form for the drift and diffusion functions (Nielsen et al., 2000)—however,
these techniques have the drawback of having to know the parametric form of
the SDE _a priori_.

The neural-network-based approach used in this study was first introduced by
(Dietrich et al., 2021), and was subsequently applied to derive effective SDE
descriptions from simulated Brownian dynamics (Evangelou et al., 2022). This
is a _‘physics-inspired, grey-box’_ approach, where partial knowledge and
assumptions about the underlying physics is used to make the neural network
model at least partly interpretable.

Alternative approaches for data-driven neural SDE modelling also exist in
literature, such as (Song et al., 2020; Li et al., 2020), which build on
previous work on neural ODEs (Chen et al., 2018a; Liu et al., 2020).
Alternative approaches work by matching moments (O’Leary et al., 2021) or
ensemble distributions (O’Leary et al., 2021), or use VAEs to learn an SDE in
a latent space (Hasan et al., 2021). We chose the approach of (Dietrich et
al., 2021) for its simplicity and flexibility because of the way the drift and
diffusion are explicitly modelled and incorporated into a straightforward
loss-function.

##  2 Background & Methods

Figure 1: Schematic of the estimation procedure. (a) Individual trajectories
are recorded for a real or simulated collective motion. (b) From the
individual trajectories, the polarization order parameter ($\mathbf{m}$) is
computed, quantifying the degree of alignment in the school. (c) Neural
networks are trained, using a likelihood framework, to estimate the drift
($\mathbf{f}$) and diffusion ($\mathbf{g}$) functions. (d) For interpretation,
the drift and diffusion functions are visualized as vector and tensor fields
respectively.

###  2.1  Datasets

We use both simulated as well as real-world datasets of collective movement in
this study. The simulated data comes from an idealized model of collective
movement, for which the mesoscopic dynamics are well-studied (Jhawar et al.,
2020, 2019). The real-world dataset is an open-access dataset ((jiteshjhawar &
ashwinkk23, 2020)), consisting of trajectories of fish schools of different
group sizes.

#### Agent based models for collective movement.

We use simulations from an idealized agent based model of collective movement,
with $N=30$ individuals. Each individual is described solely by its velocity
vector $\mathbf{v}_{i}$. The magnitude of $\mathbf{v}_{i}$ is assumed to be
constant across all agents and across time, only the direction changes.
Direction changes can happen asynchronously at random time-points, in one of
two ways:

  1. 1.

An agent $i$ can spontaneously turn and choose a new direction, i.e.
$\theta_{i}(t)\leftarrow\eta$ for $\eta\sim\mathsf{Unif}[-\pi,+\pi)$.

  2. 2.

An agent $i$ can copy the direction of one or more other agents, chosen
randomly from the entire group.

The turn and copy events themselves happen as a Poisson process. The models
can be classified as _pairwise_ or _ternary_ interaction models based on how
many neighbours are copied. The mesoscopic SDEs for these models are already
known in literature (Jhawar et al., 2020, 2019)—also see Appendix A.
Therefore, the model simulations act as a test-bed to verify the efficacy of
the data-driven estimation procedure.

#### Real-world datasets of schooling fish.

We used an openly available dataset of fish schooling (Jhawar et al., 2020;
jiteshjhawar & ashwinkk23, 2020). This dataset consists of the position
$\mathbf{x}_{i}(t)$ and velocity $\mathbf{v}_{i}(t)$ trajectories of fish
(species _Etroplus suratensis_), swimming in a circular arena. From the
individual fish trajectories, the polarization order parameter is computed as
detailed in the following section. We used datasets with $N=15,30$ and $60$
fish. The data was available at a sampling interval of $0.12s$.

###  2.2  Mesoscopic descriptions of collective dynamics

While modelling collective movement, the emergent dynamics of the group can
often be characterized in terms of an _order parameter_ , which characterizes
the degree of order in the group. An oft-used order parameter is the group
polarization, which captures the level of alignment among the individuals of
the group.

For a group of $N$ individuals, each with positions $\mathbf{x}_{i}$ and
velocities $\mathbf{v}_{i}$, the polarization can be computed as the mean of
the normalized velocity vectors, i.e.,

$\displaystyle\mathbf{m}=\frac{1}{N}\sum_{i=1}^{N}\frac{\mathbf{v}_{i}}{|\mathbf{v}_{i}|}$
(1)

Modelling the time-evolution of $\mathbf{m}$ allows us to gain an
understanding of the overall emergent dynamics of the group. For relatively
small group sizes, idiosyncrasies in the individual $\mathbf{v}_{i}$’s can
have a significant effect on the group dynamics. These result in stochastic
fluctuations which cannot be ignored while modelling the dynamics of
$\mathbf{m}$. Therefore, we use the framework of _stochastic differential
equations (SDE)_ to model the $\mathbf{m}$, which models both deterministic as
well as stochastic aspects of the time-evolution of $\mathbf{m}$.

###  2.3 Data-driven discovery of stochastic differential equations

Given the time series of the 2-dimensional polarization vector $\mathbf{m}(t)$
sampled with a finite time-interval $\Delta t$, our goal is to discover an Itô
stochastic differential equation (SDE) model that explains the time-series, of
the following form:

$\displaystyle\frac{d\mathbf{m}}{dt}=\mathbf{f}(\mathbf{m})+\mathbf{g}(\mathbf{m})\cdot\boldsymbol{\eta}(t)$
(2)

where $\boldsymbol{\eta}$ is 2-dimensional white noise process. Here,
$\mathbf{f}$ is called the _drift_ function and $\mathbf{g}$ is called the
_diffusion_ function. The goal in the SDE discovery problem is to discover
representations for $\mathbf{f}$ and $\mathbf{g}$. In this work, we use a
likelihood-based framework, first introduced by (Dietrich et al., 2021). Here,
$\mathbf{f}$ and $\mathbf{g}$ are represented using neural networks, trained
to maximize a likelihood function based on the finite-time transition
probability.

Specifically, let $p(\cdot,t_{1}|\mathbf{m}_{0},t_{0})$ be the probability
density function of $\mathbf{m}(t_{1})$ conditional to
$\mathbf{m}(t_{0})=\mathbf{m}_{0}$. Then, for a small time-step $\Delta t$,
$p(\cdot,t+\Delta t|\mathbf{m}_{0},t)$ can be approximated based on the
stochastic Euler approximation as:

$\displaystyle p(\cdot,t+\Delta
t|\mathbf{m}_{0},t)\sim\mathcal{N}\left(\mathbf{m}_{0}+\mathbf{f}(\mathbf{m}_{0})\cdot\Delta
t,\mathbf{g}(\mathbf{m}_{0})\cdot\Delta t\right)$ (3)

Parameterizing $p$ in terms of the drift and diffusion functions as
$p_{\mathbf{f},\mathbf{g}}$, this gives rise to the following log-likelihood
loss function which can be maximized to fit $\mathbf{f}$ and $\mathbf{g}$:

$\displaystyle\mathcal{L}\left(\mathbf{f},\mathbf{g}\middle|\mathbf{m}_{0},\mathbf{m}_{1}\right)=\log
p_{\mathbf{f},\mathbf{g}}\left(\mathbf{m}_{1},t+\Delta
t\middle|\mathbf{m}_{0},t\right)$
$\displaystyle=\frac{(\mathbf{m}_{1}-\mathbf{m}_{0}-\mathbf{f}(\mathbf{m})\Delta
t)^{2}}{\mathbf{g}(\mathbf{m}_{0})^{2}\Delta
t}+\log\left|\mathbf{g}(\mathbf{m}_{0})^{2}\Delta t\right|+\log 2\pi$ (4)

where $\mathbf{m}_{0}$ and $\mathbf{m}_{1}$ are successive points in the time
series dataset, sampled $\Delta t$ apart. The loss function $\mathcal{L}$ can
be minimized to fit an appropriate representation of $\mathbf{f}$ and
$\mathbf{g}$.

We represent and $\mathbf{g}$ using neural networks, which were trained to
minimize the loss function in Equation 4–see Section 2.4 for further details
about the implementation. Figure 1 illustrates the full data-driven model
discovery procedure, starting from the individual trajectories, computing the
polarization order parameter, estimating drift and diffusion using neural
networks, and visualizing them as fields.

###  2.4  Implementation details

#### Dataset preparation:

The agent-based schooling model was implemented and simulated in MATLAB, using
the stochastic simulation algorithm (SSA) (Gillespie et al., 2007). From the
trajectories of individual velocity vectors produced by the simulation, the
group polarization $\mathbf{m}(t)$ was computed (Section 2.2), which was used
for training the model. Similarly, for the fish schooling dataset, the tracked
individual trajectories were tracked (see Section 2.1) were used to compute
the group polarization $\mathbf{m}(t)$.

In both the agent-based models as well as real-world datasets, there is
rotational and mirror symmetry in the dynamics (the simulations has no angular
bias, while the fish schools were swimming in a circular arena). Therefore, we
expect these symmetries to carry over to the mean-field models as well. To
encourage the neural network to learn these symmetries, we augment the time-
series of $\mathbf{m}(t)$ with rotated (16 rotations equally spaced between
$-\pi$ and $\pi$) and flipped (horizontal and vertical flips) versions of
itself.

#### Architecture and training details:

The drift and diffusion functions $\mathbf{f}$ and $\mathbf{g}$ were
represented using fully-connected neural networks with 5 layers and 150
neurons per layer, with ELU activation function (Clevert et al., 2015). The
network was trained with a batch size 512 and test-train split of 90:10, with
the Adamax optimizer. To validate the discovered models, we generate simulated
trajectories using the model, and compare the probability density (histogram)
and autocorrelation function of the simulated trajectory to that of the input
trajectory. The similarity between the histograms was quantified using the
Wasserstein metric, $W_{1}$ (Ramdas et al., 2017). The similarity of the
autocorrelation was quantified using the relative timescale discrepancy,
$T_{\text{rel}}=|\tau-\hat{\tau}|/\tau$, where $\tau$ and $\hat{\tau}$ are the
autocorrelation times of the real and simulated time series respectively. A
good fit according to these metrics ensures that the discovered dynamical
model matches the actual dynamics at least in the weak (distributional) sense
(Kloeden et al., 1992).

Figure 2: Goodness-of-fit metrics for mesoscale models for agent-based
simulations. (a, c) Comparison of histograms of $|\mathbf{m}|$ from the actual
time series as well as a time series simulated with the discovered neural
stochastic model; for the pairwise (a) and ternary (c) interaction models. (b,
d) A similar comparison between the autocorrelation functions of the actual
and simulated $|\mathbf{m}|$ time series from the two models.

###  2.5 Visualizing the drift and diffusion fields

Since the drift and diffusion functions are represented as neural networks
(and not equations that are immediately interpretable), a proper visualization
is crucial to enable easy interpretation of the discovered dynamical model. To
do this, we visualize the drift and diffusion functions as fields.

Visualizing $\mathbf{f}$ is relatively straightforward. Since
$\mathbf{f}:\mathbb{R}^{2}\to\mathbb{R}^{2}$ is a vector function, it can be
visualized as a vector field—see Figure 3 (a) for an example. The arrows
indicate the strength and direction of the drift field at each value of
$\mathbf{m}$. The vector field representation makes it easy to identify
features equilibrium points (the vectors will have zero length at the
equilibrium point) and their stability (near an equilibrium point, the vectors
will point towards it if it is stable, and away from it if it is unstable).

Visualizing the diffusion field is more intricate: since
$\mathbf{g}:\mathbb{R}^{2}\to\mathbb{R}^{2\times 2}$ is a matrix function, the
diffusion field is a tensor field. First, notice that for a given
$\mathbf{m}$, the covariance matrix of $\mathbf{m}$ is given by
$G(\mathbf{m})=\mathbf{g}(\mathbf{m})\mathbf{g}(\mathbf{m})^{T}$. For any
value of $\mathbf{m}$, $G(\mathbf{m})$ can be represented by an ellipse
centered at $\mathbf{m}$, with its axes parallel to the eigenvectors of
$G(\mathbf{m})$ and axis lengths proportional to the eigenvalues of
$G(\mathbf{m})$. This is a representation of the diffusion tensor field, and
gives an overview of the strength and directionality of the diffusive force
for each value of $\mathbf{m}$. For readers familiar with tensor field
visualizations, this is a visualization of the $G$-tensor with elliptical
glyphs.

##  3 Results

###  3.1 Mesoscopic equations for simulation models of collective behaviour

We first used our approach to discover the mesoscopic dynamics from simulated
trajectories, for a class of models where the mesoscopic SDEs are already
known—see Section 2.1. This serves as a test-bed for the efficacy of the
approach. With pairwise as well as ternary interactions, the simulated flocks
show a high degree of polarization. However, the mechanism by which order is
created is fundamentally different in these two models.

Figure 2 compares the histogram and autocorrelation of a simulated time series
from the discovered models, to that of the input time series used to train the
model (also see Table 1). There is a close match between the input and the
generated data (_pairwise:_ $W_{1}=0.0686,T_{\text{rel}}=0.2770$, _ternary:_
$W_{1}=0.0270,T_{\text{rel}}=0.0437$). We emphasize that this is despite the
model not being directly trained to match the histogram or the
autocorrelation—suggesting that the discovered models are good approximations
for the actual mesoscopic dynamics.

Figure 3: Mesoscopic dynamics of simulated fish schools. The drift and
diffusion fields estimated by the neural model agrees with the theoretically
predicted drift and diffusion fields. (a) Drift field for a simulated school
with only pairwise interactions, showing a single attractor at $\mathbf{m}=0$.
(b) Diffusion field for the pairwise model simulation, showing that the
strength of diffusion is maximum at $\mathbf{m}=0$ and decreases outwards. (c)
Drift field for a simulated school with higher-order (ternary) interactions,
showing a ring-shaped attractor for a high value of $|\mathbf{m}|$. (d)
Diffusion field for the ternary model simulation, showing the same pattern of
diffusion as in the pairwise model.

The estimated drift and diffusion fields are shown in Figure 3. For a model
with pairwise interactions, the drift field shows a single attractor at the
origin, denoting that the deterministic dynamics has the effect of
depolarising the group. However, the strength of diffusion is maximum at the
origin and decreases as $|\mathbf{m}|$: this means that the stochastic forces
push the system away from the $\mathbf{m}=0$ stable state, creating order in
the system. This counter-intuitive phenomenon of high polarisation arising as
a consequence of stochastic effects is called _noise induced order_ in
literature.

Contrast this with the case of the ternary interaction model, where the the
drift field shows ring-shaped attractor, at a high value of $|\mathbf{m}|$.
The diffusion field is similar to the pairwise model, with diffusion being
strongest near the origin and decreasing outwards. This means that the order
in the ternary interaction model is primarily deterministic.

This demonstrates that our approach is able to distinguish between the
fundamentally different dynamics in the two models considered, despite the
time series and the histogram of the datasets looking qualitatively similar.
Indeed, the estimated models are also in qualitative agreement with
theoretical predictions (Appendix A).

###  3.2  Mesoscopic equations for real-world fish schools

Figure 4: Mesoscopic dynamics of a real-world school of 30 _Etroplus
suratensis_ fish. The observed drift and diffusion fields show evidence for
noise-induced order. (a) Drift field for the _Etroplus school_ , showing a
single attractor at $\mathbf{m}=0$. (b) Diffusion field for school, showing
that the strength of diffusion is maximum at $\mathbf{m}=0$ and decreases
outwards.

We now move on to modelling the stochastic mesoscopic dynamics of a real-world
data set of collective motion, based on a recently published data on schools
of fish (Jhawar et al., 2020). Similar to the simulation models, we seek to
discover mesoscopic SDEs for the dynamics of the polarization order-parameter
$\mathbf{m}$ of real schools of fish.

We observe that the deterministic mean-field dynamics of the school,
characterized by the drift field, drives the system towards a disordered state
($\mathbf{m}=0$). However, the diffusion is maximum at $\mathbf{m}=0$, driving
the system away from the stable state and resulting in noise induced order,
similar to the model system with pairwise interactions. This surprising
phenomenon was first observed in real-world fish schools in (Jhawar et al.,
2020); our neural approach is able to reproduce this result. This is a non-
trivial result: despite the neural networks being able to represent a much
larger hypothesis class, the drift and diffusion fields converged to the
simple forms shown in Figure 4.

We also derived mesoscopic models for fish schools of group sizes 15 and 60.
The results are described in Appendix B and Figure 7, and are qualitatively
similar to the results for the 30-fish school. The net strength of diffusion
goes down as the school size increases, an observation that matches
theoretical predictions (Appendix A).

The fit metrics for the fish schooling dataset (30 fish) is shown in Figure 5
(also see Table 1). Like the agent-based models, there is a good agreement
between the histograms of the actual and the simulated $|\mathbf{m}|$ time
series ($W_{1}=0.0823$). However, autocorrelation function shows a deviation
from the actual time-series ($T_{\text{rel}}=0.6652$). This deviation is
present in other group sizes also (15 and 60 fish). This deviation is because
of a spurious periodicity in the time series of $|\mathbf{m}|$—caused by the
fish swimming along the boundary of the tank—which is not captured by the
mesoscopic SDE model. The deviation is also present with other SDE estimation
techniques (see Table 1), and is a constraint of the model itself, and not the
estimation procedure.

Table 1 summarizes the quantitative performance of the neural effective SDE
estimation technique, on the five different datasets considered in this study
(2 simulated agent-based models, 3 real-world fish schools of different
sizes). For comparison, we use a package PyDaDDy (Nabeel et al., 2022), which
estimates SDEs from data using classical techniques, viz. jump-moment
estimation and equation-learning with sparse-regression. While we reemphasize
that the focus of a study like this is to obtain interpretable models;
overall, the neural SDE estimation approach performs quantitatively better
than PyDaddy, in both the Wasserstein metric $W_{1}$ and relative timescale
discrepancy $T_{\text{rel}}$.

Figure 5: Goodness-of-fit metrics for mesoscale models for a real-world fish
school. (a) Comparison of histograms of $|\mathbf{m}|$ from the actual time
series as well as a time series simulated with the discovered neural
stochastic model; for the schooling dynamics of a school of 30 fish. (b)
Comparison of the autocorrelation functions between the real and simulated
time series for the same dataset.

##  4 Discussion

Compared to conventional techniques for SDE identification, neural networks
have the advantage of being less restrictive in the hypothesis space they
encompass, and hence are more general. It is worth noting that, even with the
added flexibility, the neural approach still discovers qualitatively the same
model for fish schools as (Jhawar et al., 2020), lending further credence to
the observation of noise-induced order in the original study.

Indeed, we expect the neural-network based approach to be more powerful—even
necessary—in scenarios with more complex movement, which lack these
symmetries. For example, organisms could be moving under the influence of an
inhomegenous light or chemical field (Tian et al., 2021; Puckett et al.,
2018). In such cases, enforcing user-defined parametric models, or relying on
simple hypothesis classes (such as polynomials), may impose strong and often
undesirable inductive biases on the SDE discovery process, which the neural-
net-based approach avoids.

Further, the likelihood-based estimation procedure is readily generalizable to
more complex stochastic dynamics. While such extensions are out of the scope
of the current work, it is worth highlighting that the framework is quite
general, and can be extended to accommodate jump discontinuities, extrinsic
driver variables, non-stationarity, etc., all of which are commonly observed
in real-world datasets (Carpenter et al., 2020; Salmaso, 2000).

Finally, we emphasize the need for proper visualization tools to facilitate
interpretability of the discovered models. When the goal of model discovery is
scientific enquiry (and not merely to obtain a predictive model), tools to
interpret and analyze the discovered models become crucial. Our proposed way
of visualizing the drift and diffusion fields achieves exactly this, and
enables one to easily understand the deterministic phase portrait and the
stochasticity landscape of the discovered model.

Table 1: Comparison of model performance with a classical technique for data-
driven SDE discovery
Dataset | Neural SDE | PyDaddy
---|---|---
| $W_{1}$ | $T_{\text{rel}}$ | $W_{1}$ | $T_{\text{rel}}$
ABM (Pairwise) | 0.0686 | 0.2770 | 0.767 | 0.3020
ABM (Ternary) | 0.0270 | 0.0437 | 0.0283 | 0.1583
Etroplus (15 fish) | 0.0823 | 0.6652 | 0.0818 | 0.7787
Etroplus (30 fish) | 0.0357 | 0.5895 | 0.0263 | 0.7144
Etroplus (60 fish) | 0.0155 | 0.2217 | 0.0302 | 0.4409

##  5 Conclusion

In this study, we utilize a physics-informed, grey-box neural network model to
discover mesoscale SDEs for collective movement (Dietrich et al., 2021;
Evangelou et al., 2022). Using this technique, we discover SDEs for both
simulated and real-world datasets of collective movement. With simulated
datasets, this approach is able to identify the difference between different
models, even when the time series and histograms look qualitatively similar.
Furthermore, the SDEs discovered by this approach in each case are are
qualitatively similar to the ones predicted by theory. For a real-world
dataset of schooling fish, the neural-net approach discovers an SDE that
predicts _noise-induced order_ , reproducing the result of a recent study
(Jhawar et al., 2020) that used conventional SDE estimation techniques.

To conclude, we demonstrate the applicability and versatility of neural
network models for studying stochastic dynamics of collective animal movement,
and of biological systems in general. The combination of a physics-informed,
grey-box estimation process along with proper visualization techniques enabled
us to discover models that are quantitatively more accurate than traditional
approaches (Table 1), while still remaining readily interpretable. As high-
quality datasets of complex phenomena become commonplace in biology and
related fields, these techniques—and clear ways to interpret them—become
indispensable.

## Code & Data Availability

The fish schooling dataset used in study is an openly available dataset,
available at (jiteshjhawar & ashwinkk23, 2020).

The neural SDE identification code was a modified version of (Dietrich, 2021).
Our modified version of the code is available at
https://github.com/utkarshp1161/Neu_sde.

## Acknowledgements

We thank Vivek Jadhav for his help with the agent-based model simulations, and
Somaditya Santra for discussions.

## References

  * Alert & Trepat (2020) Alert, R. and Trepat, X.  Physical models of collective cell migration.  _Annual Review of Condensed Matter Physics_ , 11:77–101, 2020.
  * Be’er & Ariel (2019) Be’er, A. and Ariel, G.  A statistical physics view of swarming bacteria.  _Movement ecology_ , 7(1):1–17, 2019.
  * Biancalani et al. (2014) Biancalani, T., Dyson, L., and McKane, A. J.  Noise-induced bistable states and their mean switching time in foraging colonies.  _Physical review letters_ , 112(3):038101, 2014\.
  * Boninsegna et al. (2018) Boninsegna, L., Nüske, F., and Clementi, C.  Sparse learning of stochastic dynamic equations.  _The Journal of Chemical Physics_ , 148(24):241723, June 2018.  ISSN 0021-9606, 1089-7690.  doi: 10.1063/1.5018409.  URL http://arxiv.org/abs/1712.02432.  arXiv: 1712.02432.
  * Brückner et al. (2019) Brückner, D. B., Fink, A., Schreiber, C., Röttgermann, P. J., Rädler, J. O., and Broedersz, C. P.  Stochastic nonlinear dynamics of confined cell migration in two-state systems.  _Nature Physics_ , 15(6):595–601, 2019.
  * Buhl et al. (2006) Buhl, J., Sumpter, D. J., Couzin, I. D., Hale, J. J., Despland, E., Miller, E. R., and Simpson, S. J.  From disorder to order in marching locusts.  _Science_ , 312(5778):1402–1406, 2006.
  * Camazine et al. (2020) Camazine, S., Deneubourg, J.-L., Franks, N. R., Sneyd, J., Theraula, G., and Bonabeau, E.  Self-organization in biological systems.  In _Self-Organization in Biological Systems_. Princeton university press, 2020.
  * Carpenter et al. (2020) Carpenter, S. R., Arani, B. M., Hanson, P. C., Scheffer, M., Stanley, E. H., and Van Nes, E.  Stochastic dynamics of cyanobacteria in long-term high-frequency observations of a eutrophic lake.  _Limnology and Oceanography Letters_ , 5(5):331–336, 2020.
  * Chen et al. (2018a) Chen, R. T., Rubanova, Y., Bettencourt, J., and Duvenaud, D. K.  Neural ordinary differential equations.  _Advances in neural information processing systems_ , 31, 2018a.
  * Chen et al. (2018b) Chen, X., Treiber, M., Kanagaraj, V., and Li, H.  Social force models for pedestrian traffic–state of the art.  _Transport reviews_ , 38(5):625–653, 2018b.
  * Clevert et al. (2015) Clevert, D.-A., Unterthiner, T., and Hochreiter, S.  Fast and accurate deep network learning by exponential linear units (elus).  _arXiv preprint arXiv:1511.07289_ , 2015.
  * Dietrich (2021) Dietrich, F.  SDE Identification.  https://gitlab.com/felix.dietrich/sde-identification, 2021.
  * Dietrich et al. (2021) Dietrich, F., Makeev, A., Kevrekidis, G., Evangelou, N., Bertalan, T., Reich, S., and Kevrekidis, I. G.  Learning effective stochastic differential equations from microscopic simulations: combining stochastic numerics and deep learning.  _arXiv:2106.09004 [physics]_ , 2021.  URL http://arxiv.org/abs/2106.09004.  arXiv: 2106.09004.
  * Dinet et al. (2021) Dinet, C., Michelot, A., Herrou, J., and Mignot, T.  Linking single-cell decisions to collective behaviours in social bacteria.  _Philosophical Transactions of the Royal Society B_ , 376(1820):20190755, 2021.
  * Dyson et al. (2015) Dyson, L., Yates, C. A., Buhl, J., and McKane, A. J.  Onset of collective motion in locusts is captured by a minimal model.  _Physical Review E_ , 92(5):052708, 2015.
  * Evangelou et al. (2022) Evangelou, N., Dietrich, F., Bello-Rivas, J. M., Yeh, A., Stein, R., Bevan, M. A., and Kevekidis, I. G.  Learning Effective SDEs from Brownian Dynamics Simulations of Colloidal Particles.  _arXiv:2205.00286 [cs, math]_ , April 2022.  URL http://arxiv.org/abs/2205.00286.  arXiv: 2205.00286.
  * Gardiner (2009) Gardiner, C.  _Stochastic methods_ , volume 4.  springer Berlin, 2009.
  * Gillespie et al. (2007) Gillespie, D. T. et al.  Stochastic simulation of chemical kinetics.  _Annual review of physical chemistry_ , 58(1):35–55, 2007.
  * Gorjão & Meirinhos (2019) Gorjão, L. R. and Meirinhos, F.  kramersmoyal: Kramers–moyal coefficients for stochastic processes.  _Journal of Open Source Software_ , 4(44):1693, 2019.
  * Hasan et al. (2021) Hasan, A., Pereira, J. M., Farsiu, S., and Tarokh, V.  Identifying latent stochastic differential equations.  _IEEE Transactions on Signal Processing_ , 70:89–104, 2021\.
  * Jhawar & Guttal (2020) Jhawar, J. and Guttal, V.  Noise-induced effects in collective dynamics and inferring local interactions from data.  _Philosophical Transactions of the Royal Society B_ , 375(1807):20190381, 2020.
  * Jhawar et al. (2019) Jhawar, J., Morris, R. G., and Guttal, V.  Deriving Mesoscopic Models of Collective Behavior for Finite Populations.  In _Handbook of Statistics_ , volume 40, pp. 551–594. Elsevier, 2019.  ISBN 978-0-444-64152-6.  doi: 10.1016/bs.host.2018.10.002.
  * Jhawar et al. (2020) Jhawar, J., Morris, R. G., Amith-Kumar, U., Danny Raj, M., Rogers, T., Rajendran, H., and Guttal, V.  Noise-induced schooling of fish.  _Nature Physics_ , 16(4):488–493, 2020.
  * jiteshjhawar & ashwinkk23 (2020) jiteshjhawar and ashwinkk23.  tee-lab/schooling_fish: More updates on the codes and much elaborated readme file, January 2020.  URL https://doi.org/10.5281/zenodo.3632470.
  * Kloeden et al. (1992) Kloeden, P. E., Platen, E., Kloeden, P. E., and Platen, E.  _Stochastic differential equations_.  Springer, 1992.
  * Li et al. (2020) Li, X., Wong, T.-K. L., Chen, R. T., and Duvenaud, D.  Scalable gradients for stochastic differential equations.  In _International Conference on Artificial Intelligence and Statistics_ , pp. 3870–3882. PMLR, 2020.
  * Liu et al. (2020) Liu, J., Long, Z., Wang, R., Sun, J., and Dong, B.  Rode-net: learning ordinary differential equations with randomness from data.  _arXiv preprint arXiv:2006.02377_ , 2020.
  * McKane & Newman (2004) McKane, A. J. and Newman, T. J.  Stochastic models in population biology and their deterministic analogs.  _Physical Review E_ , 70(4):041902, 2004.
  * Nabeel et al. (2022) Nabeel, A., Karichannavar, A., Palathingal, S., Jhawar, J., Danny Raj, M., and Guttal, V.  Pydaddy: A python package for discovering stochastic dynamical equations from timeseries data.  _arXiv preprint arXiv:2205.02645_ , 2022.
  * Nielsen et al. (2000) Nielsen, J. N., Madsen, H., and Young, P. C.  Parameter estimation in stochastic differential equations: an overview.  _Annual Reviews in Control_ , 24:83–94, 2000.
  * O’Leary et al. (2021) O’Leary, J., Paulson, J. A., and Mesbah, A.  Stochastic physics-informed neural networks (spinn): A moment-matching framework for learning hidden physics within stochastic differential equations.  _arXiv preprint arXiv:2109.01621_ , 2021.
  * Puckett et al. (2018) Puckett, J. G., Pokhrel, A. R., and Giannini, J. A.  Collective gradient sensing in fish schools.  _Scientific reports_ , 8(1):1–11, 2018.
  * Ramaswamy (2010) Ramaswamy, S.  The mechanics and statistics of active matter.  _Annual Review of Condensed Matter Physics_ , 1(1):323–345, 2010.  doi: 10.1146/annurev-conmatphys-070909-104101.  URL https://doi.org/10.1146/annurev-conmatphys-070909-104101.
  * Ramaswamy (2017) Ramaswamy, S.  Active matter.  _Journal of Statistical Mechanics: Theory and Experiment_ , 2017(5):054002, 2017.
  * Ramdas et al. (2017) Ramdas, A., García Trillos, N., and Cuturi, M.  On wasserstein two-sample testing and related families of nonparametric tests.  _Entropy_ , 19(2):47, 2017.
  * Rørth (2009) Rørth, P.  Collective cell migration.  _Annual Review of Cell and Developmental Biology_ , 25:407–429, 2009.
  * Salmaso (2000) Salmaso, N.  Factors affecting the seasonality and distribution of cyanobacteria and chlorophytes: a case study from the large lakes south of the alps, with special reference to lake garda.  _Hydrobiologia_ , 438(1):43–63, 2000.
  * Song et al. (2020) Song, Y., Sohl-Dickstein, J., Kingma, D. P., Kumar, A., Ermon, S., and Poole, B.  Score-based generative modeling through stochastic differential equations.  _arXiv preprint arXiv:2011.13456_ , 2020.
  * Sumpter (2010) Sumpter, D. J.  Collective animal behavior.  In _Collective Animal Behavior_. Princeton University Press, 2010\.
  * Tabar (2019) Tabar, R.  _Analysis and data-based reconstruction of complex nonlinear dynamical systems_ , volume 730.  Springer, 2019.
  * Tian et al. (2021) Tian, M., Zhang, C., Zhang, R., and Yuan, J.  Collective motion enhances chemotaxis in a two-dimensional bacterial swarm.  _Biophysical journal_ , 120(9):1615–1624, 2021\.
  * Van Kampen (1992) Van Kampen, N. G.  _Stochastic processes in physics and chemistry_ , volume 1.  Elsevier, 1992.
  * Vicsek & Zafeiris (2012) Vicsek, T. and Zafeiris, A.  Collective motion.  _Physics reports_ , 517(3-4):71–140, 2012.
  * Vicsek et al. (1995) Vicsek, T., Czirók, A., Ben-Jacob, E., Cohen, I., and Shochet, O.  Novel Type of Phase Transition in a System of Self-Driven Particles.  _Physical Review Letters_ , 75(6):1226–1229, August 1995.  ISSN 0031-9007, 1079-7114.  doi: 10.1103/PhysRevLett.75.1226.  URL https://link.aps.org/doi/10.1103/PhysRevLett.75.1226.
  * Wang et al. (2022) Wang, Y., Fang, H., Jin, J., Ma, G., He, X., Dai, X., Yue, Z., Cheng, C., Zhang, H.-T., Pu, D., et al.  Data-driven discovery of stochastic differential equations.  _Engineering_ , 2022.
  * Yates et al. (2009) Yates, C. A., Erban, R., Escudero, C., Couzin, I. D., Buhl, J., Kevrekidis, I. G., Maini, P. K., and Sumpter, D. J. T.  Inherent noise can facilitate coherence in collective swarm motion.  _Proceedings of the National Academy of Sciences_ , 106(14):5464–5469, April 2009.  ISSN 0027-8424, 1091-6490.  doi: 10.1073/pnas.0811195106.  URL https://pnas.org/doi/full/10.1073/pnas.0811195106.

##  Appendix A  Analytically derived mesoscopic SDEs for agent-based models

Figure 6: Analytically derived mesoscale dynamics of agent based models. (a)
Analytically derived drift field for a model with only pairwise interactions,
showing a single attractor at $\mathbf{m}=0$. (b) Analytically derived
diffusion field for the pairwise model, showing that the strength of diffusion
is maximum at $\mathbf{m}=0$ and decreases outwards. (c) Analytically derived
drift field for a model with higher-order (ternary) interactions, showing a
ring-shaped attractor for a high value of $|\mathbf{m}|$. (d) Analytically
derived diffusion field for the ternary model simulation, showing the same
pattern of diffusion as in the pairwise model.

The agent-based models of collective behaviour used in this study have been
well-studied, and mesoscopic SDEs for these models have been analytically
derived (Jhawar et al., 2019, 2020).

Recap the model description:

  * •

An agent $i$ can spontaneously turn and choose a new direction, i.e.
$\theta_{i}(t)\rightarrow\eta$ for $\eta\sim\mathsf{Unif}[-\pi,\pi)$. The
spontaneous turns happen as a Poisson process with a rate $r_{1}$.

  * •

An agent $i$ can copy the direction of one (pairwise interaction) or the
average direction of two (ternary interaction) other agents, chosen randomly
from the group. The copy interactions can happen at a rates $r_{2}$ and
$r_{3}$ respectively.

For the pairwise interaction model, the mesoscale SDE can be analytically
derived to be of the following form:

$\displaystyle\frac{d\mathbf{m}}{dt}$
$\displaystyle=-r_{1}\mathbf{m}+\sqrt{\frac{r_{1}+r_{2}\left(1-|\mathbf{m}|^{2}\right)}{N}}I\cdot\boldsymbol{\eta}(t)$
(5)

where $N$ is the total number of agents in the group.

Similarly, for the ternary interaction model, the derived mesoscale SDE has
the form:

$\displaystyle\frac{d\mathbf{m}}{dt}$
$\displaystyle=-r_{1}\mathbf{m}+r_{3}\left(1-|\mathbf{m}|^{2}\right)\mathbf{m}+\sqrt{\frac{r_{1}+(r_{2}+r_{3})\left(1-|\mathbf{m}|^{2}\right)}{N}}I\cdot\boldsymbol{\eta}(t)$
(6)

The diffusion term is a diagonal matrix, and has a $1-|\mathbf{m}|^{2}$ term.
This means that the noise is maximum at $|\mathbf{m}|=0$ and decreases with
increasing $|\mathbf{m}|$. The drift term is linear in $\mathbf{m}$ for the
pairwise interaction model, and the deterministic stable equilibrium for this
system is at $|\mathbf{m}|=0$. Therefore, the order in the pairwise
interaction model is noise-induced, and arises solely due to the high strength
of diffusion near $|\mathbf{m}|=0$. On the other hand, for the ternary
interaction model, drift term is cubic in $\mathbf{m}$, and has a
$\left(1-|\mathbf{m}|^{2}\right)\mathbf{m}$ term. This contributes to a ring-
shaped attractor near $|\mathbf{m}|=1$, i.e. the order in the ternary
interaction model originates from a combination of both deterministic and
noise-induced effects.

The drift and diffusion fields generated from equations 5 and 6 are shown in
Figure 6. There is a close correspondence between these plots and the ones in
Figure 3, suggesting that the neural effective SDE approach is able to
correctly infer the underlying mesoscale SDEs.

##  Appendix B  Mesoscopic models for fish schools of different group sizes

Figure 7: Estimated drift and diffusion fields for polarization dynamics of
fish schools of different group sizes. (a, b) Drift and diffusion fields for
$\mathbf{m}$ for a 15-fish school (c, d) Drift and diffusion fields for
$\mathbf{m}$ for a 15-fish school, same as in Figure 4. (e, f) Drift and
diffusion fields for $\mathbf{m}$ for a 60-fish school.

We repeated the analysis in Section 3.2 on two other group sizes, _viz._
$N=15$ and $N=60$. The fit metrics and drift/diffusion fields are shown in
Figure. 7. Overall, the fields look qualitatively similar across different
group sizes. However, the net strength of diffusion decreases with increasing
group size (notice the range of the color axis). This is in accordance with
theoretical predictions: in theory, we expect the diffusion term $\mathbf{g}$
to scale as $1/\sqrt{N}$ (see Equation 5).


# A goodness-of-fit test based on a recursive product of spacings

Philipp Eller 11footnotetext<EMAIL_ADDRESS>and Lolian Shtembari
22footnotetext<EMAIL_ADDRESS>
###### Abstract

We introduce a new statistical test based on the observed spacings of ordered
data. The statistic is sensitive to detect non-uniformity in random samples,
or short-lived features in event time series. Under some conditions, this new
test can outperform existing ones, such as the well known Kolmogorov-Smirnov
or Anderson-Darling tests, in particular when the number of samples is small
and differences occur over a small quantile of the null hypothesis
distribution. A detailed description of the test statistic is provided
including a detailed discussion of the parameterization of its distribution
via asymptotic bootstrapping as well as a novel per-quantile error estimation
of the empirical distribution. Two example applications are provided, using
the test to boost the sensitivity in generic "bump hunting", and employing the
test to detect supernovae. The article is rounded off with an extended
performance comparison to other, established goodness-of-fit tests.

##  1 Introduction

Assessing the goodness-of-fit of a distribution given a number of random
samples is an often-encountered problem in data analysis. Such statistical
hypothesis tests find applications in many fields, ranging from the natural
and social sciences over engineering to quality control. Several non-
parametric tests exist, some of which have become standard tools, including
the Kolgogorov-Smirnov (KS) test [1, 2] or the Anderson-Darling (AD) test [3].
[4] provides a comprehensive overview of existing tests, and a comparison of
their performance for the case of detecting non-uniformity for a set of
alternative distributions.

In this work, we are in contrast interested in the case where the bulk of
samples are actually distributed according to the null hypothesis, and only
few additional samples are introduced that are following a different
distribution, representing a narrow excess over a known background. We present
the new test statistic "recursive product of spacings", or short RPS, that is
based on the spacings between ordered samples, and introduced in Sec. 2. In
Sec. 3 we provide a parametrization of its distribution based on simulations,
introducing techniques to estimate the asymptotic result of infinite
bootstrapping steps in order to improve the quality of our fits. Subsequently
we discuss the quality of the approximation deriving a per-quantile error
estimate up to a desired confidence level. The rest of the article focuses on
some illustrations and example applications, as well as a detailed performance
comparison to several other test statistics.

###  1.1 Goodness-of-fit Tests

Suppose that we have obtained $n$ samples $x_{i}$, and want to quantitatively
test the hypothesis of those samples being random variates of a known
distribution $f(x)$, i.e. independent and identically distributed (i.i.d.)
according to $f(x)$. Here, we consider only continuous distributions $f(x)$
with cumulative $F(x)$, and hence can transform samples onto the unit interval
$[0,1]$ via $y_{i}=F(x_{i})$. This reduces the task at hand to test
transformed samples $y_{i}$ being distributed according to the standard
uniform distribution $\mathcal{U}(0,1)$.

First, let us briefly introduce other, existing test statistics to which we
will compare the RPS statistic. We consider in particular two groups of
statistics, those based on the empirical distribution (EDF Statistics), and
those based on the spacings between ordered samples (Spacings Statistics). An
comprehensive overview of existing test statistics can be found in [4].

####  1.1.1 EDF Statistics

This class of test statistics compares the empirical distribution function
(EDF) $F_{n}(x)$ to the cumulative distribution function (CDF) $F(x)$, (here
$F(x)=x$). In particular the following are widely used and considered here:

  * •

Kolmogorov-Smirnov (KS) [1, 2]: $D_{n}=\sup_{x}|F_{n}(x)-F(x)|$

  * •

Cramer-von-Mises (CvM) [5, 6]:
$T=n\int_{-\infty}^{\infty}(F_{n}(x)-F(x))^{2}dF(x)$

  * •

Anderson-Darling (AD) [3]:
$A^{2}=n\int_{-\infty}^{\infty}\frac{(F_{n}(x)-F(x))^{2}}{F(x)\;(1-F(x))}\,dF(x)$

Similar are a type of statistics defined on the ordered set. Given the $n$
samples $\\{x_{1},x_{2},\dots,x_{n}\\}$, we define the ordered set of samples
as $\\{x_{(1)},x_{(2)},\dots,x_{(n)}\\}$, where $x_{(i)}<x_{(i+1)}\,\forall
i$. The expected value of ordered sample $i$ is $i/(n+1)$, and we define the
deviation to the expected values as $\delta_{i}=x_{(i)}-i/(n+1)$ for each
sample $i$. Based on this we can write out the following two statistics:

  * •

Pyke’s Modified KS (C) [7, 8]:
$C_{n}=\max(\max(\delta_{i}),-\min(\delta_{i}))$

  * •

Brunk’s Modified KS (K) [9]: $K_{n}=\max(\delta_{i})-\min(\delta_{i})$

####  1.1.2 Spacings Statistics

Based on the ordered set, we can further define the $n+1$ spacings $s$ as
$s_{i}=x_{(i)}-x_{(i-1)}$, with $x_{(0)}=0$ and $x_{(n+1)}=1$. Several test
statistics built from these spacings are considered in literature, including:

  * •

Moran (M) [10]: $M=-\sum_{i=1}^{n+1}\log s_{i}$

  * •

Greenwood (G) [11]: $G=\sum_{i=1}^{n+1}s^{2}_{i}$

In the context of a fixed rate Poisson process, these spacings can also be
interpreted as interarrival times or waiting times. In some other areas,
spacings are also referred to as gaps.

So-called higher order spacings can be defined by summing up neighbouring
spacings. Here we consider the overlapping $m$-th order spacings
$s_{i}^{(m)}=x_{(i+m)}-x_{(i)}$. With those, we can define generalisations of
Moran and Greenwood, respectively, as discussed by Cressie:

  * •

Logarithms of higher order spacings (Lm) [12]:
$L_{n}^{(m)}=-\sum_{i=0}^{n-m+1}\log s_{i}^{(m)}$

  * •

Squares of higher order spacings (Sm) [13]:
$S_{n}^{(m)}=\sum_{i=0}^{n-m+1}(s_{i}^{(m)})^{2}$

For our comparisons presented later, we choose $m=2$ and $m=3$, respectively,
to limit ourselves to a finite list of tests.

Other statistics based on spacings exist and are being actively developed and
used, such as, for example, tests based on the $k$ smallest or largest
spacings [14].

##  2 Recursive Product of Spacings (RPS)

In this work, our goal is to construct a new test statistic, that is sensitive
to narrow features or clusters in an otherwise uniform distribution of
samples. The tell-tale sign we are looking for is a localized group of
uncommonly small spacings of the ordered data. For this purpose, we propose a
new class of test statistics, that are including higher order spacings in a
recursive way.

The recursive product of spacings (RPS) can be thought of as an extension of
the Moran statistic, and is defined as

$RPS(n)=M^{n+1}+M^{n}+\dots+M^{1},$ (2.1)

where the term $M^{n+1}$ is the simple sum of negative log spacings equivalent
to the Moran statistic

$M^{n+1}=-\sum_{i=1}^{n+1}\log\left(s_{i}\right).$ (2.2)

All following terms are computed in the same way

$M^{j}=-\sum_{i=1}^{j}\log\left(s_{i}^{j}\right),$ (2.3)

but with modified spacings $s_{i}^{j}$, defined as:

$s^{j}_{i}=\frac{s^{j+1}_{i}+s^{j+1}_{i+1}}{\sum_{i}s^{j}_{i}}$ (2.4)

which there are $j$ of, and that depend on the spacings $s^{j+1}$ used to
compute the previous term $M^{j+1}$ (hence the recursiveness). In order to
better understand Eq. 2.4 we can turn to Fig. 1, where we show how to
transition from $s^{j+1}$ (top) to $s^{j}$ (bottom): in the top plot we show a
list of events (blue), where we also highlight the boundaries 0 and 1 since
they contribute to defining spacings; in the middle plot the middle points of
the top row spacings are shown, forming a reduced set of "events", which is
then transformed in order to ensure that the spacings of the new set sum up to
1, as shown in the bottom plot; the number of spacings going from the top plot
to the bottom one is reduced by one, showing how we have a finite number of
reduction steps in the definition of the $RPS$.

Figure 1: Example of the reduction step included in the $RPS$ calculation.
Given an initial set of events (top; blue), the middle points are calculated
(top and middle; green), which are then scaled in order to fill the $[0,1]$
interval, forming a new set of data (bottom; red). The evolution of sample
positions on the $[0,1]$ interval are annotated via the arrows.

We can see that term $M^{n}$ is identical to $L_{n}^{(2)}$ up to a
normalization factor $1/{\sum_{i}s_{i}}$. It is important to include such a
normalization for the spacings of each level, as this ensures that the case of
equidistant spacings—the most regular and uniform case—yields the smallest
possible RPS value. This minimum value of $RPS(n)$, given by the configuration
of equidistant samples, can be expressed easily, as each spacing $s^{j}$ is
equal to $\frac{1}{j}$, and thus:

$\min(RPS(n))=-\sum_{j=1}^{n+1}\sum_{i=1}^{j}\log\left(\frac{1}{j}\right)=\sum_{j=i}^{n+1}j\cdot\log(j).$
(2.5)

At the other extreme, very small spacings will yield a large contribution to
the sum of Eq. 2.3, thus $\max(RPS(n))=\infty$ for any given number of samples
$n$. These extrema show that RPS measures the irregularity in sample
positions. The RPS statistic increases the more samples aggregate into local
clusters.

The RPS quantity calculated so far has an infinite support. Since
approximating the tails of distributions is often easier when dealing with a
bounded quantity instead, we transform the $RPS$ in order to bound its support
in the range $[0,1]$, via

$RPS^{*}(n)=\frac{\min(RPS(n))}{RPS(n)}$ (2.6)

into a new quantity $RPS^{*}$ that we will consider when using our test-
statistic.

The following pseudo code illustrates how the computation of the RPS value can
be implemented:

Algorithm 1 Calculates the recursive product of spacings $rps$ from ordered
samples $x_{(i)}$

$x=[0,x_{(1)},x_{(2)},\dots,x_{(n)},1]$

$rps=0$

$s=x[\texttt{first}+1:\texttt{last}]-x[\texttt{first}:\texttt{last}-1]$
$\triangleright$ initial spacings

while $len(s)>1$ do

$rps=rps-sum(log(s))$

$s=s[\texttt{first}:\texttt{last}-1]+s[\texttt{first}+1:\texttt{last}]$
$\triangleright$ spacings for next iteration

$s=s/sum(s)$ $\triangleright$ normalize

$normalized\\_rps=min\\_rps(n)/rps$

This algorithm has a computational complexity of $\mathcal{O}(n^{2})$, and can
become inefficient for very large sample sizes $n$. In this work we limit
ourselves to $n\leq 1000$.

In an analogue way, we can also define an extension to Greenwood $G(n)$, that
instead of logarithms of spacings, sums over the squares of spacings. This
means that we substitute Eq. 2.3 with
$M^{j}=\sum_{i=1}^{j}\left(s_{i}^{j}\right)^{2}$, while keeping the definition
of $s_{i}^{j}$ from Eq. 2.4. We call this recursive form the "RSS" test
statistic in the following comparison.

###  2.1 Illustration

To illustrate better how our test statistic works, and to highlight
differences to other tests, we use the example set of samples drawn from a
uniform (null hypothesis $H^{0}$) and a non-uniform distribution,
respectively, shown in Fig. 2.

Figure 2: Example of 15 standard uniformly distributed samples (left) and 10
standard uniformly + 5 normally ($\mu=0.5,\sigma=0.1$) distributed samples
(right). The sample positions on the $[0,1]$ interval are annotated via the
arrows + text.

The Moran test is based on the spacings between samples, and the smallest and
largest spacings in the specific example are present in the uniform case. This
leads to a more extreme test statistic value $t$ and hence p-value $p=P(T\geq
t|H^{0})$ of 0.117 for the uniform case, while it evaluates to $p=0.335$ in
the non-uniform case. The feature of samples clustering locally together, as
in the non-uniform case, is—by construction—completely missed by Moran’s test.

The KS test can detect such clustering via the CDF, however in our chosen
example it is challenged by the fact that samples trend towards the left in
the uniform case, while they are more balanced in the non-uniform case. This
leads to p-values of 0.048 for uniform, and 0.356 for non-uniform,
respectively.

The RPS test, however, taking into account also spacing between spacings,
finds a p-value of 0.532 for the uniform case, and a much lower p-value of
0.057 for the non-uniform samples. The behaviour of RPS is further illustrated
in Fig. 3, that shows the individual contribution of spacings of all recursion
levels that build up the test statistic value. The Moran statistic corresponds
to the sum over the first row ($M^{16}$), while all subsequent levels are
added for RPS. The uniform samples exhibit the most extreme values in the
first level, but then even out rapidly. Contrary to that, the non-uniform
samples with their clustering give rise to larger values at later levels,
explaining the lower observed p-value.

Figure 3: Illustration of the test statistic contributions from all recursion
levels for the uniformly distributed samples (left) and the non-uniform
samples (right). The sum over the first level only ($M^{16}$) is equivalent to
the Moran statistic.

##  3 Cumulative Distribution of RPS

In order to use RPS as a statistical test yielding p-values, we need its
cumulative distribution $F$. In the case of $n=1$ that has only two
spacings—the simplest non-trivial case we can encounter—we can easily derive
the distribution of $RPS^{*}(1)$, which is:

$F_{RPS^{*}}(x;n=1)=1-\sqrt{1-4^{\frac{x-1}{x}}}$ (3.1)

For $n\geq 2$, however, it is not simple to derive this distribution.
Therefore, we resort to numerically approximating the distribution of
$RPS^{*}$ discussed in the following section.

###  3.1 Approximate Distribution

We have built an approximation for the cumulative distribution
$F_{RPS^{*}}(x;n)$ precise enough to compute meaningful p-values up to
relatively extreme values of up to $10^{-7}$, and large sample sizes $n$ of up
to 1000. Figure 4 shows some examples of $RPS^{*}$ distributions for a few
values of $n$.

Figure 4: Example of CDFs of the $RPS^{*}$ distribution for a few different
values of $n$. N.B.: the x-axis is displayed in inverted logarithm.

We base our approximation on simulation, drawing events with uniform
distribution in the range $[0,1]$ for a given $n$, and collecting $N=2\cdot
10^{8}$ samples of $RPS^{*}(n)$. Such simulation could be directly used to
calculate p-value estimates by counting the fraction of trials below or above
an observed $RPS^{*}$ value $x$ for a fixed $n$. However, we want to provide a
continuous and smooth function valid for any $n\leq 1000$. For this, we use
simulated data to infer the values $x$ of our test statistics corresponding to
a discrete list of specific quantiles $p\in[10^{-7},1-10^{-7}]$. Taking the
$i$-th element in the sorted simulation set gives an estimate for the value of
$x(p=i/N)$. In order to improve this estimate, we could use bootstrapping
[15], collecting different realisations of $x$ by resampling the original
dataset with replacement, resulting in a distribution of values of $x$ for
each $p$, from which we can then extract the mean and the standard deviation,
indicative of the error (see Fig 5). Instead of manually performing the
bootstrapping, we can calculate the probability of each sample $x$ to
represent a specific quantile $p$ if we were to sample randomly with
replacement. For simplicity, let us consider rational quantiles that can be
expressed in the form $p=\frac{k}{N}$; the probability that the $i$-th sample
could end up representing the $k$-th quantile is:

$\pi_{k,i}=F_{B}\left(k,N+1-k;\frac{i}{N}\right)-F_{B}\left(k,N+1-k;\frac{i-1}{N}\right)$
(3.2)

where $F_{B}(a,b;t)$ is the cumulative function of the Beta distribution with
parameters [$a$, $b$] estimated at $t$. The distribution
$\mathrm{Beta}(k,N+1-k)$ represents the $k$-th order statistic of the uniform
distribution [16], i.e. the $k$-th largest element of a set on $N$ uniformly
distributed random variable. Eq. 3.2 corresponds to the limiting case of
performing an infinite number of bootstrapping steps and can be used to
quickly estimate the mean and standard deviation of all $x(p)$ for a choice on
$n$, especially when dealing with large datasets:

$\mathrm{E}\left[x\left(\frac{k}{N}\right)\right]=\sum_{i=1}^{N}x_{i}\cdot\pi_{k,i}$
(3.3)

$\mathrm{Std}\left[x\left(\frac{k}{N}\right)\right]=\sqrt{\sum_{i=1}^{N}\left(x_{i}-\mathrm{E}\left[x\left(\frac{k}{N}\right)\right]\right)^{2}\cdot\pi_{k,i}}$
(3.4)

It would be inefficient to produce such simulation for any $n$, and hence we
repeat the above procedure for only 180 different choices of $n$ between 2 and
1000 following approximately a logarithmic spacing.

###  3.2 Fitting procedure

Using Eq. 3.3 and Eq. 3.4 we are able to define a grid of points with mean
$\mu(n,p)$ and standard deviation $\sigma(n,p)$. Our goal is to estimate a set
of points $\hat{x}(n,p)$, which will be the basis to interpolate and infer the
distribution of the test statistic for all values of $n$ and $p$ defined
above. The points $\hat{x}(n,p)$ is allowed to deviate from the means
$\mu(n,p)$ within the uncertainties $\sigma(n,p)$, and can thereby provide a
more accurate approximation by smoothing out stochastic noise. Additionally,
points from the analytic solution for $n=1$ (Eq. 3.1) are added to the list as
anchor points at the boundary.

Given a trial set $\tilde{x}(n,p)$, we interpolate a cubic spline polynomial
across the values of $n$ for each value of $p$, similarly to the fits shown in
Fig 5. Given one such cubic spline, we evaluate the third derivative on both
sides of each node, calculating the square of their difference and summing up
across all nodes. Since we are using cubic splines, the third derivative is
not continuous, and the "size" of the discontinuity is indicative of the
smoothness of the interpolation. Summing up the contributions form all nodes
of all cubic splines construct the smoothing cost function. The construction
of this cost function is based on [17, 18, 19], where smoothness is treated
very similarly. The estimation of the cubic spline coefficients and the
evaluation of the smoothness cost function can be represented as a quadratic
objective function, which we want to minimize:

$G(\tilde{x})\propto\frac{1}{2}\tilde{x}^{T}\cdot
Q\cdot\tilde{x}+\bar{h}^{T}\cdot\tilde{x}$ (3.5)

In addition to obtaining a smooth fit, there are also some additional
constraints that need to be considered: monotonicity and sum of squared
residuals.

Since the samples $\tilde{x}(p|n)$ should represent a cumulative density
function, then it is important they are properly ordered, ensuring that
$\tilde{x}(p_{i}|n)\leq\tilde{x}(p_{j}|n)$ for $i\leq j$. This is ensured
including a number of linear inequality constraints modelled as a linear
constraint matrix:

$A\cdot\tilde{x}\leq b$ (3.6)

Lastly, we assume that the values $\tilde{x}(n,p)$ are normally distributed
with means $\mu(n,p)$ and standard deviations $\sigma(n,p)$. Since we want to
move away from the initial values $\mu(n,p)$ in order to obtain a smoother
fit, it is important to limit this movement the further away we get and we do
so by considering the sum of squared residuals, which is a typical measure to
account for the global deviation from the mean. Since we assume gaussian
deviations, the sum of all squared residuals can be modelled by a $\chi^{2}$
distribution with $m$ degrees of freedoms, where $m$ is the total number of
parameters, i.e. the number of nodes. Given this distribution, we can estimate
the value of the cost function to be limited to the mean ($m$) plus one
standard deviation ($\sqrt{2m}$) of the $\chi^{2}$ distribution, thus:

$\sum_{i=1}^{m}\frac{\left(\tilde{x}_{i}-\mu_{i}\right)^{2}}{\sigma_{i}^{2}}\leq
m+\sqrt{2\cdot m}$ (3.7)

Figure 5 shows a fitted spline representation of $\hat{x}(n|p)$ for different
values of $p$. Based on the resulting list of corresponding $p$ and $\hat{x}$
values, that we obtained for any $n$, we generate another spline interpolation
as the approximation of the desired cumulative distribution $F(\hat{x};n)$ for
a given $n$. As the cumulative distribution function $F$ is strictly
monotonous in $\hat{x}$, we use the [20] monotonic spline interpolation on the
points $[\hat{x}(p|n),p]$ to produce the final CDFs, shown in Figure 4 for a
few values of $n$.

Figure 5: Example of spline fitted $x$-values across $n$ for a few extreme
p-values. The colored bands show the 1, 2 and 3 sigma bands estimated via
bootstrapping, the black, dashed lines show the approximations by the spline
fits.

###  3.3 Error estimation

Finally, we are also able to estimate the precision of our approximation.
Given any set of i.i.d. random variables, such as $x$, the corresponding list
of estimated quantiles $p$ represents a random set of uniform variates. For
any rational quantile $p_{test}=\frac{k}{N}$ we can estimate the 98% credible
interval $\left(p_{0.01},p_{0.99}\right)$ using the distribution of the $k$-th
order statistic $\mathrm{Beta}(k,N+1-k)$. Given the credible interval, we
calculate the relative error of $p_{test}$ against the extrema of the
interval, considering the largest value representative of the relative error
of a random EDF up to a specified credible level. The results of the estimated
relative error for our choice of $N=2\cdot 10^{8}$ and for quantiles as low as
$p=10^{-7}$ are shown in Figure 6(a).

(a) Estimated relative error of empirical p-value with respect to the 98 %
credible interval and $2\cdot 10^{8}$ samples. The vertical axis reports the
scale of the relative error in percent for two extremes, the 1% and the 99%
quantile of the order statistic distribution.

(b) Estimated relative error of fitted p-value with respect to p-values
obtained via bootstrapping. The vertical axis reports the scale of the
relative error in percent for two extremes, the 1% and the 99% quantile of the
bootstrapping distribution. Results for $n=75$.

Figure 6: Per-quantile relative error stimation of the approximate RPS
distribution.

As expected, the errors are increasing towards smaller p-values and exhibit an
approximately linear behaviour in the log-log plot. We see that the estimated
upper bound of the relative error for a p-value of $10^{-3}$ is below $1\%$,
while for a p-value of $10^{-5}$ it increases to $<10\%$ and ultimately to
$<100\%$ for p-values of $10^{-7}$. Such a "large" relative error for small
p-values may sound alarming at first, but estimating a p-value of $10^{-7}$
and knowing it could actually be closer to $2\cdot 10^{-7}$ would hardly
change the statistical interpretation of a result.

In order to show the validity of these results, we compute the relative error
of our approximate distributions against a test dataset containing 10 times
more samples using bootstrapping. We do so for a few choices of number of
events $n$, and in Figure 6(b) it can bee seen that the behavior of the
relative error is in complete agreement with our analytic estimates of Figure
6(a).

So defined, the relative error $\delta(p|N)$ is a function of the quantile $p$
and number of samples $N$, but this relationship can also be inverted in order
to determine the number of samples necessary to achieve a desired relative
error for a specific quantile: $N(p|\delta)$. Our choice of $N=2\cdot 10^{8}$
was in fact guided by the requirement of having a relative error lower than
100% for a p-value of $10^{-7}$ in at least 99% of cases.

It is worth stressing that these estimates of the relative error are accurate
with respect to the EDF that was sampled for each independent $n$, but might
be subject to small changes after the smoothing fit we performed in order to
regularize and infer the distributions for all missing values of $n$.

###  3.4 Implementation

The RPS test is made available as open-source packages for
Python333https://pypi.org/project/spacings/ and
Julia444https://juliapackages.com/p/spacingstatistics, respectively, with the
p-value parametrizations initially available up to 1000 samples.

Below we give a minimal example to evaluate the RPS test for an array $x$ in
both language implementations, with x being:

⬇

x = [0.1, 0.4, 0.76]

The python library can be used like the following:

⬇

>>> from spacings import rps

>>> rps(x, "uniform")

RPStestResult(statistic=0.9547378863245608, pvalue=0.8865399970192409)

and the Julia equivalent giving identical results in the following:

⬇

>>> using SpacingsTests

>>> rps(x, Uniform())

(statistic=0.9547378863245608, pvalue=0.8865399970192409)

##  4 Example Application 1: Bump Hunting

In this section, we illustrate how the RPS test could be used in a physics
scenario. We consider a detector that collects a number of events in an
observable $x$, where $x$ could for example be the energy of an event, the
detection time, or a reconstructed quantity like an invariant mass. We expect
some or all of the observed events to follow a known background distribution
$f_{B}(x)$, but there may be an additional contribution of events from an
unknown signal distribution $f_{S}(x)$—such as a rare, exotic particle decay
with unknown mass. Hence we want to quantify the goodness-of-fit of the
background only model to our data. A resulting low p-value could indicate the
presence of events distributed according to an additional, unknown signal
distribution.

In the example here, we use an exponential distribution $f_{B}(x)=e^{-x}$ for
the background model (null-hypothesis). In order to illustrate how the
presence of an actual signal (alternative hypothesis) would affect the
outcome, we also inject additional events following a normal distribution
centred at $x=1$ and width $\sigma=0.05$. The number of events is Poisson
fluctuated for both background and signal, with expected values of $\langle
n_{b}\rangle=100$ and $\langle n_{s}\rangle$ varied as specified. In Fig. 7,
an example distribution of observed events is shown, together with the assumed
background distribution, and the distribution with injected signal (here
$\langle n_{s}\rangle=5$).

The example case chosen is similar to that, for instance, of a search for an
exotic particle with unknown mass —a problem sometime referred to as "bump
hunting". In this case, $x$ would represent an invariant mass.

N.B., we do not assume that we know the rate of the underlying processes,
meaning that the number of observed counts is not included in our analysis
other than for the calculation of the test statistic. This means that we test
for the "shape" of the distribution, not its normalization. The conversion of
events via the CDF of the distribution under test $f_{B}$ transforms the
problem into a test of uniformity.

Figure 7: Example physics problem, with observed events distributed in $x$. We
test the goodness-of-fit of the background only model (blue) to the samples.
Here the samples have been generated according to a different distribution
with an injected signal (orange).

The p-value distributions under the assumption of $H^{0}$ (i.e. only
background is present) for repeated trials with $\langle n_{b}\rangle=100$,
and various injected $\langle n_{s}\rangle=[0,3,6,9,12,15]$ are shown in Fig.
8. All distributions with no signal ($\langle n_{s}\rangle=0$) show a flat
p-value distribution as expected, since in that case all events are drawn from
the background distribution $p_{B}$. For trials with injected signal, the
distributions are trending towards smaller p-values, indicating the worsened
goodness-of-fit for the background only model. In the example, all tests
exhibit this behaviour, while the RPS test offers the largest rejection
probability of the null hypothesis.

Figure 8: p-value distributions for background only samples ($\langle
n_{s}\rangle=0$) and background plus randomised signal injections comparing to
the background model for several choices of test statistics.

We quantify the sensitivity of the analysis to reject the background only
model at different significance levels under the assumption of the presence of
a signal. Therefore we check the median p-value of repeated trials, and at
what value of $\langle n_{s}\rangle$ it crosses specific critical values (See
left panel of Fig. 9). In our chosen example, for a signal of strength
$\langle n_{s}\rangle=10$ we expect to reject the background only model using
RPS at the $2\sigma$ significance level555A significance level in terms of
numbers of $k$ standard deviations $\sigma$ can be translated to a p-value as
one minus the integral over a unit normal distribution form $-k$ to $+k$.,
whereas for the other tests, a signal of at least $\langle n_{s}\rangle=20$ is
needed to achieve the same. Such a large signal of $\langle n_{s}\rangle=20$
would allow to reject the background only model at $>4\sigma$ significance
with the RPS test.

Figure 9: The expected significance level at which the background model can be
excluded under the assumption of a signal, as a function of $\langle
n_{s}\rangle$ for the different tests.

##  5 Example Application 2: Trigger for Transient Neutrino Emission

This section summarizes one of the first applications of the RPS test in
astrophysics, namely for triggering transient events in cryogenic neutrino
detectors [21] such as the RES-NOVA experiment [22]. While the technical
details about the experimental setup, the simulation and the application of
the RPS test can be found in the aforementioned references, here we will
summarize some highlights.

Cryogenic neutrino detectors can be described as counting experiments, that
output a temporal data stream of observed neutrino interactions. Without the
presence of a transient neutrino source, we only expect some activity from
background events. If a source of neutrinos is placed at an observable
distance, such as a core-collapse supernova (CC-SN) at 10 kpc, we expect a
short burst of neutrinos resulting in an excess in the observed counts over
the background only expectation. The sources of such transient neutrinos can
vary in their overall duration, temporal distribution and amplitude. Figure 10
shows as examples the expected counts of two different neutrino sources, a CC-
SN and a failed CC-SN, respectively, together with a constant background
expectation at a rate of 0.18 Hz.

To issue alerts in near real time about the presence of such sources, one
needs a triggering system with a chosen false alarm rate (FAR), which is set
to 1 per week for SNEWS [23]. The standard approach for building such triggers
is the usage of Poisson statistics, that analyse the data stream in windows of
a fixed length, and check the level of observed counts compared to the
expectation from background, see for example Ref. [24]. The Poisson approach
works well if the window size is chose optimally for a given signal. However,
if the chosen window size does not match the signal, the performance is
affected as either the signal is not contained in the window (window too
small), or the window is too large and the signal is washed out by background
events. Performance of the Poisson test, as a function of the background rate
and for the case of optimal window choice for the CC-SN and the failed CC-SN
signals is shown in Fig. 11.

The RPS test can likewise be used to analyze data streams to look for
transient phenomena. Here we do not make any explicit assumption on the
background rate, but rather assume that the background is constant in rate,
which means the distribution of background events in the time dimension is
following a uniform distribution. With RPS we can test for this uniformity,
which can be used to detect short additional contribution of events in the
data. The performance of the RPS test as trigger is also shown in Fig. 11.

In the case where the window size for the Poisson test is optimal, the
performance can not be matched with RPS (the panels in the upper left and
lower right, respectively, in Fig. 11) and results in up to 10% lower
sensitivity. However, the more interesting case is when using the window size
optimized for one signal for the analysis of a different signal (the panels in
the upper right and lower left, respectively, in Fig. 11). The RPS test is
more robust to such changes, and in the example of searching for a failed CC-
SN signal with a window optimized on a particular CC-SN scenario, we find up
to 20% increase in sensitivity.

In general, what we find is that the RPS test being non parametric and able to
deal with much larger analysis windows is more robust to changing conditions.
Less assumptions about the background rate and the expect signals have to be
made at the trade off of being non-optimal to one specific signal choice, but
good performance for the more agnostic case of unknown signal distributions.
This makes RPS an interesting choice for a general-purpose, agnostic trigger
algorithm for the search of transient events.

Figure 10: Example of observed counts at a neutrino detector for signals from
a core-collapse SN (at time $t=15$ s) and a failed core-collapse SN (at $t=5$
s) for progenitors stars with 27 $M_{\odot}$ and 40 $M_{\odot}$ respectively,
both at a distance of 10 kpc. (Modified version of a Figure from Ref. [21])
Figure 11: Maximum distance probed at a 95% success rate as a function of time
with respect to two sample signals, the Core-Collapse SN and the failed Core-
Collapse SN, obtained using analysis windows optimised on each of the tested
signals. (Figure from Ref. [21])

##  6 Performance Comparison

This section presents an in-depth performance comparison of the RPS test to
several other tests referenced in the introduction (KS, AD, CvM and Moran—all
those that allow to compute p-values). We are interested in detecting small
changes in an otherwise uniform distribution, and therefore construct the
following generic benchmark scenario: For one simulation of a specific test
case $H^{K}(n,s,w)$ we generate $(1-s)\cdot n$ random variates666Numbers of
samples are rounded to the closest integer from a standard uniform
distribution $\mathcal{U}(0,1)$, where $s$ is a signal fraction. In addition,
we include $s\cdot n$ samples distributed according to
$\Delta+\mathcal{U}(0,w)$ with the offset $\Delta=\mathcal{U}(0,1-w)$, i.e. a
more narrow uniform distribution of width $w$ over a random interval within
$(0,1)$. In our comparison, we vary all three parameters of $H^{K}(n,s,w)$,
i.e. the number of samples $n$, as well as the fraction $s$ and width $w$ of
the injected signal events. A sensitive test should be able to detect the
presence of the added, narrower signal samples by reporting a low p-value.

Figure 12: Comparison of the performance (median p-value of repeated trials,
individual panel’s xy-axis) as a function of the number of total samples
(large y-axis), the width of the signal (large x-axis) samples, and the
fraction of signal samples (individual panel’s x-axis). The number of signal
samples is rounded to the closest integer, hence the "step"-like features
visible mostly in the first few rows.

Figure 12 show the performance of our choice of tests as a function of the
above three parameters. As a metric, we show the median p-value obtained from
repeated trials, and we interpret a lower reported median p-value as a more
powerful test. This number can be interpreted as the median significance at
which we expect to be able to reject the null hypothesis. What can be observed
is, that for all the tested scenarios the RPS test is performing either on par
or significantly better than the Moran test. The EDF based tests (KS, AD or
CvM) start to dominate in terms of performance only for relatively wide
signals of around 25% total width or more. When analysing the goodness-of-fit
given a large number of samples, i.e. order of several hundreds, the
differences between RPS and the EDF-style tests start to become smaller.
Overall, the outcome of this performance study suggests that when signals are
expected of widths that span over less than a 25% percentile of the null
hypothesis distribution, and if the number of samples is $n<1000$, the RPS
tests compares very favourably against all others considered.

We also investigated other metrics to judge the test’s performance, such as
the area under the receiver operating characteristics (ROC) curve between
signal and null hypothesis trials. The overall picture does not change
substantially.

##  7 Conclusions

The RPS test statistic is a sensitive measure to detect deviations of samples
from a continuous distribution with known CDF. The analytic distribution of
the RPS statistic is not available for $n>1$, but a high accuracy
parameterization valid up to sample sizes of $n=1000$ is provided in order to
use RPS as a goodness-of-fit test. In the presented test scenarios, the RPS
test outperforms other tests significantly under certain circumstances, in
particular when the observed sample is small ($n<1000$) and introduced
deviations are narrow, i.e. concentrated over a small quantile. Two example
physics analysis cases were presented, we show that the sensitivity of a "bump
hunting" experiment could be boosted by up to a factor of two by choosing the
RPS test over others. And we show how RPS can be used to build a robust and
agnostic trigger algorithm for a SN experiment.

## Acknowledgements

We would like to thank Allen Caldwell, Oliver Schulz and Johannes Buchner for
helpful discussions and comments.
This research was supported by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) under Germany´s Excellence Strategy – EXC-2094 –
390783311

## References

  * [1] A. Kolmogorov, _Sulla determinazione empirica di una legge di distribuzione_ , _G. Ist. Ital._ 4 (1933) .
  * [2] N. Smirnov, _Table for estimating the goodness of fit of empirical distributions_ , _Annals of Mathematical Statistics_ 19 (1948) 279.
  * [3] T.W. Anderson and D.A. Darling, _A test of goodness of fit_ , _Journal of the American Statistical Association_ 50 (1954) 765.
  * [4] Y. Marhuenda, D. Morales and M.C. Pardo, _A comparison of uniformity tests_ , _Statistics_ 39 (2005) 315.
  * [5] H. Cramer, _On the composition of elementary errors_ , _Scandinavian Actuarial Journal_ 1928 (1928) 13.
  * [6] R.V. Mises, _Wahrscheinlichkeit Statistik und Wahrheit_ , Springer-Verlag Berlin Heidelberg (1928), 10.1007/978-3-662-36230-3.
  * [7] R. Pyke, _The Supremum and Infimum of the Poisson Process_ , _The Annals of Mathematical Statistics_ 30 (1959) 568 .
  * [8] J. Durbin, _Tests for serial correlation in regression analysis based on the periodogram of least-squares residuals_ , _Biometrika_ 56 (1969) 1.
  * [9] H.D. Brunk, _On the Range of the Difference between Hypothetical Distribution Function and Pyke’s Modified Empirical Distribution Function_ , _The Annals of Mathematical Statistics_ 33 (1962) 525 .
  * [10] R.C.H. Cheng and M.A. Stephens, _A goodness-of-fit test using moran’s statistic with estimated parameters_ , _Biometrika_ 76 (1989) 385.
  * [11] M. Greenwood, _The statistical study of infectious diseases_ , _Journal of the Royal Statistical Society_ 109 (1946) 85.
  * [12] N. Cressie, _On the logarithms of high-order spacings_ , _Biometrika_ 63 (1976) 343.
  * [13] N. Cressie, _An optimal statistic based on higher order gaps_ , _Biometrika_ 66 (1979) 619.
  * [14] L. Shtembari and A. Caldwell, _On the sum of ordered spacings_ , 2020.
  * [15] B. Efron, _Bootstrap Methods: Another Look at the Jackknife_ , _Annals Statist._ 7 (1979) 1.
  * [16] H.A. David and H.N. Nagaraja, _Order statistics_ , Wiley (2003).
  * [17] P. Dierckx, _An algorithm for smoothing, differentiation and integration of experimental data using spline functions_ , _Journal of Computational and Applied Mathematics_ 1 (1975) 165.
  * [18] P. Dierckx, _A fast algorithm for smoothing data on a rectangular grid while using spline functions_ , _SIAM Journal on Numerical Analysis_ 19 (1982) 1286 [https://doi.org/10.1137/0719093].
  * [19] P. Dierckx, _Curve and surface fitting with splines_ , in _Monographs on numerical analysis_ , 1996.
  * [20] F.N. Fritsch and J. Butland, _A method for constructing local monotone piecewise cubic interpolants_ , _SIAM Journal on Scientific and Statistical Computing_ 5 (1984) 300 [https://doi.org/10.1137/0905021].
  * [21] P. Eller, N. Iachellini Ferreiro, L. Pattavina and L. Shtembari, _Online triggers for supernova and pre-supernova neutrino detection with cryogenic detectors_ , 2205.03350.
  * [22] RES-NOVA collaboration, _RES-NOVA sensitivity to core-collapse and failed core-collapse supernova neutrinos_ , _JCAP_ 10 (2021) 064 [2103.08672].
  * [23] SNEWS collaboration, _SNEWS 2.0: a next-generation supernova early warning system for multi-messenger astronomy_ , _New J. Phys._ 23 (2021) 031201 [2011.00035].
  * [24] N.Y. Agafonova et al., _On-line recognition of supernova neutrino bursts in the LVD detector_ , _Astropart. Phys._ 28 (2008) 516 [0710.0259].