NewReferences NewReferences

# A magnetar giant flare in the nearby starburst galaxy M82

Sandro Mereghetti1,∗ Michela Rigoselli1 Ruben Salvaterra1 Dominik Patryk
Pacholski1,2 James Craig Rodi3 Diego Gotz4 Edoardo Arrigoni5,1 Paolo D’Avanzo6
Christophe Adami7 Angela Bazzano3 Enrico Bozzo8,6 Riccardo Brivio9,6 Sergio
Campana6 Enrico Cappellaro10 Jerome Chenevez11 Fiore De Luise12 Lorenzo
Ducci13,8 Paolo Esposito14,1 Carlo Ferrigno8,6 Matteo Ferro9,6 Gian Luca
Israel15 Emeric Le Floc’h4 Antonio Martin-Carrillo16 Francesca Onori12 Nanda
Rea17,18 Andrea Reguitti6,10 Volodymyr Savchenko19 Leonardo Tartaglia12
William Thuillot20 Andrea Tiengo14,1 Lina Tomasella10 Martin Topinka21 Damien
Turpin4 Pietro Ubertini3

###### Abstract

Giant flares, short explosive events releasing up to 1047 erg of energy in the
gamma-ray band in less than one second, are the most spectacular manifestation
of magnetars, young neutron stars powered by a very strong magnetic field,
1014-15 G in the magnetosphere and possibly higher in the star interior[1, 2].
The rate of occurrence of these rare flares is poorly constrained, as only
three have been seen from three different magnetars in the Milky Way[3, 4] and
in the Large Magellanic Cloud[5] in $\sim$50 years since the beginning of
gamma-ray astronomy. This sample can be enlarged by the discovery of
extragalactic events, since for a fraction of a second giant flares reach peak
luminosities above 1046 erg s-1, which makes them visible by current
instruments up to a few tens of Mpc. However, at these distances they appear
similar to, and difficult to distinguish from, regular short gamma-ray bursts
(GRBs). The latter are much more energetic events, 1050-53 erg, produced by
compact binary mergers and originating at much larger distances[6]. Indeed,
only a few short GRBs have been proposed[7, 8, 9, 10, 11], with different
levels of confidence, as magnetar giant flare candidates in nearby galaxies.
Here we report the discovery of a short GRB positionally coincident with the
central region of the starburst galaxy M82[12]. Its spectral and timing
properties, together with the limits on its X-ray and optical counterparts
obtained a few hours after the event and the lack of an associated
gravitational wave signal, qualify with high confidence this event as a giant
flare from a magnetar in M82.

INAF - Istituto di Astrofisica Spaziale e Fisica Cosmica di Milano, Via A.
Corti 12, 20133 Milano, Italy

Università degli Studi di Milano Bicocca, Dipartimento di Fisica G.
Occhialini, Piazza della Scienza 3, 20126 Milano, Italy

INAF - Istituto di Astrofisica e Planetologia Spaziali di Roma, Via del Fosso
del Cavaliere 100, 00133 Roma, Italy

Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-
sur-Yvette, France

Università degli Studi di Milano, Dipartimento di Fisica, Via Celoria 16,
20133 Milano, Italy

INAF - Osservatorio Astronomico di Brera, Via E. Bianchi 46, 23807 Merate,
Italy

Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

University of Geneva, Department of Astronomy, Chemin d’Ecogia 16, Versoix
1290, Switzerland

Università dell’Insubria, Dipartimento di Scienza e Alta Tecnologia, Via
Valleggio 11, 22100 Como, Italy

INAF - Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122
Padova, Italy

DTU Space, Technical University of Denmark, Elektrovej 327, 2800 Kongens
Lyngby, Denmark

INAF - Osservatorio Astronomico d’Abruzzo, Via M. Maggini snc, 64100 Teramo,
Italy

Institut fuer Astronomie und Astrophysik Tuebingen, Sand 1, Tuebingen, Germany

Scuola Universitaria Superiore IUSS Pavia, Piazza della Vittoria 15, 27100
Pavia, Italy

INAF - Osservatorio Astronomico di Roma, 00078 Monte Porzio Catone, Italy

School of Physics and Centre for Space Research, University College Dublin,
Belfield, Dublin 4, Dublin, Ireland

Institute of Space Sciences (CSIC-ICE), Campus UAB, Carrer de Can Magrans s/n,
08193 Barcelona, Spain

Institut d’Estudis Espacials de Catalunya (IEEC), Carrer Gran Capità 2–4,
08034 Barcelona, Spain

École polytechnique fédérale de Lausanne, Lausanne, Switzerland.

Institut de mecanique celeste et de calcul des ephemerides (IMCCE) UMR 8028 du
CNRS - Observatoire de Paris 77 Av. Denfert Rochereau 75014 Paris, France

INAF - Osservatorio Astronomico di Cagliari, Via della Scienza 5, 09047
Selargius (CA), Italy

∗ Corresponding author<EMAIL_ADDRESS>
Instruments on board the INTEGRAL satellite detected GRB 231115A on 2023
November 15 at 15:36:20.7 UT[13]. This short burst (duration $T_{90}$ = 93 ms)
was seen also by other satellites[14, 15, 16, 17, 18], but only INTEGRAL could
promptly associate it to M82 thanks to a localization at few arcmin level,
publicly distributed by the INTEGRAL Burst Alert System (IBAS[19]) only 13 s
after the burst detection. The burst occurred in the field of view of the IBIS
coded mask instrument[20], at coordinates R.A.=149.0205 deg, Dec.=+69.6719 deg
(J2000, 2 arcmin 90% c.l. radius). This position is consistent with, and
supersedes, the one derived using the preliminary satellite attitude available
in near real time[13]. The position of GRB 231115A coincides with the nearby
starburst galaxy M82[12] (see Figure 1). Notably, the central region of the
galaxy, where most star formation activity occurred, is inside the
INTEGRAL/IBIS error region. Considering that the total angular size of all the
galaxies with apparent luminosity brighter than M82 (excluding the Magellanic
Clouds and M31) is $\sim$6000 arcmin2, the chance alignment with GRB 231115A
is 4$\times 10^{-5}$. A more elaborate analysis[21] indicates that the Bayes
factor favouring a giant flare in M82 with respect to a chance alignment of a
short GRB in the background is a factor 30 larger than that of the previous
best candidate (GRB 200415A possibly associated to NGC 253[9]), making GRB
231115A by far the most compelling case for a magnetar giant flare outside the
local group of galaxies.

The light curves plotted in Figure 2 show that GRB 231115A was clearly visible
in the two IBIS detectors, which operate in different energy ranges: ISGRI[22]
(20 keV – 1 MeV) and PICsIT[23] (175 keV – 15 MeV). The ISGRI and PICsIT
spectra of the burst, extracted from the time interval with a significant
detection in both detectors (from T0+0.687 s to T0+0.748 s, with T0=15:36:20.0
UT), have been jointly fit and are well described by an exponentially cut-off
power law with photon index $\alpha=0.04^{+0.27}_{-0.24}$, peak energy
$E_{p}=551^{+81}_{-59}$ keV, and flux $F_{30-2600}=(7.2^{+0.6}_{-0.7})\times
10^{-6}$ erg cm-2 s-1 in the 30 keV – 2.6 MeV range (see Methods sect. 1).
Assuming this spectral shape for the whole duration of the burst, the fluence
in the same energy range is $(6.3\pm 0.5)\times 10^{-7}$ erg cm-2.

M82 was in the field of view of IBIS starting from 9.2 hours before up to 11.8
minutes after the burst, as well as during a target of opportunity (ToO)
observation performed in the following three days, but no other bursting or
persistent emission besides GRB 231115A was detected. Integrating all the data
of the ToO observation, we obtained an upper limit of $2.6\times 10^{-11}$ erg
cm-2 s-1 on the 30–100 keV persistent emission at the position of GRB 231115A.
Also, no other bursts from M82 were detected in all the available ISGRI
archival data obtained since 2003, which provided a total exposure of more
than 16 Ms (see Methods sect. 1).

At the distance of M82 (3.6 Mpc[24]), the GRB 231115A fluence measured with
IBIS implies an emitted isotropic energy $E_{\rm iso}$ = $10^{45}$ erg (1 keV
– 10 MeV), well below the typical value for a short GRB, but consistent with
the one of the initial pulses of the three giant flares securely associated
with magnetars[3, 4, 5]. The short duration and hard spectrum of GRB 231115A
are also in agreement with the properties seen in the initial pulses of
magnetar giant flares. In the three giant flares securely associated with
magnetars, the initial short and hard pulse was followed by a softer tail,
with a duration of a few minutes and a maximum luminosity of $\sim\\!10^{42}$
erg s-1, characterized by a periodic modulation induced by the neutron star
rotation. A similar feature would be too faint to be visible by IBIS at the
distance of M82. This is shown in Figure 3, where we compare the limits
obtained for GRB 231115A with the light curve of the most energetic Galactic
giant flare ever observed (the one emitted by SGR 1806–20 in 2004[4]) rescaled
to a distance of 3.6 Mpc.

The field of GRB 231115A was promptly observed with the Neil Gehrels Swift
Observatory starting 9 ks after the burst[25] and a deeper X-ray observation
was carried out with the XMM-Newton satellite, starting 0.7 days after the
burst[26] (see Methods sect. 2). The comparison of the Swift and XMM-Newton
images with those obtained in previous observations of M82 does not show
evidence for new X-ray sources inside the GRB 231115A error region. Due to the
presence of unresolved X-ray emission from the central part of the galaxy, the
corresponding upper limit is position-dependent. The limit of $4\times
10^{-14}$ erg cm-2 s-1 derived from the XMM-Newton data, which applies to 60%
of the error region, is inconsistent with the flux at T0+1 day of the large
majority of X-ray afterglows of short GRBs detected with the Swift/XRT
instrument (see Methods Figure 7). Although we can not exclude a rapid
afterglow drop, as observed in the class of so-called short lived short
GRBs[27], this X-ray upper limit disfavors the GRB interpretation. It is
instead consistent with the quiescent X-ray emission of a magnetar in M82. In
fact the X-ray luminosity of known magnetars, also considering the peak values
reached during outbursts[28], never exceeds 1036 erg s-1, corresponding to
$6\times 10^{-16}$ erg cm-2 s-1 at the M82 distance.

Holding true the association of GRB 231115A with M82 on the basis of the
remarkable spatial coincidence, the low value of $E_{\rm iso}$ could be
reconciled with a short GRB origin only if the jet axis is pointed away from
our direction, similar to the case of the gravitational wave event GRB
170817A[29], which produced the bright kilonova AT2017gfo in the optical/NIR
band. However, the deep upper limits of m $>$ 20.0 – 24.0 mag (depending on
the assumed position in the INTEGRAL error circle) on the optical counterpart
that we obtained starting at T0+0.2 days (see Methods sect. 3) exclude this
possibility, both for an AT2017gfo-like event and for even fainter kilonovae
(see Methods sect. 3). We finally note that the binary merger of two compact
objects at 3.6 Mpc would have produced a strong signal in gravitational waves,
at variance with the non-detection reported by the LIGO/Virgo/KAGRA
Collaboration[30].

The discovery of a young, active magnetar in M82, a starburst galaxy
characterized by a high star formation rate[31] is consistent with the origin
of magnetars in core collapse supernova explosions[32]. The volumetric rate of
giant flares with $E_{\rm iso}>4\times 10^{44}$ erg has been recently
estimated[10] as $(3.8^{+4.0}_{-3.1})\times 10^{5}$ Gpc-3 yr-1. Assuming a
direct link with the young stellar population of core collapse progenitors and
considering the total star formation rate of 4000 $M_{\odot}$ yr-1 within 50
Mpc, this corresponds to a rate of $R(E_{iso}>4\\!\times
10^{44})=0.05^{+0.05}_{-0.04}$ yr-1 ($M_{\odot}$ yr-1)-1. Given the star
formation rate of 7.1 $M_{\odot}$ yr-1 in M82[31] and assuming a power-law
distribution of the giant flares energies with slope 1.7, the expected rate of
giant flares with $E_{\rm iso}>10^{45}$ erg in M82 is $R_{\rm
M82}(E_{iso}>10^{45})=0.19^{+0.19}_{-0.15}$ yr-1. Recent calculations based on
relativistic hydrodynamical simulations[33] show that such ejecta are
efficient sources for the nucleosynthesis of heavy elements through the
r-process. Thus the giant flares with $E_{\rm iso}>10^{46}$ erg (each
producing up to $10^{26}$ g of ejecta) can give a yield of $\sim 2\times
10^{-9}$ $M_{\odot}$ yr-1 of r-process elements in M82 through this channel.

GRB 051103 is another short GRB that has been proposed as a possible magnetar
giant flare in the M81 group of galaxies[7]. According to a detailed
statistical analysis[10], M82 is the most likely host of GRB 051103, despite
being slightly outside its 260 arcmin2 localization region. The rate of giant
flares derived above is fully consistent with the detection of two events from
M82 in 20 years. Therefore, starburst galaxies such as M82, which produce a
significant population of young magnetars, appear as promising targets to
constrain the energy distribution function of giant flares by means of long
dedicated observations with future high sensitivity instruments.

## References

## References

  * [1] Mereghetti, S. The strongest cosmic magnets: soft gamma-ray repeaters and anomalous X-ray pulsars. _A &A Rev._, 15, 225 (2008).
  * [2] Kaspi, V. M., & Beloborodov, A. M. Magnetars. _ARA &A_, 55, 261 (2017).
  * [3] Hurley, K., et al. A giant periodic flare from the soft $\gamma$-ray repeater SGR1900+14. _Nature_ , 397, 41 (1999).
  * [4] Palmer, D. M., et al. A giant $\gamma$-ray flare from the magnetar SGR 1806-20. _Nature_ , 434, 1107 (2005).
  * [5] Mazets, E. P., Golentskii, S. V., Ilinskii, V. N., Aptekar, R. L., & Guryan, I. A. Observations of a flaring X-ray pulsar in Dorado. _Nature_ , 282, 587 (1979).
  * [6] D’Avanzo, P. Short gamma-ray bursts: A review. _Journal of High Energy Astrophysics_ , 7, 73 (2015).
  * [7] Frederiks, D. D., et al. On the possibility of identifying the short hard burst GRB 051103 with a giant flare from a soft gamma repeater in the M81 group of galaxies. _Astronomy Letters_ , 33, 19 (2007).
  * [8] Mazets, E. P., et al. A Giant Flare from a Soft Gamma Repeater in the Andromeda Galaxy (M31). _ApJ_ , 680, 545 (2008).
  * [9] Svinkin, D., et al. A bright $\gamma$-ray flare interpreted as a giant magnetar flare in NGC 253. _Nature_ , 589, 211 (2021).
  * [10] Burns, E., et al. Identification of a Local Sample of Gamma-Ray Bursts Consistent with a Magnetar Giant Flare Origin. _ApJ_ , 907, L28 (2021).
  * [11] Trigg, A. C., et al. GRB 180128A: A Second Magnetar Giant Flare Candidate from the Sculptor Galaxy. _arXiv e-prints_ , arXiv:2311.09362 (2023).
  * [12] Förster Schreiber, N. M., Genzel, R., Lutz, D., & Sternberg, A. The Nature of Starburst Activity in M82. _ApJ_ , 599, 193 (2003).
  * [13] Mereghetti, S., et al. GRB 231115A: a short hard GRB detected by IBAS, poistionally coincident with M82. _GRB Coordinates Network_ , 35037, 1 (2023).
  * [14] Dalessi, S., Roberts, O. J., Veres, P., Meegan, C., & Fermi Gamma-ray Burst Monitor Team GRB 231115A: Fermi Observations of a probable Magnetar Giant Flare from M82. _GRB Coordinates Network_ , 35044, 1 (2023).
  * [15] Cheung, C. C., et al. GRB 231115A (short): Glowbug gamma-ray detection. _GRB Coordinates Network_ , 35045, 1 (2023).
  * [16] Xue, W. C., Xiong, S. L., Li, X. B., Li, C. K., & Insight-HXMT Team GRB 231115A: Insight-HXMT/HE detection. _GRB Coordinates Network_ , 35060, 1 (2023).
  * [17] Frederiks, D., et al. Konus-Wind detection of GRB 231115A (a probable Magnetar Giant Flare from M82). _GRB Coordinates Network_ , 35062, 1 (2023).
  * [18] Wang, Y., et al. Study the origin of GRB 231115A, short gamma-ray burst or magnetar giant flare?. _arXiv e-prints_ , arXiv:2312.02848 (2023).
  * [19] Mereghetti, S., Götz, D., Borkowski, J., Walter, R., & Pedersen, H. The INTEGRAL Burst Alert System. _A &A_, 411, L291 (2003).
  * [20] Ubertini, P., et al. IBIS: The Imager on-board INTEGRAL. _A &A_, 411, L131 (2003).
  * [21] Burns, E. GRB 231115A: significance of INTEGRAL localization alignment with M82. _GRB Coordinates Network_ , 35038, 1 (2023).
  * [22] Lebrun, F., et al. ISGRI: The INTEGRAL Soft Gamma-Ray Imager. _A &A_, 411, L141 (2003).
  * [23] Labanti, C., et al. The Ibis-Picsit detector onboard Integral. _A &A_, 411, L149 (2003).
  * [24] Freedman, W. L., et al. The Hubble Space Telescope Extragalactic Distance Scale Key Project. I. The Discovery of Cepheids and a New Distance to M81. _ApJ_ , 427, 628 (1994).
  * [25] Osborne, J. P., et al. GRB 231115A: Swift-XRT and Swift-UVOT observations. _GRB Coordinates Network_ , 35064, 1 (2023).
  * [26] Rigoselli, M., Pacholski, D. P., Mereghetti, S., Salvaterra, R., & Campana, S. GRB 231115A: XMM-Newton observation. _GRB Coordinates Network_ , 35175, 1 (2023).
  * [27] Sakamoto, T., & Gehrels, N. Indication of Two Classes in the Swift Short Gamma-Ray Bursts from the XRT X-Ray Afterglow Light Curves. _Gamma-ray Burst: Sixth Huntsville Symposium_ , 1133, 112 (2009).
  * [28] Coti Zelati, F., Rea, N., Pons, J. A., Campana, S., & Esposito, P. Systematic study of magnetar outbursts. _MNRAS_ , 474, 961 (2018).
  * [29] Abbott, B. P., et al. Multi-messenger Observations of a Binary Neutron Star Merger. _ApJ_ , 848, L12 (2017).
  * [30] Ligo Scientific Collaboration, VIRGO Collaboration, & Kagra Collaboration GRB 231115A: Non-detection in low-latency of gravitational waves with LIGO/Virgo/KAGRA. _GRB Coordinates Network_ , 35049, 1 (2023).
  * [31] Leroy, A. K., et al. A z = 0 Multiwavelength Galaxy Synthesis. I. A WISE and GALEX Atlas of Local Galaxies. _ApJS_ , 244, 24 (2019).
  * [32] Duncan, R. C., & Thompson, C. Formation of Very Strongly Magnetized Neutron Stars: Implications for Gamma-Ray Bursts. _ApJ_ , 392, L9 (1992).
  * [33] Cehula, J., Thompson, T. A., & Metzger, B. D. Dynamics of baryon ejection in magnetar giant flares: implications for radio afterglows, r-process nucleosynthesis, and fast radio bursts. _arXiv e-prints_ , arXiv:2311.05681 (2023).
  * [34] Frederiks, D. D., et al. Giant flare in SGR 1806-20 and its Compton reflection from the Moon. _Astronomy Letters_ , 33, 1 (2007).
  * [35] Mereghetti, S., et al. The First Giant Flare from SGR 1806-20: Observations Using the Anticoincidence Shield of the Spectrometer on INTEGRAL. _ApJ_ , 624, L105 (2005).

Figure 1: Optical image of M82 (from data obtained with the TNG). RGB scale:
z-band (red), i-band (green), r-band (blue). The 90% c.l. error circle of GRB
231115A has a radius of 2 arcmin.  Figure 2: Light curves of GRB 231115A.
Background-subtracted light curves (errors are at 1$\sigma$) obtained with the
ISGRI detector in the $30-250$ keV energy range (top panel) and with the
PICsIT detector in the $312-2600$ keV energy range (bottom panel). Time is
referred to T0 = 2023-11-15 15:36:20 UT (time at the INTEGRAL position; the
burst reached instruments on low Earth orbit satellites about 0.4 s later).
Figure 3: Comparison between GRB 231115A and the 2004 giant flare from SGR
1806–20. Errors are at 1$\sigma$, upper limits are at 3$\sigma$. The light
curve of SGR 1806–20 (red) has been rescaled to the distance of M82 (data from
Helicon-Coronas-F[34], Konus-Wind[34] and INTEGRAL SPI-ACS[35]). Note that the
5.2 s periodicity in the SGR 1806–20 light curve is not visible due to the
used bin size.

##  1 INTEGRAL data analysis and results

IBIS is a $\gamma$-ray telescope based on two position-sensitive detection
planes working in different overlapping energy ranges: ISGRI[22] (20 keV – 1
MeV) and PICsIT[23] (175 keV – 15 MeV). A coded mask placed about 3.2 m above
ISGRI provides imaging capabilities over a field of view of $\sim$30$\times$30
deg2. GRB 231115A was detected at an off-axis angle of 4∘, inside the fully
coded field of view. For the data reduction and analysis, we used version 11.2
of the Offline Science Analysis (OSA) software[1], with the most recent
calibration files, that properly account for the time evolution of the
instrument response.

###  1.1 Burst spectrum

We extracted an ISGRI spectrum, in 14 logarithmically spaced energy channels
between $30-500$ keV, for the time interval from 15:36:20.687 to 15:36:20.748
UT. The PICsIT spectrum for the same time interval was extracted from the
spectral-timing data mode in 8 predefined energy channels spanning $212-2600$
keV. Data in spectral-timing mode (counts integrated over the whole detector
in 3.9 ms time bins) do not provide imaging information (contrary to the ISGRI
data). Therefore, the background count-rate in each PICsIT energy channel was
calculated as the median count rate during the INTEGRAL pointing (5013 s
duration) containing the burst. A proper correction factor was applied to the
PICsIT effective area to account for the angular distance between the IBIS
pointing direction and GRB 231115A.

We added 5% systematic uncertainties to both spectra. The ISGRI and PICsIT
spectra were fitted simultaneously using XSPEC[2] version 12.12.1. A fit with
a power law was unacceptable ($\chi^{2}=59.4$ for 17 degrees of freedom),
while a power law with exponential cut-off, defined as
$F(E)=KE^{\alpha}\,\exp(-E(2+\alpha)/E_{p})$, gave a good fit with photon
index $\alpha=0.04^{+0.27}_{-0.24}$, peak energy $E_{p}=551^{+81}_{-59}$ keV,
and $30-2600$ keV flux $F$ = $(6.4\pm 0.6)\times 10^{-6}$ erg cm-2 s-1.

The ISGRI spectrum for the time interval corresponding to the whole GRB, from
15:36:20.687 UT to 15:36:20.833 UT, is less well constrained, but the best fit
parameters are within the errors consistent with those of the ISGRI plus
PICsIT spectrum. All the fit parameters are summarized in Table 1, where we
give also the results obtained with a blackbody model.

Table 1: Results of spectral fits of GRB 231115A. Errors are at 1$\sigma$, fluxes are in the $30-2600$ keV energy range. Time intervals are in seconds after T0 = 15:36:20 UT.  Detectors | ISGRI+PICsIT | ISGRI
---|---|---
Time interval | 0.687–0.748 | 0.687–0.833
Blackbody
$kT$ (keV) | $108^{+7}_{-5}$ | $68^{+12}_{-11}$
Flux ($10^{-6}$ erg cm-2 s-1) | $7.9^{+0.6}_{-0.7}$ | $2.1^{+0.6}_{-0.7}$
$\chi^{2}$ / d.o.f. | $31.9/17$ | $16.4/12$
Cut-off power law
$\alpha$ | $0.04^{+0.27}_{-0.24}$ | $0.55^{+0.50}_{-0.42}$
$E_{p}$ (keV) | $551^{+81}_{-59}$ | $355^{+158}_{-89}$
Flux ($10^{-6}$ erg cm-2 s-1) | $7.2^{+0.6}_{-0.7}$ | $2.7\pm 1.1$
$\chi^{2}$ / d.o.f. | $17.2/16$ | $14.9/11$

###  1.2 Follow up observation.

Immediately after the discovery of GRB 231115A we requested an INTEGRAL Target
of Opportunity (ToO) observation that started only 21 hours after the burst
and continued in the next satellite revolution. The first part of the ToO was
done between 12:38 UT and 19:46 UT of November 16. The next one from 07:21 UT
of November 17 to 23:09 UT of November 18, resulting in a total exposure of
162 ks divided in 57 pointings. Using the OSA 11.2 software, we produced the
mosaics of the ISGRI images in the $30-100$ keV energy range for three
periods: November 16, November 17-18, and November 16 to 18. No sources were
detected at the position of GRB 231115A. The 3-$\sigma$ upper limits on the
$30-100$ keV flux, derived assuming a thermal bremsstrahlung spectrum with
temperature $kT=30$ keV, are $6.5\times 10^{-11}$ erg cm-2 s-1 for November
16, $2.9\times 10^{-11}$ erg cm-2 s-1 for November 17-18, and $2.6\times
10^{-11}$ erg cm-2 s-1 for the sum of the two periods. We also derived an
upper limit of $3.8\times 10^{-10}$ erg cm-2 s-1 on the flux in the time
interval from T0 \+ 1.7 s to T0 \+ 11.8 min, during which the GRB 231115A
position was in the IBIS field of view. All these upper limits are plotted in
Figure 3, where it can be seen that, even a giant flare as energetic as the
one emitted from SGR 1806–20, at the distance of M82 would have been
undetectable by IBIS after the initial bright pulse.

###  1.3 Search for past activity from M82.

The region of M82 has been repeatedly observed with INTEGRAL, starting in
November 2003. After eliminating the time intervals with strong and variable
background (typically due to solar activity or to particles trapped in the
Earth radiation belts when the satellite is close to perigee) we obtained
about 16 millions of seconds of useful time with M82 in the IBIS field of
view. We extracted ISGRI lightcurves in the $30-150$ keV energy range, using
only detector pixels illuminated for more than 50% from a source at the GRB
231115A position and searched for excess counts in the light curves over eight
logarithmically spaced timescales between 0.01 and 1.28 s (see [3] for more
details). All the burst candidates found in the light curves were then
examined by making the sky images of the corresponding time intervals. None of
them showed the presence of a source at the M82 position. Through simulations
we found that bursts down to a factor $\sim$5 fainter than GRB 231115A would
have been detected in these data, assuming the same time profile and spectrum.

##  2 X-ray observations

Swift pointed at M82, starting to collect data 9.0 ks after the burst
occurrence[25]. Swift/XRT observations were carried out in Photon Counting
mode in the $9.0-39.2$ ks time interval, collecting 4.4 ks of data. No new
X-ray sources were detected within the INTEGRAL error region of GRB 231115A.
The 3-$\sigma$ upper limit is position-dependent due to the diffuse X-ray
emission from the M82 galaxy. Outside the galaxy, the upper limit on the count
rate is $\sim 2-3\times 10^{-3}$ counts s-1. Assuming a power-law spectrum
with photon index $\Gamma=2$ and the Galactic column density of $6.5\times
10^{20}$ cm-2 (the total Galactic absorption column in the line of sight of
M82[4]) this corresponds to a $0.3-10$ keV unabsorbed flux of $\sim 10^{-13}$
erg cm-2 s-1.

A 47 ks long ToO observation of GRB 231115A was carried out with the XMM-
Newton satellite, starting on 2023 November 16 at 08:29:16 UT, about 16.9
hours after the burst (Obs.ID 093239)[26]. The EPIC-pn[5] and the two EPIC-
MOS[6] cameras were operated in Full window mode, with the thin optical-
blocking filter. We processed the data with the epproc and emproc pipelines of
version 18 of the Science Analysis System (sas)[7] and the most recent
calibration files and selected the EPIC events with standard filtering
expressions. Time intervals of high background were removed with the espfilt
task with standard parameters, thus yielding net exposure times of 7.56 ks
(pn), 18.78 ks (MOS1) and 23.41 ks (MOS2). We selected single- and multiple-
pixel events (pattern$\leq$4 for the pn and $\leq$12 for the MOS cameras);
out-of-time events were also removed following the standard
procedureaaahttp://www.cosmos.esa.int/web/xmm-newton/sas-thread-epic-oot.

Figure 4: EPIC-pn images of M82. The exposure-corrected images refer to the
$0.3-2$ keV (left panel) and $2-10$ keV (right panel) energy ranges. The 90%
c.l. error circle of GRB 231115A has a radius of 2 arcmin.

Figure 4 shows the pn exposure-corrected images in the soft ($0.3-2$ keV) and
hard ($2-10$ keV) energy ranges. The images are dominated by the diffuse X-ray
emission of the M82 galaxy, especially in the soft X-ray band. We compared the
images obtained in this observation with those of all the previous pointings
on M82 available in the XMM-Newton archive (Obs.ID 011229, 020608, 056059,
065780, 087094, 089106). We did not find any evidence for the appearance of a
new source. In our observation, the total emission between $2-10$ keV from a
circular region of radius 1.5 arcmin at the center of the galaxy (coordinates
R.A. = ${09}^{\rm h}$ ${55}^{\rm m}$ ${50}^{\rm s}{\\!.9}$, Dec. =
${+69}^{\circ}$ ${40}^{\prime}$ ${47}^{\prime\prime}{}$; J2000) can be
described by an absorbed power law with $N_{\rm H}=(1.0\pm 0.1)\times 10^{22}$
cm-2, photon index $\Gamma=2.01\pm 0.04$ and unabsorbed flux $F_{2-10\,\rm
keV}=(2.88\pm 0.04)\times 10^{-11}$ erg cm-2 s-1. The flux measured in
previous XMM-Newton observations from the same extraction region and the same
spectral model varied between $\sim\\!1.2\times 10^{-11}$ erg cm-2 s-1 and
$\sim\\!5.4\times 10^{-11}$ erg cm-2 s-1, indicating the presence of one or
more variable sources that cannot be individually resolved with the XMM-Newton
spatial resolution.

To compute the upper limit for a new point source in the GRB 231115A error
box, we applied the eupper task to the $2-10$ keV images. We used a circle of
radius $15^{\prime\prime}$ for the source extraction and a concentric annulus
of radii $22.5^{\prime\prime}$ and $30^{\prime\prime}$ for the background. The
derived upper limits have been corrected for the fraction of source counts
falling outside the extraction circle. Figure 5 shows the 3-$\sigma$ upper
limits obtained by applying this procedure on a grid of positions with a step
of $5^{\prime\prime}$.

We used the appropriate EPIC response matrices to convert the count rates of
each camera to fluxes in the $0.3-10$ keV and $2-10$ keV energy ranges, with
the assumption of an absorbed power-law spectrum with photon index $\Gamma=2$
and $N_{\rm H}=6.5\times 10^{20}$ cm-2 Finally we combined the upper limits of
the three cameras to obtain an EPIC upper limit for each energy band. The
results are shown in Figure 6. About 60% of the error circle has a 3-$\sigma$
upper limit of $F_{2-10\,\rm keV}=1.2\times 10^{-14}$ erg cm-2 s-1 and
$F_{0.3-10\,\rm keV}=4.0\times 10^{-14}$ erg cm-2 s-1. In less than 10% of the
error circle, the upper limit is worse than $7\times 10^{-14}$ erg cm-2 s-1
($2-10$ keV) and $1.6\times 10^{-13}$ erg cm-2 s-1 ($0.3-10$ keV). At the
distance of M82, these fluxes correspond to luminosities of $1.9\times
10^{36}$ and $6.2\times 10^{36}$ erg s-1, respectively.

Figure 5:  Maps of count rate upper limits. The figure gives 3-$\sigma$ upper
limits on the $2-10$ keV count rates of the MOS1 (left panel), MOS2 (central
panel) and pn (right panel) cameras.  Figure 6: Maps of flux upper limits. The
figure gives 3-$\sigma$ upper limits on the fluxes in the $2-10$ keV (left
panel) and $0.3-10$ keV (right panel) energy range, obtained by combining the
three maps of Figure 5.  Figure 7: X-ray light curves of short GRB afterglows.
The Swift/XRT (black square) and XMM-Newton/EPIC (black diamond) 3-$\sigma$
upper limits of GRB 231115A are indicated.

Figure 7 shows the comparison between the X-ray light curves of short GRBs
(defined such as T${}_{90}\leq 2$ s), obtained from the Swift GRB
catalogbbbhttp://swift.gsfc.nasa.gov/archive/grb_table and the 3-$\sigma$
upper limits we derived from the Swift and XMM-Newton observations of GRB
231115A.

##  3 Optical observations

Multi-filter follow-up optical observations of GRB 231115A were carried out
with the wide-field Schmidt telescopes sited in INAF observatories of Padova
(Asiago, Italy), Abruzzo (Campo Imperatore, Italy) and with the 3.6-m
Telescopio Nazionale Galileo (TNG, Canary Islands, Spain) between about 5 and
12 hours from the event T0. We also took additional VRI images about 7 hours
after T0 with the 120 cm Newton telescope located at the Observatoire de Haute
Provence (OHP, France). The log of the observations is reported in Table 2.
Image reduction was carried out following the standard procedures: subtraction
of an averaged bias frame and division by a normalized flat frame. Astrometry
was performed using the Pan-
STARRSccchttp://outerspace.stsci.edu/display/PANSTARRS/ catalogue. We carried
out image subtraction with respect to SDSSdddhttp://www.sdss.org and telescope
archival templates using the HOTPANTS (High Order Transform of Psf ANd
Template Subtraction code[8]) package to find and pinpoint variable sources
within the INTEGRAL error circle. Aperture and PSF-matched photometry were
performed using the DAOPHOT package[9] and the STDPipe[10] and the
photutilseeehttp://github.com/astropy/photutils [11] packages for OHP/T120
images. To minimize any systematic effect, we performed differential
photometry with respect to a selection of local isolated and non-saturated
reference stars from the Pan-
STARRSfffhttp://outerspace.stsci.edu/display/PANSTARRS/. The presence of the
M82 galaxy makes the background highly non-uniform within the INTEGRAL error
circle. Magnitude upper limits have been computed for two regions (see Table
2): one (most conservative) corresponding to the centre of the M82 galaxy
(R.A. = ${09}^{\rm h}$ ${55}^{\rm m}$ ${53}^{\rm s}{\\!.75}$, Dec. =
${+69}^{\circ}$ ${40}^{\prime}$ ${53}^{\prime\prime}{\\!.9}$; J2000) and
another outside the core of the galaxy light (R.A. = ${09}^{\rm h}$ ${56}^{\rm
m}$ ${11}^{\rm s}{\\!.31}$, Dec. = ${+69}^{\circ}$ ${39}^{\prime}$
${09}^{\prime\prime}{\\!.6}$; J2000).

Table 2: Log of optical observations of GRB 231115A. Magnitudes are in the AB system, not corrected for Galactic extinction. Upper limits are given at 3-$\sigma$ confidence level for a source within (outside) the M82 galaxy.  UT observation | Exposure | T - T0 | Telescope | Magnitude | Filter
---|---|---|---|---|---
(start - stop) | (s) | (days) |  |  |
2023-11-15 20:40:03 - 20:54:15 |  $5\times 180$ s | 0.211 | Asiago | $>19.1(>20.1)$ | $g$
2023-11-15 20:54:59 - 21:12:11 |  $5\times 180$ s | 0.221 | Asiago | $>19.0(>20.1)$ | $r$
2023-11-15 21:12:45 - 21:27:45 |  $5\times 180$ s | 0.234 | Asiago | $>18.3(>19.8)$ | $i$
2023-11-15 21:38:00 - 22:16:30 |  $7\times 300$ s | 0.251 | Campo Imperatore | $>14.8(>18.9)$ | $z$
2023-11-15 22:38:25 - 21:31:30 |  $9\times 300$ s | 0.293 | Campo Imperatore | $>16.5(>19.8)$ | $i$
2023-11-15 23:33:32 - 21:38:00 |  $7\times 300$ s | 0.331 | Campo Imperatore | $>18.8(>20.7)$ | $g$
2023-11-15 22:35:32 - 22:57:23 |  $1\times 1400$ s | 0.299 | OHP | $>17.3(>21.6)$ | $r$
2023-11-15 23:04:02 - 23:14:02 |  $1\times 600$ s | 0.314 | OHP | $>18.0(>22.0)$ | $g$
2023-11-15 23:14:30 - 23:19:30 |  $1\times 300$ s | 0.320 | OHP | $>17.1(>20.1)$ | $i$
2023-11-16 03:25:04 - 03:38:16 |  $5\times 120$ s | 0.492 | TNG | $>20.0(>24.0)$ | $r$
2023-11-16 03:40:24 - 03:59:09 |  $7\times 120$ s | 0.503 | TNG | $>18.9(>23.8)$ | $i$
2023-11-16 04:00:28 - 04:19:17 |  $7\times 120$ s | 0.530 | TNG | $>18.7(>22.7)$ | $z$

Our upper limits in the optical bands (Table 2) exclude the majority of
cosmological short GRB afterglows. If occurred at the distance of M82, in a
position located inside (outside) the galaxy, a short GRB would have produced
an optical afterglow more than 12 (16) magnitudes brighter than our limits
(see Fig. 9). These upper limits are about 9 (13) magnitudes fainter than the
expected brightness of a kilonova like AT2017gfo occurring at the distance of
M82 inside (outside) the galaxy. Since kilonovae associated to short GRBs can
display significant differences in their luminosity, spectral properties, and
temporal evolution[12], we simulated a set of light curves with the POSSIS
code[13] using the parameter range reported in [14]. We found that even the
faintest event, placed at the distance of M82, is still about 5 (9) magnitudes
brighter than our upper limits.

Figure 8:  Optical light curves of kilonovae. The r-band light curves of
AT2017gro and of the faintest red+blue kilonova (simulated with the POSSIS
code) are shown with dashed and dashed-dotted lines, respectively, assuming
the M82 distance (3.6 Mpc). The magnitude 3-$\sigma$ upper limits obtained for
a position inside (outside) the M82 galaxy are shown as blue (red) arrows.

Figure 9:  Optical light curves of short GRB afterglows. The observed light
curves are shown in the left panel, while the right panel shows the same light
curves rescaled to the M82 distance (3.6 Mpc). The 3-$\sigma$ upper limits
obtained for a position inside (outside) the M82 galaxy are shown as blue
(red) arrows.

## References

## References

  * [1] Goldwurm, A., et al. The INTEGRAL/IBIS scientific data analysis. _A &A_, 411, L223 (2003).
  * [2] Arnaud, K. A. XSPEC: The First Ten Years. _Astronomical Data Analysis Software and Systems V_ , 101, 17 (1996).
  * [3] Mereghetti, S., Topinka, M., Rigoselli, M., & Götz, D. INTEGRAL Limits on Past High-energy Activity from FRB 20200120E in M81. _ApJ_ , 921, L3 (2021).
  * [4] HI4PI Collaboration, et al. HI4PI: A full-sky H I survey based on EBHIS and GASS. _A &A_, 594, A116 (2016).
  * [5] Strüder, L., et al. The European Photon Imaging Camera on XMM-Newton: The pn-CCD camera. _A &A_, 365, L18 (2001).
  * [6] Turner, M. J. L., et al. The European Photon Imaging Camera on XMM-Newton: The MOS cameras. _A &A_, 365, L27 (2001).
  * [7] Gabriel, C., et al. The XMM-Newton SAS - Distributed Development and Maintenance of a Large Science Analysis System: A Critical Analysis. _Astronomical Data Analysis Software and Systems (ADASS) XIII_ , 314, 759 (2004).
  * [8] Becker, A. HOTPANTS: High Order Transform of PSF ANd Template Subtraction. _Astrophysics Source Code Library_ , ascl:1504.004 (2015).
  * [9] Stetson, P. B. DAOPHOT: A Computer Program for Crowded-Field Stellar Photometry. _PASP_ , 99, 191 (1987).
  * [10] Karpov, S. STDPipe: Simple Transient Detection Pipeline. _arXiv e-prints_ , arXiv:2112.006 (2021).
  * [11] Bradley, L., Sipőcz, B. Robitaille, T., et al. astropy/photutils: 1.1.0, _Zenodo_ , 10.5281/zenodo.4624996 (2021).
  * [12] Rossi, A., et al. A comparison between short GRB afterglows and kilonova AT2017gfo: shedding light on kilonovae properties. _MNRAS_ , 493, 3379 (2020).
  * [13] Bulla, M. POSSIS: predicting spectra, light curves, and polarization for multidimensional models of supernovae and kilonovae. _MNRAS_ , 489, 5037 (2019).
  * [14] Ferro, M., et al. A search for the afterglows, kilonovae, and host galaxies of two short GRBs: GRB 211106A and GRB 211227A. _A &A_, 678, A142 (2023).

We thank the ESA Mission Scientists Jan-Uwe Ness and Norbert Schartel for
approving and quickly implementing the INTEGRAL and XMM-Newton ToO
observations. This work is based on observations with INTEGRAL and XMM-Newton
, ESA missions with instruments and science data centres funded by ESA member
states, and with the participation of the Russian Federation and the USA. It
is also based on observations made with the Italian Telescopio Nazionale
Galileo (TNG) operated on the island of La Palma by the Fundación Galileo
Galilei of the INAF (Istituto Nazionale di Astrofisica) at the Spanish
Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de
Canarias. This paper includes optical data taken with the Schmidt 67/92
telescope operated by INAF Osservatorio Astronomico di Padova (Mt. Ekar,
Asiago). This work received financial support from INAF through the Magnetars
Large Program Grant (PI S.Mereghetti) and from the GRAWITA Large Program Grant
(PI P. D’Avanzo). JCR, AB, SM and PU acknowledge financial support from ASI
under contract n. 2019-35-HH.0. FO acknowledges support from MIUR, PRIN 2020
(grant 2020KB33TP) “Multimessenger astronomy in the Einstein Telescope Era
(METE).

All authors reviewed the manuscript and contributed to the source
interpretation. SM coordinated the work and the interpretation of the results,
contributed to the analysis of the INTEGRAL and XMM-Newton data, and wrote
most of the manuscript. RS and EA contributed to write the main part of the
paper. DP and JCR carried out most of the INTEGRAL data analysis. DG, CF, EB,
LD and VS routinely contribute to the operation of the IBAS software and
participated to the near real time INTEGRAL analysis. PDA coordinated the
analysis of the optical data from Italian telescopes. MR analysed the XMM-
Newton data and contributed to the INTEGRAL spectral analysis. SC analysed the
Swift data. MT contributed to the software for the burst search in archival
data. DT, WT and CA coordinated the observation and the analysis of the
optical data taken at OHP. LT analysed the optical data taken with the Schmidt
67/92 telescope in Asiago under the Large Program ”Search and characterisation
of optical counterparts of GW triggers” (P.I. Tomasella). AR and EC triggered,
reduced and analysed the observations at the Asiago Schmidt telescope. RB and
MF provided the short GRB afterglows and kilonovae observed and simulated
optical light curves.

The authors declare no competing interests.

Correspondence and requests for materials should be addressed to Sandro
Mereghetti (email: sandro.mereghetti@inaf.it).

The data of the INTEGRAL, XMM-Newton and Swift satellites are publicly
available in the respective online archives
(https://www.isdc.unige.ch/integral/archive,
https://www.cosmos.esa.int/web/xmm-newton/xsa,
https://swift.gsfc.nasa.gov/archive/ )
. Optical data are available upon reasonable request.

The software used for the data analysis is public and can be retrieved at
https://www.cosmos.esa.int/web/xmm-newton/sas,
https://www.isdc.unige.ch/integral/analysis Software,
https://heasarc.gsfc.nasa.gov/xanadu/xspec/


# Efficient One Sided Kolmogorov Approximation

Liat Cohen1 Tal Grinshpoun2 Gera Weiss3
1University of Basel 2Ariel University 3Ben-Gurion University of the Negev

###### Abstract

We present an efficient algorithm that, given a discrete random variable $X$
and a number $m$, computes a random variable whose support is of size at most
$m$ and whose Kolmogorov distance from $X$ is minimal, also for the one-sided
Kolmogorov approximation. We present some variants of the algorithm, analyse
their correctness and computational complexity, and present a detailed
empirical evaluation that shows how they performs in practice. The main
application that we examine, which is our motivation for this work, is
estimation of the probability missing deadlines in series-parallel schedules.
Since exact computation of these probabilities is NP-hard, we propose to use
the algorithms described in this paper to obtain an approximation.

##  1 Introduction

Various approaches for approximation of probability distributions are studied
in the literature [17, 15, 20, 3, 16, 2]. These approaches vary in the types
random variables considered, how they are represented, and in the criteria
used for evaluation of the quality of the approximations. In this paper we
propose an approach for compressing the probability mass function of a random
variable $X$ such that the errors added to queries such as $Pr(X\leq t)$, for
any $t>0$, is minimal. In other words, we minimise the Kolmogorov distance
between the approximation and the original variable, see alternative
definition in Equation (1).

Our main motivation for this work is estimation of the probability for missing
deadlines, as described, e.g., in Cohen et al. [3, 2] and in [13].
Specifically, when $X$ represents the probability distribution of the time to
complete some complex schedule and we cannot afford to maintain the full table
of its probability mass function, we propose an algorithm for producing a
smaller table, whose size can be specified, such that probabilities for
missing deadlines are preserved as much as possible.

The main contribution of this paper is an efficient algorithm for computing
the best possible approximation of a given random variable with a random
variable whose size is not above a prescribed threshold, where the measures of
the quality of the approximation and of its size are as specified in the
following two paragraphs.

We measure the quality of an approximation scheme by the distance between
random variables and their approximations. Specifically, we use the Kolmogorov
distance which is commonly used for comparing random variables in statistical
practice and literature. Given two random variables $X$ and $X^{\prime}$ whose
cumulative distribution functions (cdf) are $F_{X}$ and $F_{X^{\prime}}$,
respectively, the Kolmogorov distance between $X$ and $X^{\prime}$ is
$d_{K}(X,X^{\prime})=\sup_{t}|F_{X}(t)-F_{X^{\prime}}(t)|$ (see, e.g., [7]).
We say that $X^{\prime}$ is a good approximation of $X$ if
$d_{K}(X,X^{\prime})$ is small. This distance is the basis for the often used
Kolmogorov-Smirnoff test for comparing a sample to a distribution or two
samples to each other.

The size of a random variable is measured by the size of its support, the set
of possible outcomes, $|X|{=}|\\{x\colon Pr(X{=}x)\neq 0\\}|$. When
probability mass functions are maintained as tables, as done in many
implementations of statistical software, the support size is proportional to
the memory needed to store the variable and to the complexity of the
computations that manipulate it. The exact notion of optimality of the
approximation targeted in this paper is:

######  Definition 1.

A random variable $X^{\prime}$ is an optimal $m$-approximation of a random
variable $X$ if $|X^{\prime}|\leq m$ and there is no random variable
$X^{\prime\prime}$ such that $|X^{\prime\prime}|\leq m$ and
$d_{K}(X,X^{\prime\prime})<d_{K}(X,X^{\prime})$.

In these terms, the main contribution of the paper is an efficient (linear
time and constant memory) algorithm that takes $X$ and $m$ as parameters and
constructs an optimal $m$-approximation of $X$.

The rest of the paper is organised as follows. In Section we describe how our
work relates to other algorithms and problems studied in the literature. In
Section we detail the proposed algorithm, analyse its properties, and prove
the main theorems. In Section we demonstrate how the proposed approach
performs on the problem of estimating the probability of missing deadlines in
series-parallel schedules on randomly generated random variables and compare
it to alternative approximation approaches from the literature. The paper is
concluded with a discussion and with ideas for future work in Section.

##  2 Related work

The most relevant work related to this paper is the papers on approximations
of random variables in the context of estimating deadlines [3, 2]. In these
papers, $X^{\prime}$ is defined to be a good approximation of $X$ if
$F_{X^{\prime}}(t)>F_{X}(t)$ for any $t$ and
$\sup_{t}F_{X^{\prime}}(t)-F_{X}(t)$ is small. Note that this measure is not a
proper distance measure because it is not symmetric. The motivation given in
these papers for using this type of approximation is for cases where
overestimation of the probability of missing a deadline is acceptable but
underestimation is not. We consider in this paper the same case-studies
examined by Cohen et al. and show how the algorithm proposed in this paper
performs relative to the algorithms proposed there when both over- and under-
estimations are allowed. As expected, the Kolmogorov distance between the
approximated and the original random variable is considerably smaller when
using the algorithm proposed in this paper.

In the technical level, the problem we study in this paper is similar to the
problem of approximating a set of 2-D points by a step function. The study of
this problem was motivated by query optimisation and histogram constructions
in database management systems [8, 9, 10, 11, 12, 6] and computational
geometry [4, 5]. There are, however, two technically significant differences
between the problem studied in the context of databases and the problem we
analyse in this paper. The first difference is that in the context of
approximation of random variables, the step function (which is the cumulative
distribution function in our context) must end with a value one, since we are
dealing with random variables which sums to one. The second difference is that
the first step is not counted because there is no need to put a value in the
support of the approximated random variable to generate this first step. These
cannot be addressed by adding a constant (two) to $m$ because the first step
is always present and because the requirement to end with the value one,
restricts the set of eligible step functions.

Another relevant prior work is the theory of Sparse Approximation (aka Sparse
Representation) that deals with sparse solutions for systems of linear
equations, as follows. Given a matrix $D\in\mathbb{R}^{n\times p}$ and a
vector $x\in\mathbb{R}^{n}$, the most studied sparse representation problem is
finding

$\min_{\alpha\in\mathbb{R}^{p}}\|\alpha\|_{0}\text{ subject to }x=D\alpha$

where $\|\alpha\|_{0}=|\\{i\in[p]:\alpha_{i}\neq 0\\}|$ is the $\ell_{0}$
pseudo-norm, counting the number of non-zero coordinates of $\alpha$. This
problem is known to be NP-hard with a reduction to NP-complete subset
selection problems. In these terms, using also the $\ell_{\infty}$ norm that
represents the maximal coordinate and the $\ell_{1}$ norm that represents the
sum of the coordinates, our problem can be phrased as:

$\min_{\alpha\in[0,\infty)^{p}}\|x-D\alpha\|_{\infty}\text{ subject to
}\|\alpha\|_{0}=m\text{ and }\|\alpha\|_{1}=1$ (1)

where $D$ is the lower unitriangular matrix, $x$ is related to $X$ such that
the $i$th coordinate of $x$ is $F_{X}(x_{i})$ where
$\operatorname{support}(X)=\\{x_{1}<\cdots<x_{n}\\}$ and $\alpha$ is related
to $X^{\prime}$ such that the $i$th coordinate of $\alpha$ is
$f_{X^{\prime}}(x_{i})$. The functions $F_{X}$ and $f_{X^{\prime}}$ represent,
respectively, the cumulative distribution function of $X$ and the mass
distribution function of $X^{\prime}$, i.e., the coordinates of $x$ are
positive and monotonically increasing and its last coordinate is one.

The presented work is also related to the research on binning in statistical
inference. Consider, for example, the problem of credit scoring [21] that
deals with separating good applicants from bad applicants where the
Kolmogorov–Smirnov statistic KS is a standard measure. The KS comparison is
often preceded by a procedure called binning where small values in the
probability mass function are moved to nearby values. There are many methods
for binning [14, 18, 1, 19]. In this context, our algorithm can be considered
as a binning strategy that provides optimality guarantees with respect to the
Kolmogorov distance.

Our study is also related to the work of Pavlikov and Uryasev [16], where a
procedure for producing a random variable $X^{\prime}$ that optimally
approximates a random variable $X$ is presented. Their approximation scheme,
achieved using linear programming, is designed for a different notion of
distance called CVaR. The contribution of the present work in this context is
that our method is direct, not using linear programming, thus allowing tighter
analysis of time and memory complexities. Also, our method is designed for
minimising the Kolmogorov distance that is more prevalent in applications. For
comparison, in Section we briefly discuss the performance of linear
programming approach similar to the one proposed in [16] for the Kolmogorov
distance and compare it our algorithm.

A problem very similar to ours is termed “order reduction” by Vidyasagar in
[20]. There, the author defines an information-theoretic based distance
between discrete random variables and studies the problem of finding a
variable whose support is of size $m$ and its distance from $X$ is as small as
possible (where $X$ and $m$ are given). The main difference between this and
the problem studied in this paper, is that Vidyasagar examines a different
notion of distance. Vidyasagar proves that computing the distance (that he
considers) between two probability distributions, and computing the optimal
reduced order approximation, are both NP-hard problems, because they can both
be reduced to nonstandard bin-packing problems. He then develops efficient
greedy approximation algorithms. In contrast, our study shows that there are
efficient solutions to these problems when the Kolmogorov distance is
considered.

######  Definition 2.

A random variable $X^{\prime}$ is an optimal one-sided $m$-approximation of a
random variable $X$ if $|X^{\prime}|\leq m$ and $F_{X}(t)<F_{X^{\prime}}(t)$
for all $t$ and there is no random variable $X^{\prime\prime}$ such that
$|X^{\prime\prime}|\leq m$ and $F_{X}(t)<F_{X^{\prime\prime}}(t)$ for all $t$
and $d_{K}(X,X^{\prime\prime})<d_{K}(X,X^{\prime})$.

1$f\leftarrow 1$

2 while  _$c_{f}\leq\varepsilon$_ do $f\leftarrow f+1$

3 $S\leftarrow\emptyset$, $s\leftarrow 0$, $e\leftarrow f$

4 while  _$e <n+1$_ do

5 while  _$c_{e+1}-c_{f}\leq\varepsilon\wedge e <n+1$_ do

6 $e\leftarrow e+1$

7 $S\leftarrow S\cup\\{(x_{f},c_{e})\\}$

8 $s\leftarrow c_{e}$, $f\leftarrow e\leftarrow e+1$

9

10if  _$s <1$_ then

11 $S\leftarrow S\cup\\{(x_{f},1-s)\\}$

12return _A r.v. $X^{\prime}$ as CDF such that $Pr(X^{\prime}{\leq}x)=c$ if
there is $c$ such that $(x,c)\in S$._

Algorithm 1 $dual(\\{(x_{i},c_{i})\\}_{i=1}^{n},\varepsilon)$

######  Example 3.

When $dual$ is invoked with the parameters
$X{=}\\{(1,0.3),(2,0.7),(3,0.9),(4,1)\\}$ and $\varepsilon{=}0.1$: Line 1
$f=1$. After an iteration of the main loop (line 4), $S=\\{(1,0.3)\\}$,
$s=0.3$, and $f=e=2$. After a second iteration, $S=\\{(1,0.3),(2,0.7)\\}$,
$s=0.7$, and $f=e=4$. At the end, $S=\\{(1,0.3),(2,0.7),(4,1)\\}$.

######  Proposition 4.

If $p=\\{(x_{i},c_{i})\\}_{i=1}^{n}$ is such that $c_{i}=Pr(X\leq x_{i})$ and
$\operatorname{support}(X)=\\{x_{1}<\cdots<x_{n}\\}$ then

$dual(\\{(x_{i},c_{i})\\}_{i=1}^{n},\varepsilon)\in\underset{X^{\prime}\in\bar{\mathcal{B}}(X,\varepsilon)}{\operatorname{arg\,min}}\,|X^{\prime}|$

where $\bar{\mathcal{B}}(X,\varepsilon)=\\{X^{\prime}\colon
d_{K}(X,X^{\prime})\leq\varepsilon\text{ and }F_{X}(x)<F_{X^{\prime}}(x)\text{
for all }x\\}$.

###### Proof.

The number of level sets of the CDF of $X^{\prime}$ is minimal: (1) The first
set is of maximal length; (2) By construction, an extension of any of the
other sets to the right will generate a random variable whose Kolmogorov
distance from $X$ is bigger than $\varepsilon$. Thus, there is no random
variable whose Kolmogorov distance from $X$ is smaller or equal to
$\varepsilon$ and its support is smaller than $|X^{\prime}|$. In addition, we
keep the invariant $f\leq e$ and since $X$ is given as CDF, $c_{f}\leq c_{e}$
consistent throughout the algorithm. ∎

######  Proposition 5.

$dual(\\{(x_{i},c_{i})\\}_{i=1}^{n},\varepsilon)$ runs in time $O(n)$, using
$O(n)$ memory.

###### Proof.

The algorithm describes a single pass over $\\{(x_{i},c_{i})\\}_{i=1}^{n}$.
Line 1 is easy to follow, and takes $O(n)$ in the worst case. Lines 3-7 also
describe a single pass since the counter $e$ is updated to $e+1$ at most $n$
times. All together we get run-time complexity of $O(n)$. We are constructing
the set $S$ which is of size $n$ in the worst case, therefore, memory
complexity is $O(n)$. ∎

1Let $\\{(x_{i},c_{i})\\}_{i=1}^{n}$ be such that $c_{i}=Pr(X\leq x_{i})$ and
$\operatorname{support}(X)=\\{x_{1}<\cdots<x_{n}\\}$.

2 $E\leftarrow\emptyset$

3 for  _$i\leftarrow 1$ to $n-1$_ do

4 for  _$j\leftarrow i+1$ to $n-1$_ do

5 if  _$i=1$_ then

6 $E\leftarrow E\cup\\{c_{j}\\}$

7 $E\leftarrow E\cup\\{c_{j}-c_{i}\\}$

8

9Let $e_{1}<\cdots<e_{n^{\prime}}$ be such that
$E=\\{e_{1},\dots,e_{n^{\prime}}\\}$

10

11for  _$i\leftarrow 1\textbf{ to }2$_ do

12 $l\leftarrow 1$, $r\leftarrow n^{\prime}$, $k\leftarrow 1$,
$k^{\prime}\leftarrow 0$

13 while  _$k\neq k^{\prime}$_ do

14 $k^{\prime}\leftarrow k$,$k\leftarrow\lceil(l+r)/2\rceil$

15 $m^{\prime}\leftarrow|dual(X,e_{k})|$

16 if  _$m^{\prime} <m$_ then $l\leftarrow k$

17 if  _$m^{\prime} >m$_ then $r\leftarrow k-1$

18 if  _$m^{\prime}=m$_ then $r\leftarrow k$

19

20 $m\leftarrow m^{\prime}$

21

22return _$dual(X,e_{k})$_

Algorithm 2 $binsApprox(X,m)$

######  Proposition 6.

$binsApprox(X,m)\in\underset{|X^{\prime}|\leq
m}{\operatorname{arg\,min}}\,d_{K}(X,X^{\prime})$

###### Proof.

Let $E=\\{Pr(x_{1}{<}X{\leq}x_{2})\colon
x_{1},x_{2}\in\operatorname{support}(X)\cup\\{-\infty,\infty\\},x_{1}<x_{2}\\}$
and let $E^{\prime}=\\{d_{K}(X,X^{\prime})\colon|X^{\prime}|\leq
m,F_{X}(x)\leq F_{X^{\prime}}(x)\text{ for all }x\\}.$ Since lines 10-18 are
describing a binary search to find the smallest $e\in E$ such that
$|dual(X,e)|\leq m$, we only need to prove that $\min E^{\prime}\in E$.

Let $\varepsilon=\min E^{\prime}$ and assume by contradiction that exists
$\varepsilon\notin E$. Let $X^{\prime}=dual(X,\varepsilon)$. By definition,
there is $x\in\operatorname{support}(X)$ in which
$F_{X^{\prime}}(x)-F_{X}(x)=\varepsilon$. Let $x_{-}=\max\\{x^{\prime}\colon
x^{\prime}\leq x,x^{\prime}\in\operatorname{support}(X^{\prime})\\}$. If
$x\in\operatorname{support}(X^{\prime})$ we have that $x_{-}=x$. Otherwise we
have that $Pr(x_{-}<X^{\prime}\leq x)=0$. In both cases we get that
$F_{X^{\prime}}(x_{-})=F_{X^{\prime}}(x)$. Plugging this in the previous
inequality gives us that $F_{X^{\prime}}(x_{-})-F_{X}(x)=\varepsilon$. Since
$X^{\prime}$ is constructed by the $dual$ algorithm,
$x_{-}\in\operatorname{support}(X)$ and $\exists
t\in\operatorname{support}(X)\colon F_{X^{\prime}}(x_{-})=F_{X}(t)$.
Therefore, $F_{X^{\prime}}(x_{-})-F_{X}(x)\in E$ in contradiction to
$\varepsilon\notin E$.

∎

######  Proposition 7.

The $binsApprox(X,m)$ algorithm runs in time $O(n^{2}\log(n))$, using
$O(n^{2})$ memory where $n=|X|$.

###### Proof.

In the first part of the algorithm, lines 2-8, we construct the set $E$ which
takes $O(n^{2})$ run-time. In the second part of the algorithm, line 9, we
sort the set $E$ which takes $O(n^{2}\log(n^{2}))=O(n^{2}\log(n))$ run-time.
The third part of the algorithm, lines 10-19, describes a binary sort over the
set $E$ where in each step of the sorting we run the $dual$ algorithm, which
takes $O(n\log(n^{2}))=O(n\log(n))$ run-time. We run this part twice then we
get $O(2n\log(n))$. All together, the run-time complexity is
$O(n^{2}+n^{2}\log(n)+2n\log(n))=O(n^{2}\log(n))$ and the memory complexity is
for storing the set $E$ which is $O(n^{2})$. ∎

Towards an improved algorithm let us introduce the matrix
$E=(e_{i,j})_{i,j=1}^{\infty}$ defined by:

$e_{i,j}=\begin{cases}c_{i}&\text{if }i\leq n\wedge j=n+1;\\\
c_{i}-c_{n+1-j}&\text{if }i<n\wedge j\leq n\wedge i{+}j\geq n;\\\ 0&\text{if
}i{+}j<n;\\\ 1&\text{otherwise}.\end{cases}$

Is is sorted??? I don’t think so..

Towards an improved algorithm let us introduce the matrix
$E=(e_{i,j})_{i,j=1}^{\infty}$ defined by:

$e_{i,j}=\begin{cases}1-c_{n+1-j}&\text{if }j\leq n\wedge i=n;\\\
c_{i}&\text{if }i\leq n\wedge j=n+1;\\\ (c_{i}-c_{n+1-j})/2&\text{if
}i<n\wedge j\leq n\wedge i{+}j\geq n;\\\ 0&\text{if }i{+}j<n;\\\
1&\text{otherwise}.\end{cases}$

Where $c_{1},\dots,c_{n}$ are as in the first line of Algorithm 2. It is easy
to see that the set of values in this matrix are the elements of the set $E$
in Algorithm 2. An additional useful fact is that $E$ is a sorted matrix:

######  Lemma 8.

If $i\leq i^{\prime}$ and $j\leq j^{\prime}$ then $e_{i,j}\leq
e_{i^{\prime},j^{\prime}}$.

###### Proof.

Since the $c_{i}$s are monotonically increasing, $\forall i,j$, $0\leq
c_{i}-c_{n+1-j}\leq c_{i}-c_{n+1-(j+1)}$ and $\forall i,j$, $0\leq
c_{i}-c_{n+1-j}\leq c_{i+1}-c_{n+1-j}$, moreover, 0 is the minimal value of
$E$. The elements in the last row and the last column in $E$ are keeping the
terms of sorted matrix. It is suffice to compare column $n$ with only column
$n+1$ and show that $c_{i}\geq(c_{i}-c_{n+1-n})$. After some manipulation we
get that $c_{1}\geq 0$ which is true. ∎

The fact that this matrix is sorted allows us to use the saddleback search
algorithm listed as Algorithm 3. The algorithm starts at the top right entry
of the matrix $(e_{i,j})_{i=1..n,j=1..(n+1)}$ and traverses it as follows. If
it hits an entry $e$ such that $|dual(X,e)|\leq m$ it goes left, otherwise it
goes down. This assures that the minimal $e$ such that $|dual(X,e)|\leq m$ is
visited after at most $n+1$ steps. The optimal random variable is found by
brute-force search over the visited entries.

1Let $\\{(x_{i},c_{i})\\}_{i=1}^{n}$ be such that $c_{i}=Pr(X\leq x_{i})$ and
$\operatorname{support}(X)=\\{x_{1}<\cdots<x_{n}\\}$.

2 $i\leftarrow 1$, $j\leftarrow n$+1, $S\leftarrow\emptyset$

3

4while  _$i <n\wedge j\geq 1$_ do

5 $m^{\prime}\leftarrow|dual(X,e_{i,j})|$

6

7 if  _$m^{\prime}\leq m$_ then $j\leftarrow j-1$

8 if  _$m^{\prime} >m$_ then $i\leftarrow i+1$

9 if  _$m^{\prime}\geq m$_ then $S\leftarrow S\cup(m^{\prime},e_{i,j})$

10

11$e\leftarrow\min\\{e\colon(m^{\prime},e)\in\underset{(m^{\prime},e)\in
S}{\operatorname{arg\,max}}\,m^{\prime}\\}$

12 return _$dual(X,e)$_

Algorithm 3 $sdlbkApprox(X,m)$

######  Proposition 9.

$sdlbkApprox(X,m){\in}\underset{|X^{\prime}|\leq
m}{\operatorname{arg\,min}}\,d_{K}(X,X^{\prime})$

###### Proof.

The algorithm traverses all the frontier between those $e$s that satisfy
$|dual(X,e)|\leq m$ and those that do not satisfy it. Since all the entries
that satisfy the condition are recorded in $S$ and considered in the brute-
force phase in line 7, the minimal satisfying $e$ is found. ∎

######  Proposition 10.

The $sdlbkApprox(X,m)$ algorithm runs in time $O(n^{2})$, using $O(n)$ memory
where $n=|X|$.

###### Proof.

At each step of the loop in lines 3-6, either $j$ decreases or $i$ increases,
thus the loop can be executed at most $2n$ times. Since we execute $dual$ once
in a loop round, the total time complexity is $O(n^{2})$. Storing visited
states in $S$ on line 6 requires $O(n)$ memory. ∎

The problem with the saddleback algorithm is that it needs to run $dual$ at
every step so it has a quadratic time complexity. Since we cannot find the
required entry of the matrix in less than $n$ steps, we can only reduce the
complexity by proposing an algorithm that does not execute $dual$ in all the
steps, only in $\log(n)$ of them. Such an algorithm, based on Section 2.1 of
[6], is listed as Algorithm 4. The algorithm maintains a set $S$ of sub-
matrices of
$(e_{i,j})_{i=1..2^{\lceil\log_{2}(n)\rceil},j=1..2^{\lceil\log_{2}(n+1)\rceil}}$.
At each round of its execution, each sub-matrix is split to four and then
about three quarters of the matrix are discarded. At the end, at most four
scalar matrices remain containing the index of the entry we seek. Note that
this algorithm runs on a matrix that, in the worst case, can be almost four
times bigger than the matrix traversed by the saddleback algorithm. This, of
course, does not affect the asymptotic complexity, but it may matter when
dealing with relatively small random variables.

1Let $\\{(x_{i},c_{i})\\}_{i=1}^{n}$ be such that $c_{i}=Pr(X\leq x_{i})$ and
$\operatorname{support}(X)=\\{x_{1}<\cdots<x_{n}\\}$.

2

3$S\leftarrow\\{((1,1),(2^{\lceil\log_{2}(n)\rceil},2^{\lceil\log_{2}(n+1)\rceil}))\\}$,
$S^{\prime}\leftarrow\emptyset$

4

5while  _$S\neq S^{\prime}$_ do

6 $S^{\prime}\leftarrow S$, $S\leftarrow\bigcup_{s\in S}split(s)$

7

8
$e^{-}{\leftarrow}median(\\{e_{i_{1},j_{1}}\colon((i_{1},j_{1}),(i_{2},j_{2}))\in
S\\})$

9
$e^{+}{\leftarrow}median(\\{e_{i_{2},j_{2}}\colon((i_{1},j_{1}),(i_{2},j_{2}))\in
S\\})$

10

11 for  _$e\in\\{e^{-},e^{+}\\}$_ do

12 $m^{\prime}\leftarrow|dual(X,e)|$

13

14 if  _$m^{\prime}\leq m$_ then

15 $S\leftarrow$ $S\setminus\\{((i_{1},j_{1}),(i_{2},j_{2}))\in S\colon
e_{i_{1},j_{1}}>e\\}$

16

17

18 else

19 $S\leftarrow$ $S\setminus\\{((i_{1},j_{1}),(i_{2},j_{2}))\in S\colon
e_{i_{2},j_{2}}\leq e\\}$

20

21

22

23

24$S^{\prime\prime}\leftarrow\\{(|dual(X,e_{i,j})|,e_{i,j})\colon((i,j),(i,j))\in
S\\}$

25

26$e\leftarrow\min\\{e\colon(m^{\prime},e)\in\underset{(m^{\prime},e)\in
S^{\prime\prime}}{\operatorname{arg\,max}}\,m^{\prime}\\}$

27 return _$dual(X,e)$_

Algorithm 4 $linApprox(X,m)$

[htbp]  $j^{-}\leftarrow\lfloor(j_{1}+j_{2})/2\rfloor$,
$j^{+}\leftarrow\lceil(j_{1}+j_{2})/2\rceil$

2 $i^{-}\leftarrow\lfloor(i_{1}+i_{2})/2\rfloor$,
$i^{+}\leftarrow\lceil(i_{1}+i_{2})/2\rceil$

return _ $\\{((i_{1},j_{1}),(i^{-},j^{-})),((i_{1},j^{+}),(i^{-},j_{2}))$,
$((i^{+},j_{1}),(i_{2},j^{-})),((i^{+},j^{+}),(i_{2},j_{2}))\\}$_
$split(((i_{1},j_{1}),(i_{2},j_{2})))$()

######  Theorem 11.

$linApprox(X,m){\in}\underset{|X^{\prime}|\leq
m}{\operatorname{arg\,min}}\,d_{K}(X,X^{\prime})$

###### Proof.

In line 10 we only discard sub-matrices whose minimal entry (at the top left)
is larger than an entry that we prefer. In line 12 we only discard sub-
matrices whose maximal entry (at the bottom right) is smaller or equal to an
entry that does not meet our condition. ∎

######  Theorem 12.

The $linApprox(X,m)$ algorithm runs in time $O(n\log(n))$, using $O(1)$ memory
where $n=|X|$.

###### Proof.

The dimension of each matrix is halved in each round, thus the loop is
executed $O(\log(n))$ rounds. Since $dual$ is called one in each round and a
finite number of times at the end, the time complexity is $O(n\log(n))$. ∎

## References

  * [1] Christine Bolton et al.  Logistic regression and its application in credit scoring.  PhD thesis, University of Pretoria, 2010.
  * [2] L Cohen, T Grinshpoun, and G Weiss.  Optimal approximation of random variables for estimating the probability of meeting a plan deadline.  In AAAI, pages 6327–6334, 2018.
  * [3] Liat Cohen, Solomon Eyal Shimony, and Gera Weiss.  Estimating the probability of meeting a deadline in hierarchical plans.  In IJCAI, pages 1551–1557, 2015.
  * [4] José Miguel Díaz-Bánez and Juan A Mesa.  Fitting rectilinear polygonal curves to a set of points in the plane.  European Journal of Operational Research, 130(1):214–222, 2001\.
  * [5] Hervé Fournier and Antoine Vigneron.  Fitting a step function to a point set.  In European Symposium on Algorithms, pages 442–453. Springer, 2008\.
  * [6] Hervé Fournier and Antoine Vigneron.  Fitting a step function to a point set.  Algorithmica, 60(1):95–109, May 2011.
  * [7] Jean Dickinson Gibbons and Subhabrata Chakraborti.  Nonparametric statistical inference.  In International encyclopedia of statistical science, pages 977–979. Springer, 2011.
  * [8] Subharup Guha.  Posterior simulation in the generalized linear mixed model with semiparametric random effects.  J. Comput. Graph. Statist., 17(2):410–425, 2008.
  * [9] Sudipto Guha and Kyuseok Shim.  A note on linear time algorithms for maximum error histograms.  IEEE Transactions on Knowledge and Data Engineering, 19(7):993–997, 2007.
  * [10] Sudipto Guha, Kyuseok Shim, and Jungchul Woo.  Rehist: Relative error histogram construction algorithms.  In VLDB, pages 300–311, 2004.
  * [11] Hosagrahar Visvesvaraya Jagadish, Nick Koudas, S Muthukrishnan, Viswanath Poosala, Kenneth C Sevcik, and Torsten Suel.  Optimal histograms with quality guarantees.  In VLDB, volume 98, pages 24–27, 1998.
  * [12] Panagiotis Karras, Dimitris Sacharidis, and Nikos Mamoulis.  Exploiting duality in summarization with deterministic guarantees.  In SIGKDD, pages 380–389, 2007.
  * [13] M. Kashef and N. Moayeri.  Random-deadline missing probability analysis for wireless communications in industrial environments.  In WFCS, pages 1–8, 2018.
  * [14] Elizabeth Mays.  Handbook of credit scoring.  Global Professional Publishi, 2001.
  * [15] Allen C. Miller and Thomas R. Rice.  Discrete approximations of probability distributions.  Management Science, 29(3):352–362, 1983.
  * [16] Konstantin Pavlikov and Stan Uryasev.  CVaR distance between univariate probability distributions and approximation problems.  Technical Report 2015-6, University of Florida, 2016.
  * [17] A. N. Pettitt and M. A. Stephens.  The kolmogorov-smirnov goodness-of-fit statistic with discrete and grouped data.  Technometrics, 19(2):205–210, 1977.
  * [18] Mamdouh Refaat.  Credit Risk Scorecard: Development and Implementation Using SAS.  Lulu. com, 2011.
  * [19] Naeem Siddiqi.  Credit risk scorecards: developing and implementing intelligent credit scoring, volume 3.  John Wiley & Sons, 2012.
  * [20] Mathukumalli Vidyasagar.  A metric between probability distributions on finite sets of different cardinalities and applications to order reduction.  IEEE Transactions on Automatic Control, 57(10):2464–2477, 2012\.
  * [21] Guoping Zeng.  A comparison study of computational methods of kolmogorov–smirnov statistic in credit scoring.  Communications in Statistics-Simulation and Computation, 46(10):7744–7760, 2017.