Skip to content

API reference

optimum_interval

Optimum-interval upper limits in the presence of unknown background.

A small, dependency-light implementation of the method of S. Yellin, "Finding an Upper Limit in the Presence of Unknown Background", Phys. Rev. D66 (2002) 032005 (arXiv:physics/0203002).

The method sets a frequentist upper limit on a signal normalization (the Poisson mean mu of expected events) using only the known signal shape, without a model of the background and without binning. See EXPLANATION.md for a physicist-oriented walk-through and derivation.

Quick start

import numpy as np from optimum_interval import OptimumIntervalTable rng = np.random.default_rng(0) table = OptimumIntervalTable(rng=rng) events = np.sort(rng.random(8)) # 8 events, already in cumulant space mu_limit = table.upper_limit(events, confidence=0.9, n=2000) # doctest: +SKIP

Public API

k_largest_intervals / cumulant_points Pure interval geometry (:mod:optimum_interval.intervals). c0 / x0 / max_gap_upper_limit / poisson_upper_limit Analytic statistics: Yellin Eq. 2 and the counting limits (:mod:optimum_interval.analytic). OptimumIntervalTable Monte-Carlo calibration tables and the upper-limit solver (:mod:optimum_interval.montecarlo). spectrum_cdf_from_pdf / spectrum_cdf_from_samples Build a normalized spectrum_cdf (:mod:optimum_interval.spectra). scan_extremeness / excluded_interval / spectrum_from_rate Parameter scans when the spectrum shape depends on the parameter being limited (:mod:optimum_interval.scanning). ComparisonEngine Fast per-experiment limits for method comparisons (:mod:optimum_interval.comparison).

ComparisonEngine

Precomputed-grid C_max and p_max upper-limit solver.

Parameters:

Name Type Description Default
mu_grid array_like

Grid of mu values over which the calibration is tabulated. Must span the range in which the limits are expected to fall.

 required
n_cal int

Monte-Carlo trials per grid point.

 40000
rng Generator

Seed for reproducibility.

 None
confidence float

Confidence level (default 0.9).

 0.9
Source code in optimum_interval/comparison.py 
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
class ComparisonEngine:
    """Precomputed-grid ``C_max`` and ``p_max`` upper-limit solver.

    Parameters
    ----------
    mu_grid : array_like
        Grid of ``mu`` values over which the calibration is tabulated.  Must
        span the range in which the limits are expected to fall.
    n_cal : int, optional
        Monte-Carlo trials per grid point.
    rng : numpy.random.Generator, optional
        Seed for reproducibility.
    confidence : float, optional
        Confidence level (default 0.9).
    """

    def __init__(
        self,
        mu_grid,
        n_cal: int = 40000,
        rng: np.random.Generator | None = None,
        confidence: float = 0.9,
    ) -> None:
        self.mu_grid = np.sort(np.asarray(mu_grid, dtype=float))
        self.n_cal = int(n_cal)
        self.confidence = float(confidence)
        self.rng = np.random.default_rng() if rng is None else rng

        # Per grid mu: sorted k-largest reference sizes, and the sorted
        # C_max / p_max trial distributions used as calibration.
        self._ref: list[dict[int, np.ndarray]] = []
        self._opt: list[np.ndarray] = []
        self._pmax: list[np.ndarray] = []
        self._build()

    # ------------------------------------------------------------------ #
    def _build(self) -> None:
        n = self.n_cal
        for mu in self.mu_grid:
            counts = self.rng.poisson(mu, size=n)
            per_trial: list[dict[int, float]] = []
            sizes_by_k: dict[int, list[float]] = defaultdict(list)
            trials_by_k: dict[int, list[int]] = defaultdict(list)
            for t, count in enumerate(counts):
                pts = np.concatenate(([0.0], np.sort(self.rng.random(count)), [1.0]))
                sizes = k_largest_intervals(pts)
                per_trial.append(sizes)
                for k, size in sizes.items():
                    sizes_by_k[k].append(size)
                    trials_by_k[k].append(t)

            sorted_ref = {k: np.sort(v) for k, v in sizes_by_k.items()}

            # C_max per trial (empirical CDF of each k-largest, max over k).
            opt = np.zeros(n)
            for k, values in sizes_by_k.items():
                extremeness = np.searchsorted(sorted_ref[k], values, side="left") / n
                np.maximum.at(opt, np.asarray(trials_by_k[k]), extremeness)

            # p_max per trial: max over k of Poisson P(>k | mu * size).
            pmax = np.empty(n)
            for t, sizes in enumerate(per_trial):
                ks = np.fromiter(sizes.keys(), dtype=int)
                xs = mu * np.fromiter(sizes.values(), dtype=float)
                pmax[t] = poisson.sf(ks, xs).max()

            self._ref.append(sorted_ref)
            self._opt.append(np.sort(opt))
            self._pmax.append(np.sort(pmax))

    # ------------------------------------------------------------------ #
    @staticmethod
    def _sizes_array(events, spectrum_cdf) -> np.ndarray:
        """Observed k-largest sizes as a dense array indexed by k."""
        sizes = k_largest_intervals(cumulant_points(events, spectrum_cdf))
        return np.array([sizes[k] for k in range(len(sizes))])

    def _crossing(self, g: np.ndarray) -> float:
        """First ``mu`` on the grid where increasing curve ``g`` reaches confidence.

        Linear interpolation between grid points.  Clamps (with no silent
        surprise -- the grid should be chosen to contain the limit) to the grid
        ends if the crossing lies outside.
        """
        c = self.confidence
        grid = self.mu_grid
        above = g >= c
        if not above.any():
            warnings.warn(
                f"upper limit exceeds mu_grid max {grid[-1]:g}; clamping. Widen mu_grid.",
                stacklevel=3,
            )
            return float(grid[-1])
        i = int(np.argmax(above))
        if i == 0:
            warnings.warn(
                f"upper limit is at/below mu_grid min {grid[0]:g}; clamping. "
                "Lower mu_grid.",
                stacklevel=3,
            )
            return float(grid[0])
        y0, y1 = g[i - 1], g[i]
        if y1 == y0:
            return float(grid[i])
        return float(grid[i - 1] + (c - y0) * (grid[i] - grid[i - 1]) / (y1 - y0))

    # ------------------------------------------------------------------ #
    def cmax_upper_limit(self, events, spectrum_cdf=None) -> float:
        """Optimum-interval (``C_max``) upper limit on ``mu`` for one experiment."""
        sizes = self._sizes_array(events, spectrum_cdf)
        g = np.empty(self.mu_grid.size)
        for j in range(self.mu_grid.size):
            ref = self._ref[j]
            cmax = 0.0
            for k in range(min(sizes.size, max(ref) + 1 if ref else 0)):
                arr = ref.get(k)
                if arr is not None:
                    e = np.searchsorted(arr, sizes[k], side="left") / self.n_cal
                    if e > cmax:
                        cmax = e
            g[j] = np.searchsorted(self._opt[j], cmax, side="left") / self.n_cal
        return self._crossing(g)

    def pmax_upper_limit(self, events, spectrum_cdf=None) -> float:
        """``p_max`` (Poisson-probability) upper limit on ``mu`` for one experiment."""
        sizes = self._sizes_array(events, spectrum_cdf)
        ks = np.arange(sizes.size)
        g = np.empty(self.mu_grid.size)
        for j, mu in enumerate(self.mu_grid):
            pmax = poisson.sf(ks, mu * sizes).max()
            g[j] = np.searchsorted(self._pmax[j], pmax, side="left") / self.n_cal
        return self._crossing(g)

cmax_upper_limit(events, spectrum_cdf=None)

Optimum-interval (C_max) upper limit on mu for one experiment.

Source code in optimum_interval/comparison.py 
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
def cmax_upper_limit(self, events, spectrum_cdf=None) -> float:
    """Optimum-interval (``C_max``) upper limit on ``mu`` for one experiment."""
    sizes = self._sizes_array(events, spectrum_cdf)
    g = np.empty(self.mu_grid.size)
    for j in range(self.mu_grid.size):
        ref = self._ref[j]
        cmax = 0.0
        for k in range(min(sizes.size, max(ref) + 1 if ref else 0)):
            arr = ref.get(k)
            if arr is not None:
                e = np.searchsorted(arr, sizes[k], side="left") / self.n_cal
                if e > cmax:
                    cmax = e
        g[j] = np.searchsorted(self._opt[j], cmax, side="left") / self.n_cal
    return self._crossing(g)

pmax_upper_limit(events, spectrum_cdf=None)

p_max (Poisson-probability) upper limit on mu for one experiment.

Source code in optimum_interval/comparison.py 
150
151
152
153
154
155
156
157
158
def pmax_upper_limit(self, events, spectrum_cdf=None) -> float:
    """``p_max`` (Poisson-probability) upper limit on ``mu`` for one experiment."""
    sizes = self._sizes_array(events, spectrum_cdf)
    ks = np.arange(sizes.size)
    g = np.empty(self.mu_grid.size)
    for j, mu in enumerate(self.mu_grid):
        pmax = poisson.sf(ks, mu * sizes).max()
        g[j] = np.searchsorted(self._pmax[j], pmax, side="left") / self.n_cal
    return self._crossing(g)

OptimumIntervalTable

Monte-Carlo interval-size tables plus the upper-limit solver.

Parameters:

Name Type Description Default
rng Generator

Random generator. Pass a seeded np.random.default_rng(seed) for reproducible tables. Defaults to a fresh unseeded generator.

 None

Attributes:

Name Type Description
itv_sizes dict[float, dict[int, ndarray]]

itv_sizes[mu][k] is the array of k-largest interval sizes over all trials that contained at least k events.
opt_itvs dict[float, ndarray]

opt_itvs[mu] is the per-trial C_max statistic (the calibration distribution whose quantile gives bar_c_max).
n_trials dict[float, int]

Number of Monte-Carlo trials generated for each mu.
Source code in optimum_interval/montecarlo.py 
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
class OptimumIntervalTable:
    """Monte-Carlo interval-size tables plus the upper-limit solver.

    Parameters
    ----------
    rng : numpy.random.Generator, optional
        Random generator.  Pass a seeded ``np.random.default_rng(seed)`` for
        reproducible tables.  Defaults to a fresh unseeded generator.

    Attributes
    ----------
    itv_sizes : dict[float, dict[int, numpy.ndarray]]
        ``itv_sizes[mu][k]`` is the array of k-largest interval sizes over all
        trials that contained at least ``k`` events.
    opt_itvs : dict[float, numpy.ndarray]
        ``opt_itvs[mu]`` is the per-trial C_max statistic (the calibration
        distribution whose quantile gives ``bar_c_max``).
    n_trials : dict[float, int]
        Number of Monte-Carlo trials generated for each ``mu``.
    """

    def __init__(self, rng: np.random.Generator | None = None) -> None:
        self.itv_sizes: dict[float, dict[int, np.ndarray]] = {}
        self.opt_itvs: dict[float, np.ndarray] = {}
        self.n_trials: dict[float, int] = {}
        self.rng: np.random.Generator = np.random.default_rng() if rng is None else rng

    # ------------------------------------------------------------------ #
    # Monte-Carlo table construction
    # ------------------------------------------------------------------ #
    def generate_trials(self, mu: float, n: int) -> list[np.ndarray]:
        """Draw ``n`` background-free trial experiments at Poisson mean ``mu``.

        Each trial is a sorted point list on ``[0, 1]``: ``Poisson(mu)`` interior
        uniforms plus the range endpoints 0 and 1.  This mirrors exactly the
        cumulant-space representation of real data (:func:`cumulant_points`).
        """
        counts = self.rng.poisson(mu, size=n)
        trials = []
        for count in counts:
            interior = np.sort(self.rng.random(count))
            trials.append(np.concatenate(([0.0], interior, [1.0])))
        return trials

    def generate(self, mu: float, n: int) -> None:
        """Ensure tables for ``mu`` are built from at least ``n`` trials.

        Populates :attr:`itv_sizes`, :attr:`opt_itvs` and :attr:`n_trials`.
        Cached: returns immediately if ``mu`` already has ``>= n`` trials.
        """
        if self.n_trials.get(mu, 0) >= n:
            return

        log.info("generating optimum-interval table for mu=%s (n=%d)", mu, n)
        trials = self.generate_trials(mu, n)

        # k-largest sizes per trial, grouped by k with the originating trial
        # index kept so we can scatter the extremeness back per trial.
        sizes_by_k: dict[int, list[float]] = defaultdict(list)
        trials_by_k: dict[int, list[int]] = defaultdict(list)
        for t, trial in enumerate(trials):
            for k, size in k_largest_intervals(trial).items():
                sizes_by_k[k].append(size)
                trials_by_k[k].append(t)

        itv_sizes = {k: np.asarray(v) for k, v in sizes_by_k.items()}

        # C_max per trial = max over k of the empirical CDF (extremeness) of
        # that trial's k-largest size.  extremeness uses strict "<" and divides
        # by the total trial count n (Yellin's denominator convention: a trial
        # with too few events to have a k-largest interval does not count as
        # "smaller").  Vectorized per k via searchsorted, then scattered.
        opt = np.zeros(n)
        for k, obs in itv_sizes.items():
            reference = np.sort(obs)
            extremeness = np.searchsorted(reference, obs, side="left") / n
            np.maximum.at(opt, np.asarray(trials_by_k[k]), extremeness)

        self.itv_sizes[mu] = itv_sizes
        self.opt_itvs[mu] = opt
        self.n_trials[mu] = n

    # ------------------------------------------------------------------ #
    # Statistics that query the tables
    # ------------------------------------------------------------------ #
    def extremeness_of_interval(self, x: float, k: int, mu: float) -> float:
        """Fraction of MC k-largest intervals at ``mu`` smaller than ``x``.

        This is the code's empirical stand-in for :math:`C_k(x, \\mu)`.  Returns
        0 if no trial produced a k-largest interval (missing ``k``).
        """
        reference = self.itv_sizes[mu].get(k)
        if reference is None:
            return 0.0
        return float(np.count_nonzero(reference < x) / self.n_trials[mu])

    def optimum_interval_statistic(
        self, events: np.ndarray, mu: float, spectrum_cdf: SpectrumCdf | None = None
    ) -> float:
        """C_max of a run: the most extreme k-largest interval.

        ``events`` are raw positions; they are transformed to cumulant space and
        the range endpoints 0 and 1 are added (:func:`cumulant_points`) so the
        calculation matches the Monte-Carlo calibration.
        """
        points = cumulant_points(events, spectrum_cdf)
        sizes = k_largest_intervals(points)  # already in cumulant space
        return max(self.extremeness_of_interval(size, k, mu) for k, size in sizes.items())

    def extremeness_of_opt_itv_stat(self, stat: float, mu: float) -> float:
        """Fraction of MC C_max values at ``mu`` smaller than ``stat``.

        The confidence quantile of this distribution is ``bar_c_max``.
        """
        opt = self.opt_itvs[mu]
        return float(np.count_nonzero(opt < stat) / self.n_trials[mu])

    def bar_c_max(self, confidence: float, mu: float) -> float:
        """Threshold :math:`\\bar C_\\mathrm{max}(C, \\mu)` (Yellin Fig. 2).

        The value of C_max reached with probability ``confidence`` under the
        background-free hypothesis: the ``confidence`` quantile of
        :attr:`opt_itvs`.  ``generate`` must have been called for ``mu`` first.
        """
        return float(np.quantile(self.opt_itvs[mu], confidence))

    # ------------------------------------------------------------------ #
    # Upper limit
    # ------------------------------------------------------------------ #
    def upper_limit(
        self,
        events: np.ndarray,
        confidence: float = 0.9,
        spectrum_cdf: SpectrumCdf | None = None,
        n: int = 1000,
        *,
        mu_scan_start: float = 1.0,
        mu_scan_stop: float | None = None,
        mu_scan_step: float = 1.0,
    ) -> float:
        """``confidence``-level optimum-interval upper limit on ``mu``.

        Finds ``mu`` such that ``extremeness_of_opt_itv_stat(C_max(events, mu),
        mu) == confidence`` -- i.e. where the observed C_max equals
        :math:`\\bar C_\\mathrm{max}(\\text{confidence}, \\mu)`.

        The excess ``extremeness - confidence`` increases with ``mu`` (a larger
        proposed signal makes the observed emptiness look more anomalous).  It is
        evaluated on a fixed ``mu`` grid whose tables are cached on this instance,
        and the zero crossing is interpolated between the bracketing grid points,
        so the result is deterministic given the generator seed.

        Parameters
        ----------
        events : array_like
            Observed event positions (raw; transformed via ``spectrum_cdf``).
        confidence : float, optional
            Confidence level, e.g. ``0.9``.
        spectrum_cdf : callable, optional
            Signal CDF mapping ``events`` to cumulants.  Identity by default.
        n : int, optional
            Monte-Carlo trials per grid ``mu``.  Raise ``n`` for a less noisy limit.
        mu_scan_start, mu_scan_stop, mu_scan_step : float, optional
            The ``mu`` grid.  ``mu_scan_stop`` defaults to roughly twice the event
            count (never below ``mu_scan_start + 10``).

        Returns
        -------
        float
            The upper limit on ``mu``.

        Raises
        ------
        RuntimeError
            If no exclusion is bracketed within the scan range; widen
            ``mu_scan_start`` / ``mu_scan_stop`` as the message suggests.
        """
        events = np.asarray(events, dtype=float)

        def excess(mu: float) -> float:
            self.generate(mu, n)  # cached per mu -> deterministic given the seed
            stat = self.optimum_interval_statistic(events, mu, spectrum_cdf)
            return self.extremeness_of_opt_itv_stat(stat, mu) - confidence

        def interpolate(mu_lo, g_lo, mu_hi, g_hi):
            if g_hi == g_lo:
                return float(mu_hi)
            return float(mu_lo + (0.0 - g_lo) * (mu_hi - mu_lo) / (g_hi - g_lo))

        if mu_scan_stop is None:
            mu_scan_stop = max(2.0 * events.size + 5.0, mu_scan_start + 10.0)

        prev_mu, prev_g = None, None
        for mu in np.arange(mu_scan_start, mu_scan_stop + mu_scan_step, mu_scan_step):
            mu = float(mu)
            g = excess(mu)
            if g >= 0.0:
                if prev_mu is None:
                    # Limit lies below the first probed mu; step down to bracket.
                    lo, g_lo = mu, g
                    for _ in range(20):
                        lo *= 0.5
                        g_lo = excess(lo)
                        if lo < 1e-3 or g_lo < 0.0:
                            break
                    if lo < 1e-3 or g_lo >= 0.0:
                        raise RuntimeError(
                            f"upper limit is below mu_scan_start={mu_scan_start}; "
                            f"lower it."
                        )
                    prev_mu, prev_g = lo, g_lo
                log.info("bracketed limit in [%s, %s]", prev_mu, mu)
                return interpolate(prev_mu, prev_g, mu, g)
            prev_mu, prev_g = mu, g

        raise RuntimeError(
            f"no {confidence:.0%} exclusion found up to mu={mu_scan_stop}; "
            f"increase mu_scan_stop."
        )

    # ------------------------------------------------------------------ #
    # Persistence
    # ------------------------------------------------------------------ #
    def save(self, path: str | Path) -> None:
        """Pickle the tables to ``path`` (single dict, one payload)."""
        payload = {
            "itv_sizes": self.itv_sizes,
            "opt_itvs": self.opt_itvs,
            "n_trials": self.n_trials,
        }
        with open(path, "wb") as f:
            pickle.dump(payload, f)

    @classmethod
    def load(
        cls, path: str | Path, rng: np.random.Generator | None = None
    ) -> OptimumIntervalTable:
        """Load tables previously written by :meth:`save`."""
        table = cls(rng=rng)
        with open(path, "rb") as f:
            payload = pickle.load(f)
        table.itv_sizes = payload["itv_sizes"]
        table.opt_itvs = payload["opt_itvs"]
        table.n_trials = payload["n_trials"]
        return table

generate_trials(mu, n)

Draw n background-free trial experiments at Poisson mean mu.

Each trial is a sorted point list on [0, 1]: Poisson(mu) interior uniforms plus the range endpoints 0 and 1. This mirrors exactly the cumulant-space representation of real data (:func:cumulant_points).

Source code in optimum_interval/montecarlo.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
def generate_trials(self, mu: float, n: int) -> list[np.ndarray]:
    """Draw ``n`` background-free trial experiments at Poisson mean ``mu``.

    Each trial is a sorted point list on ``[0, 1]``: ``Poisson(mu)`` interior
    uniforms plus the range endpoints 0 and 1.  This mirrors exactly the
    cumulant-space representation of real data (:func:`cumulant_points`).
    """
    counts = self.rng.poisson(mu, size=n)
    trials = []
    for count in counts:
        interior = np.sort(self.rng.random(count))
        trials.append(np.concatenate(([0.0], interior, [1.0])))
    return trials

generate(mu, n)

Ensure tables for mu are built from at least n trials.

Populates :attr:itv_sizes, :attr:opt_itvs and :attr:n_trials. Cached: returns immediately if mu already has >= n trials.

Source code in optimum_interval/montecarlo.py 
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def generate(self, mu: float, n: int) -> None:
    """Ensure tables for ``mu`` are built from at least ``n`` trials.

    Populates :attr:`itv_sizes`, :attr:`opt_itvs` and :attr:`n_trials`.
    Cached: returns immediately if ``mu`` already has ``>= n`` trials.
    """
    if self.n_trials.get(mu, 0) >= n:
        return

    log.info("generating optimum-interval table for mu=%s (n=%d)", mu, n)
    trials = self.generate_trials(mu, n)

    # k-largest sizes per trial, grouped by k with the originating trial
    # index kept so we can scatter the extremeness back per trial.
    sizes_by_k: dict[int, list[float]] = defaultdict(list)
    trials_by_k: dict[int, list[int]] = defaultdict(list)
    for t, trial in enumerate(trials):
        for k, size in k_largest_intervals(trial).items():
            sizes_by_k[k].append(size)
            trials_by_k[k].append(t)

    itv_sizes = {k: np.asarray(v) for k, v in sizes_by_k.items()}

    # C_max per trial = max over k of the empirical CDF (extremeness) of
    # that trial's k-largest size.  extremeness uses strict "<" and divides
    # by the total trial count n (Yellin's denominator convention: a trial
    # with too few events to have a k-largest interval does not count as
    # "smaller").  Vectorized per k via searchsorted, then scattered.
    opt = np.zeros(n)
    for k, obs in itv_sizes.items():
        reference = np.sort(obs)
        extremeness = np.searchsorted(reference, obs, side="left") / n
        np.maximum.at(opt, np.asarray(trials_by_k[k]), extremeness)

    self.itv_sizes[mu] = itv_sizes
    self.opt_itvs[mu] = opt
    self.n_trials[mu] = n

extremeness_of_interval(x, k, mu)

Fraction of MC k-largest intervals at mu smaller than x.

This is the code's empirical stand-in for :math:C_k(x, \mu). Returns 0 if no trial produced a k-largest interval (missing k).

Source code in optimum_interval/montecarlo.py 
110
111
112
113
114
115
116
117
118
119
def extremeness_of_interval(self, x: float, k: int, mu: float) -> float:
    """Fraction of MC k-largest intervals at ``mu`` smaller than ``x``.

    This is the code's empirical stand-in for :math:`C_k(x, \\mu)`.  Returns
    0 if no trial produced a k-largest interval (missing ``k``).
    """
    reference = self.itv_sizes[mu].get(k)
    if reference is None:
        return 0.0
    return float(np.count_nonzero(reference < x) / self.n_trials[mu])

optimum_interval_statistic(events, mu, spectrum_cdf=None)

C_max of a run: the most extreme k-largest interval.

events are raw positions; they are transformed to cumulant space and the range endpoints 0 and 1 are added (:func:cumulant_points) so the calculation matches the Monte-Carlo calibration.

Source code in optimum_interval/montecarlo.py 
121
122
123
124
125
126
127
128
129
130
131
132
def optimum_interval_statistic(
    self, events: np.ndarray, mu: float, spectrum_cdf: SpectrumCdf | None = None
) -> float:
    """C_max of a run: the most extreme k-largest interval.

    ``events`` are raw positions; they are transformed to cumulant space and
    the range endpoints 0 and 1 are added (:func:`cumulant_points`) so the
    calculation matches the Monte-Carlo calibration.
    """
    points = cumulant_points(events, spectrum_cdf)
    sizes = k_largest_intervals(points)  # already in cumulant space
    return max(self.extremeness_of_interval(size, k, mu) for k, size in sizes.items())

extremeness_of_opt_itv_stat(stat, mu)

Fraction of MC C_max values at mu smaller than stat.

The confidence quantile of this distribution is bar_c_max.

Source code in optimum_interval/montecarlo.py 
134
135
136
137
138
139
140
def extremeness_of_opt_itv_stat(self, stat: float, mu: float) -> float:
    """Fraction of MC C_max values at ``mu`` smaller than ``stat``.

    The confidence quantile of this distribution is ``bar_c_max``.
    """
    opt = self.opt_itvs[mu]
    return float(np.count_nonzero(opt < stat) / self.n_trials[mu])

bar_c_max(confidence, mu)

Threshold :math:\bar C_\mathrm{max}(C, \mu) (Yellin Fig. 2).

The value of C_max reached with probability confidence under the background-free hypothesis: the confidence quantile of :attr:opt_itvs. generate must have been called for mu first.

Source code in optimum_interval/montecarlo.py 
142
143
144
145
146
147
148
149
def bar_c_max(self, confidence: float, mu: float) -> float:
    """Threshold :math:`\\bar C_\\mathrm{max}(C, \\mu)` (Yellin Fig. 2).

    The value of C_max reached with probability ``confidence`` under the
    background-free hypothesis: the ``confidence`` quantile of
    :attr:`opt_itvs`.  ``generate`` must have been called for ``mu`` first.
    """
    return float(np.quantile(self.opt_itvs[mu], confidence))

upper_limit(events, confidence=0.9, spectrum_cdf=None, n=1000, *, mu_scan_start=1.0, mu_scan_stop=None, mu_scan_step=1.0)

confidence-level optimum-interval upper limit on mu.

Finds mu such that extremeness_of_opt_itv_stat(C_max(events, mu), mu) == confidence -- i.e. where the observed C_max equals :math:\bar C_\mathrm{max}(\text{confidence}, \mu).

The excess extremeness - confidence increases with mu (a larger proposed signal makes the observed emptiness look more anomalous). It is evaluated on a fixed mu grid whose tables are cached on this instance, and the zero crossing is interpolated between the bracketing grid points, so the result is deterministic given the generator seed.

Parameters:

Name Type Description Default
events array_like

Observed event positions (raw; transformed via spectrum_cdf).

 required
confidence float

Confidence level, e.g. 0.9.

 0.9
spectrum_cdf callable

Signal CDF mapping events to cumulants. Identity by default.

 None
n int

Monte-Carlo trials per grid mu. Raise n for a less noisy limit.

 1000
mu_scan_start float

The mu grid. mu_scan_stop defaults to roughly twice the event count (never below mu_scan_start + 10).

 1.0
mu_scan_stop float

The mu grid. mu_scan_stop defaults to roughly twice the event count (never below mu_scan_start + 10).

 1.0
mu_scan_step float

The mu grid. mu_scan_stop defaults to roughly twice the event count (never below mu_scan_start + 10).

 1.0

Returns:

Type Description
float

The upper limit on mu.

Raises:

Type Description
RuntimeError

If no exclusion is bracketed within the scan range; widen mu_scan_start / mu_scan_stop as the message suggests.
Source code in optimum_interval/montecarlo.py 
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
def upper_limit(
    self,
    events: np.ndarray,
    confidence: float = 0.9,
    spectrum_cdf: SpectrumCdf | None = None,
    n: int = 1000,
    *,
    mu_scan_start: float = 1.0,
    mu_scan_stop: float | None = None,
    mu_scan_step: float = 1.0,
) -> float:
    """``confidence``-level optimum-interval upper limit on ``mu``.

    Finds ``mu`` such that ``extremeness_of_opt_itv_stat(C_max(events, mu),
    mu) == confidence`` -- i.e. where the observed C_max equals
    :math:`\\bar C_\\mathrm{max}(\\text{confidence}, \\mu)`.

    The excess ``extremeness - confidence`` increases with ``mu`` (a larger
    proposed signal makes the observed emptiness look more anomalous).  It is
    evaluated on a fixed ``mu`` grid whose tables are cached on this instance,
    and the zero crossing is interpolated between the bracketing grid points,
    so the result is deterministic given the generator seed.

    Parameters
    ----------
    events : array_like
        Observed event positions (raw; transformed via ``spectrum_cdf``).
    confidence : float, optional
        Confidence level, e.g. ``0.9``.
    spectrum_cdf : callable, optional
        Signal CDF mapping ``events`` to cumulants.  Identity by default.
    n : int, optional
        Monte-Carlo trials per grid ``mu``.  Raise ``n`` for a less noisy limit.
    mu_scan_start, mu_scan_stop, mu_scan_step : float, optional
        The ``mu`` grid.  ``mu_scan_stop`` defaults to roughly twice the event
        count (never below ``mu_scan_start + 10``).

    Returns
    -------
    float
        The upper limit on ``mu``.

    Raises
    ------
    RuntimeError
        If no exclusion is bracketed within the scan range; widen
        ``mu_scan_start`` / ``mu_scan_stop`` as the message suggests.
    """
    events = np.asarray(events, dtype=float)

    def excess(mu: float) -> float:
        self.generate(mu, n)  # cached per mu -> deterministic given the seed
        stat = self.optimum_interval_statistic(events, mu, spectrum_cdf)
        return self.extremeness_of_opt_itv_stat(stat, mu) - confidence

    def interpolate(mu_lo, g_lo, mu_hi, g_hi):
        if g_hi == g_lo:
            return float(mu_hi)
        return float(mu_lo + (0.0 - g_lo) * (mu_hi - mu_lo) / (g_hi - g_lo))

    if mu_scan_stop is None:
        mu_scan_stop = max(2.0 * events.size + 5.0, mu_scan_start + 10.0)

    prev_mu, prev_g = None, None
    for mu in np.arange(mu_scan_start, mu_scan_stop + mu_scan_step, mu_scan_step):
        mu = float(mu)
        g = excess(mu)
        if g >= 0.0:
            if prev_mu is None:
                # Limit lies below the first probed mu; step down to bracket.
                lo, g_lo = mu, g
                for _ in range(20):
                    lo *= 0.5
                    g_lo = excess(lo)
                    if lo < 1e-3 or g_lo < 0.0:
                        break
                if lo < 1e-3 or g_lo >= 0.0:
                    raise RuntimeError(
                        f"upper limit is below mu_scan_start={mu_scan_start}; "
                        f"lower it."
                    )
                prev_mu, prev_g = lo, g_lo
            log.info("bracketed limit in [%s, %s]", prev_mu, mu)
            return interpolate(prev_mu, prev_g, mu, g)
        prev_mu, prev_g = mu, g

    raise RuntimeError(
        f"no {confidence:.0%} exclusion found up to mu={mu_scan_stop}; "
        f"increase mu_scan_stop."
    )

save(path)

Pickle the tables to path (single dict, one payload).

Source code in optimum_interval/montecarlo.py 
248
249
250
251
252
253
254
255
256
def save(self, path: str | Path) -> None:
    """Pickle the tables to ``path`` (single dict, one payload)."""
    payload = {
        "itv_sizes": self.itv_sizes,
        "opt_itvs": self.opt_itvs,
        "n_trials": self.n_trials,
    }
    with open(path, "wb") as f:
        pickle.dump(payload, f)

load(path, rng=None) classmethod

Load tables previously written by :meth:save.

Source code in optimum_interval/montecarlo.py 
258
259
260
261
262
263
264
265
266
267
268
269
@classmethod
def load(
    cls, path: str | Path, rng: np.random.Generator | None = None
) -> OptimumIntervalTable:
    """Load tables previously written by :meth:`save`."""
    table = cls(rng=rng)
    with open(path, "rb") as f:
        payload = pickle.load(f)
    table.itv_sizes = payload["itv_sizes"]
    table.opt_itvs = payload["opt_itvs"]
    table.n_trials = payload["n_trials"]
    return table

c0(x, mu)

Probability the maximum gap is smaller than x (Yellin Eq. 2).

Parameters:

Name Type Description Default
x float or array_like

Maximum-gap size in expected events (0 <= x <= mu).

 required
mu float

Total expected number of events over the analysis range.

 required

Returns:

Type Description
float or ndarray

:math:C_0(x, \mu) \in [0, 1], monotonically increasing in x.
Source code in optimum_interval/analytic.py 
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def c0(x, mu):
    r"""Probability the maximum gap is smaller than ``x`` (Yellin Eq. 2).

    Parameters
    ----------
    x : float or array_like
        Maximum-gap size in **expected events** (``0 <= x <= mu``).
    mu : float
        Total expected number of events over the analysis range.

    Returns
    -------
    float or numpy.ndarray
        :math:`C_0(x, \mu) \in [0, 1]`, monotonically increasing in ``x``.
    """
    result = _c0_vec(x, mu)
    return float(result) if np.ndim(x) == 0 and np.ndim(mu) == 0 else result

max_gap_upper_limit(max_gap_fraction, confidence=0.9)

Maximum-gap (:math:C_0) upper limit on mu for one experiment.

Solve :math:C_0(\mu\,f, \mu) = C for mu, where f is the observed maximum-gap size as a fraction of the range (its expected-event size is mu * f). C_0 increases with mu here, so the root is unique.

Parameters:

Name Type Description Default
max_gap_fraction float

Largest gap between adjacent events (including the range endpoints) in cumulant space, in [0, 1].

 required
confidence float

Confidence level, e.g. 0.9.

 0.9
Source code in optimum_interval/analytic.py 
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
def max_gap_upper_limit(max_gap_fraction: float, confidence: float = 0.9) -> float:
    r"""Maximum-gap (:math:`C_0`) upper limit on ``mu`` for one experiment.

    Solve :math:`C_0(\mu\,f, \mu) = C` for ``mu``, where ``f`` is the observed
    maximum-gap size as a *fraction* of the range (its expected-event size is
    ``mu * f``).  ``C_0`` increases with ``mu`` here, so the root is unique.

    Parameters
    ----------
    max_gap_fraction : float
        Largest gap between adjacent events (including the range endpoints) in
        cumulant space, in ``[0, 1]``.
    confidence : float, optional
        Confidence level, e.g. ``0.9``.
    """
    f = float(max_gap_fraction)
    if not 0.0 < f <= 1.0:
        raise ValueError(f"max_gap_fraction must be in (0, 1]; got {f}")

    def g(mu: float) -> float:
        return _c0_scalar(mu * f, mu) - confidence

    hi = max(1.0, 1.0 / f)
    while g(hi) < 0.0 and hi < 1e6:
        hi *= 2.0
    return float(brentq(g, 1e-9, hi))

poisson_upper_limit(n_observed, confidence=0.9)

Classical Poisson upper limit on the mean given n_observed events.

The standard total-count limit: mu_up is the mean for which a fraction confidence of experiments would see more than n_observed events, i.e. :math:\sum_{k=0}^{n} e^{-\mu_\text{up}}\mu_\text{up}^k/k! = 1-C. This is the inverse incomplete Gamma function gammaincinv(n+1, confidence). (Used only for the method-comparison figures; e.g. n=0 gives 2.3026.)

Source code in optimum_interval/analytic.py 
152
153
154
155
156
157
158
159
160
161
def poisson_upper_limit(n_observed: int, confidence: float = 0.9) -> float:
    r"""Classical Poisson upper limit on the mean given ``n_observed`` events.

    The standard total-count limit: ``mu_up`` is the mean for which a fraction
    ``confidence`` of experiments would see *more* than ``n_observed`` events,
    i.e. :math:`\sum_{k=0}^{n} e^{-\mu_\text{up}}\mu_\text{up}^k/k! = 1-C`.  This
    is the inverse incomplete Gamma function ``gammaincinv(n+1, confidence)``.
    (Used only for the method-comparison figures; e.g. ``n=0`` gives ``2.3026``.)
    """
    return float(gammaincinv(n_observed + 1, confidence))

x0(confidence, mu)

Invert :func:c0: the gap size at which c0(x, mu) == confidence.

This is Yellin's :math:x_0(C, \mu). For a max-gap 90% upper limit one raises mu until the observed gap equals x0(0.9, mu).

Parameters:

Name Type Description Default
confidence float

Target confidence level, e.g. 0.9.

 required
mu float

Total expected number of events.

 required

Returns:

Type Description
float

The gap size x (in expected events) solving c0(x, mu) = confidence.

Raises:

Type Description
ValueError

If confidence exceeds the largest attainable C_0, which is C_0(mu, mu) = 1 - e^{-mu}. For a 90% level this happens when mu < 2.3026 -- exactly the regime where no 90% CL exists (Appendix B).
Source code in optimum_interval/analytic.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def x0(confidence: float, mu: float) -> float:
    r"""Invert :func:`c0`: the gap size at which ``c0(x, mu) == confidence``.

    This is Yellin's :math:`x_0(C, \mu)`.  For a max-gap 90% upper limit one
    raises ``mu`` until the observed gap equals ``x0(0.9, mu)``.

    Parameters
    ----------
    confidence : float
        Target confidence level, e.g. ``0.9``.
    mu : float
        Total expected number of events.

    Returns
    -------
    float
        The gap size ``x`` (in expected events) solving ``c0(x, mu) = confidence``.

    Raises
    ------
    ValueError
        If ``confidence`` exceeds the largest attainable ``C_0``, which is
        ``C_0(mu, mu) = 1 - e^{-mu}``.  For a 90% level this happens when
        ``mu < 2.3026`` -- exactly the regime where no 90% CL exists (Appendix B).
    """
    c_max = _c0_scalar(mu, mu)  # = 1 - e^{-mu}, the largest attainable C_0
    if confidence >= c_max:
        raise ValueError(
            f"C_0 cannot reach {confidence} for mu={mu}: its maximum is "
            f"{c_max:.4f} at x=mu. (A 90% max-gap limit requires mu > 2.3026.)"
        )
    # c0 is 0 at x->0 and c_max at x=mu, so the root is bracketed by (tiny, mu).
    return float(brentq(lambda x: _c0_scalar(x, mu) - confidence, 1e-12, mu))

cumulant_points(events, spectrum_cdf=None)

Map observed events to cumulant space and append the range endpoints.

This is the transform every real measurement must undergo before calling :func:k_largest_intervals: apply the (normalized) spectrum CDF so the signal becomes uniform on [0, 1], then add the experimental-range boundaries 0 and 1 (which act as interval delimiters exactly like events).

The Monte-Carlo path builds equivalent point lists directly (interior uniforms plus 0 and 1), so the two paths are consistent.

Parameters:

Name Type Description Default
events array_like

Raw event positions (e.g. recoil energies) inside the analysis range.

 required
spectrum_cdf callable

Monotonic CDF mapping events onto [0, 1]. Defaults to the identity (use this when events are already cumulants).

 None

Returns:

Type Description
ndarray

Sorted cumulants with 0 prepended and 1 appended.

Raises:

Type Description
ValueError

If the resulting cumulants fall outside [0, 1] or are not non-decreasing -- i.e. spectrum_cdf is not a CDF normalized so the analysis range maps onto [0, 1].
Source code in optimum_interval/intervals.py 
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
def cumulant_points(
    events: np.ndarray,
    spectrum_cdf: SpectrumCdf | None = None,
) -> np.ndarray:
    """Map observed events to cumulant space and append the range endpoints.

    This is the transform every *real* measurement must undergo before calling
    :func:`k_largest_intervals`: apply the (normalized) spectrum CDF so the
    signal becomes uniform on ``[0, 1]``, then add the experimental-range
    boundaries 0 and 1 (which act as interval delimiters exactly like events).

    The Monte-Carlo path builds equivalent point lists directly (interior
    uniforms plus 0 and 1), so the two paths are consistent.

    Parameters
    ----------
    events : array_like
        Raw event positions (e.g. recoil energies) inside the analysis range.
    spectrum_cdf : callable, optional
        Monotonic CDF mapping ``events`` onto ``[0, 1]``.  Defaults to the
        identity (use this when ``events`` are already cumulants).

    Returns
    -------
    numpy.ndarray
        Sorted cumulants with 0 prepended and 1 appended.

    Raises
    ------
    ValueError
        If the resulting cumulants fall outside ``[0, 1]`` or are not
        non-decreasing -- i.e. ``spectrum_cdf`` is not a CDF normalized so the
        analysis range maps onto ``[0, 1]``.
    """
    e = np.sort(np.asarray(events, dtype=float))
    if spectrum_cdf is not None:
        e = np.asarray(spectrum_cdf(e), dtype=float)
    if e.size:
        tol = 1e-9
        if e[0] < -tol or e[-1] > 1.0 + tol:
            raise ValueError(
                f"cumulants must lie in [0, 1]; got [{e.min():.4g}, {e.max():.4g}]. "
                "Is spectrum_cdf normalized so the analysis range maps onto [0, 1]?"
            )
        if np.any(np.diff(e) < -tol):
            raise ValueError("cumulants must be non-decreasing (spectrum_cdf monotonic).")
        e = np.clip(e, 0.0, 1.0)  # tidy tiny float excursions before adding 0/1
    return np.concatenate(([0.0], e, [1.0]))

k_largest_intervals(points, spectrum_cdf=None)

Size of the largest interval containing exactly k events, for each k.

An "interval" is a stretch of the parameter (energy) axis; its size is the expected fraction of signal events it contains. The k-largest interval of a run is the largest interval that happens to contain exactly k observed events -- i.e. an unusually empty region. Intervals are delimited by observed events (and, for a real measurement, by the range endpoints, which the caller must include; see :func:cumulant_points).

Parameters:

Name Type Description Default
points array_like

1-D event positions. If spectrum_cdf is None the points are assumed to already be in cumulant space (uniform [0, 1]). Boundary points (0 and 1) are not added here -- the caller decides, so that Monte-Carlo trials and real data are treated identically and boundaries are never double-counted.

 required
spectrum_cdf callable

Monotonic CDF mapping raw positions to cumulants. Applied after sorting. Defaults to the identity.

 None

Returns:

Type Description
dict[int, float]

Mapping k -> size where size is the largest interval containing exactly k events. k runs from 0 (the maximum gap) up to len(points) - 2.
Notes

For a sorted point list p the largest interval containing k events is max_i (p[i + k + 1] - p[i]) -- the widest span between two points that are k + 1 apart in the sorted order (hence k points strictly between them). k = 0 reduces to the maximum gap of the max-gap method.

Examples:

>>> k_largest_intervals(np.array([0.0, 0.1, 0.2, 0.84, 0.85]))[0]
0.64
Source code in optimum_interval/intervals.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
def k_largest_intervals(
    points: np.ndarray,
    spectrum_cdf: SpectrumCdf | None = None,
) -> dict[int, float]:
    """Size of the largest interval containing exactly ``k`` events, for each ``k``.

    An "interval" is a stretch of the parameter (energy) axis; its **size** is
    the expected fraction of signal events it contains.  The *k-largest*
    interval of a run is the largest interval that happens to contain exactly
    ``k`` observed events -- i.e. an unusually empty region.  Intervals are
    delimited by observed events (and, for a real measurement, by the range
    endpoints, which the caller must include; see :func:`cumulant_points`).

    Parameters
    ----------
    points : array_like
        1-D event positions.  If ``spectrum_cdf`` is ``None`` the points are
        assumed to already be in cumulant space (uniform ``[0, 1]``).  Boundary
        points (0 and 1) are **not** added here -- the caller decides, so that
        Monte-Carlo trials and real data are treated identically and boundaries
        are never double-counted.
    spectrum_cdf : callable, optional
        Monotonic CDF mapping raw positions to cumulants.  Applied *after*
        sorting.  Defaults to the identity.

    Returns
    -------
    dict[int, float]
        Mapping ``k -> size`` where ``size`` is the largest interval containing
        exactly ``k`` events.  ``k`` runs from 0 (the maximum gap) up to
        ``len(points) - 2``.

    Notes
    -----
    For a sorted point list ``p`` the largest interval containing ``k`` events
    is ``max_i (p[i + k + 1] - p[i])`` -- the widest span between two points
    that are ``k + 1`` apart in the sorted order (hence ``k`` points strictly
    between them).  ``k = 0`` reduces to the maximum gap of the max-gap method.

    Examples
    --------
    >>> k_largest_intervals(np.array([0.0, 0.1, 0.2, 0.84, 0.85]))[0]
    0.64
    """
    p = np.sort(np.asarray(points, dtype=float))  # copy => no in-place side effect
    if spectrum_cdf is not None:
        p = spectrum_cdf(p)

    n_points = p.size
    sizes: dict[int, float] = {}
    for k in range(n_points - 1):
        gap = k + 1
        sizes[k] = float(np.max(p[gap:] - p[:-gap]))
    return sizes

excluded_interval(params, extremeness, level=0.95)

(low, high) edges of the excluded interval along a 1-D parameter scan.

extremeness is the :func:scan_extremeness value at each of the (positive, increasing) params; edges are log-interpolated crossings of level. Returns (nan, nan) when nothing is excluded; an edge saturates at the end of the scan if the crossing lies beyond it.

Source code in optimum_interval/scanning.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
def excluded_interval(params, extremeness, level=0.95):
    """(low, high) edges of the excluded interval along a 1-D parameter scan.

    ``extremeness`` is the :func:`scan_extremeness` value at each of the
    (positive, increasing) ``params``; edges are log-interpolated crossings
    of ``level``. Returns ``(nan, nan)`` when nothing is excluded; an edge
    saturates at the end of the scan if the crossing lies beyond it.
    """
    params = np.asarray(params, dtype=float)
    ps = np.asarray(extremeness, dtype=float)
    above = ps >= level
    if not above.any():
        return np.nan, np.nan
    idx = np.where(above)[0]
    if idx[0] > 0:
        lo = np.interp(level, ps[idx[0] - 1: idx[0] + 1],
                       np.log10(params[idx[0] - 1: idx[0] + 1]))
    else:
        lo = np.log10(params[0])
    if idx[-1] < len(params) - 1:
        hi = np.interp(-level, -ps[idx[-1]: idx[-1] + 2],
                       np.log10(params[idx[-1]: idx[-1] + 2]))
    else:
        hi = np.log10(params[-1])
    return 10 ** lo, 10 ** hi

new_table(seed=0)

A seeded Monte-Carlo calibration table (reusable across a whole scan).

Source code in optimum_interval/scanning.py 
35
36
37
def new_table(seed=0):
    """A seeded Monte-Carlo calibration table (reusable across a whole scan)."""
    return OptimumIntervalTable(rng=np.random.default_rng(seed))

round_log(x, dex=0.02)

Round onto a dex-spaced log grid so tables are shared between points.

Source code in optimum_interval/scanning.py 
40
41
42
def round_log(x, dex=0.02):
    """Round onto a ``dex``-spaced log grid so tables are shared between points."""
    return 10 ** (np.round(np.log10(x) / dex) * dex)

scan_extremeness(table, events, x, rate, exposure, n=2500, mu_floor=0.2, mu_cap=40.0, mu_dex=0.02)

Optimum-interval extremeness of events for one scan point's spectrum.

Returns (p, mu) where p is the probability that a background-free pseudo-experiment under this hypothesis looks less extreme than the data (so the limit is where p crosses the confidence level), and mu the expected counts. Events outside the spectrum's support window are dropped automatically: they cannot be signal at this scan point.

Shortcuts: mu < mu_floor returns p = 0 (nothing expected); mu > mu_cap returns p = 1 (with few observed events the exclusion is overwhelming; no Monte Carlo needed — raise the cap if you scan with many events). mu is rounded onto a mu_dex log grid so calibration tables are reused across the scan.

Source code in optimum_interval/scanning.py 
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
def scan_extremeness(table, events, x, rate, exposure,
                     n=2500, mu_floor=0.2, mu_cap=40.0, mu_dex=0.02):
    """Optimum-interval extremeness of ``events`` for one scan point's spectrum.

    Returns ``(p, mu)`` where ``p`` is the probability that a background-free
    pseudo-experiment under this hypothesis looks *less* extreme than the
    data (so the limit is where ``p`` crosses the confidence level), and
    ``mu`` the expected counts. Events outside the spectrum's support window
    are dropped automatically: they cannot be signal at this scan point.

    Shortcuts: ``mu < mu_floor`` returns ``p = 0`` (nothing expected);
    ``mu > mu_cap`` returns ``p = 1`` (with few observed events the exclusion
    is overwhelming; no Monte Carlo needed — raise the cap if you scan with
    many events). ``mu`` is rounded onto a ``mu_dex`` log grid so calibration
    tables are reused across the scan.
    """
    spec = spectrum_from_rate(x, rate, exposure)
    if spec is None:
        return 0.0, 0.0
    mu, cdf, x_lo, x_hi = spec
    if mu < mu_floor:
        return 0.0, mu
    if mu > mu_cap:
        return 1.0, mu
    inside = np.asarray(events, dtype=float)
    inside = inside[(inside > x_lo) & (inside < x_hi)]
    mu_r = round_log(mu, mu_dex)
    table.generate(mu_r, n)
    stat = table.optimum_interval_statistic(inside, mu_r, spectrum_cdf=cdf)
    return table.extremeness_of_opt_itv_stat(stat, mu_r), mu

spectrum_from_rate(x, rate, exposure)

Expected counts and normalized spectrum CDF from a differential rate.

Parameters:

Name Type Description Default
x arrays

Differential rate dR/dx (per unit time per unit x) sampled on the grid x. Trailing zero-rate points are trimmed, so the support window ends at the last point with positive rate.

 required
rate arrays

Differential rate dR/dx (per unit time per unit x) sampled on the grid x. Trailing zero-rate points are trimmed, so the support window ends at the last point with positive rate.

 required
exposure float

Livetime in the inverse time unit of rate.

 required

Returns:

Type Description
(mu, cdf, x_lo, x_hi), or ``None`` if the rate has no support.
Source code in optimum_interval/scanning.py 
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
def spectrum_from_rate(x, rate, exposure):
    """Expected counts and normalized spectrum CDF from a differential rate.

    Parameters
    ----------
    x, rate : arrays
        Differential rate ``dR/dx`` (per unit time per unit ``x``) sampled on
        the grid ``x``. Trailing zero-rate points are trimmed, so the support
        window ends at the last point with positive rate.
    exposure : float
        Livetime in the inverse time unit of ``rate``.

    Returns
    -------
    (mu, cdf, x_lo, x_hi), or ``None`` if the rate has no support.
    """
    x = np.asarray(x, dtype=float)
    rate = np.maximum(np.asarray(rate, dtype=float), 0.0)
    if not np.any(rate > 0):
        return None
    hi = np.max(np.where(rate > 0)[0])
    x, rate = x[: hi + 1], rate[: hi + 1]
    mu = float(np.trapezoid(rate, x)) * exposure
    cum = np.concatenate(
        [[0.0], np.cumsum(0.5 * (rate[1:] + rate[:-1]) * np.diff(x))])
    if cum[-1] <= 0:
        return None
    cdf = spectrum_cdf_from_samples(x, cum / cum[-1])
    return mu, cdf, x[0], x[-1]

spectrum_cdf_from_pdf(pdf, e_min, e_max, n_grid=2048)

Normalized CDF on [e_min, e_max] from a signal density pdf.

pdf need not be normalized -- only its shape matters (the overall rate, i.e. mu, is what the method limits and must not be folded in). The density is cumulatively integrated (trapezoidal) on a fine grid and the result is normalized so cdf(e_min) == 0 and cdf(e_max) == 1.

Parameters:

Name Type Description Default
pdf callable

Non-negative density dN/dE evaluated on an array of energies.

 required
e_min float

Analysis-window bounds.

 required
e_max float

Analysis-window bounds.

 required
n_grid int

Integration/interpolation grid size.

 2048
Source code in optimum_interval/spectra.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
def spectrum_cdf_from_pdf(
    pdf: Callable[[np.ndarray], np.ndarray],
    e_min: float,
    e_max: float,
    n_grid: int = 2048,
) -> SpectrumCdf:
    """Normalized CDF on ``[e_min, e_max]`` from a signal density ``pdf``.

    ``pdf`` need not be normalized -- only its *shape* matters (the overall rate,
    i.e. ``mu``, is what the method limits and must not be folded in).  The
    density is cumulatively integrated (trapezoidal) on a fine grid and the
    result is normalized so ``cdf(e_min) == 0`` and ``cdf(e_max) == 1``.

    Parameters
    ----------
    pdf : callable
        Non-negative density ``dN/dE`` evaluated on an array of energies.
    e_min, e_max : float
        Analysis-window bounds.
    n_grid : int, optional
        Integration/interpolation grid size.
    """
    grid = np.linspace(e_min, e_max, n_grid)
    density = np.asarray(pdf(grid), dtype=float)
    if np.any(density < 0.0):
        raise ValueError("pdf must be non-negative on [e_min, e_max]")
    cumulative = np.concatenate(
        [[0.0], np.cumsum(0.5 * (density[1:] + density[:-1]) * np.diff(grid))]
    )
    total = cumulative[-1]
    if total <= 0.0:
        raise ValueError("pdf integrates to zero over [e_min, e_max]")
    cumulative /= total

    def cdf(energies):
        e = np.clip(np.asarray(energies, dtype=float), e_min, e_max)
        return np.interp(e, grid, cumulative)

    return cdf

spectrum_cdf_from_samples(energies, cdf_values)

Normalized spectrum_cdf interpolating tabulated (energy, CDF) points.

Useful when the CDF is only known at sampled energies (e.g. from a simulation). cdf_values must be non-decreasing; they are rescaled so the endpoints map to exactly [0, 1].

Source code in optimum_interval/spectra.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
def spectrum_cdf_from_samples(energies, cdf_values) -> SpectrumCdf:
    """Normalized ``spectrum_cdf`` interpolating tabulated ``(energy, CDF)`` points.

    Useful when the CDF is only known at sampled energies (e.g. from a
    simulation).  ``cdf_values`` must be non-decreasing; they are rescaled so the
    endpoints map to exactly ``[0, 1]``.
    """
    e = np.asarray(energies, dtype=float)
    c = np.asarray(cdf_values, dtype=float)
    if e.ndim != 1 or e.shape != c.shape:
        raise ValueError("energies and cdf_values must be 1-D and the same length")
    if np.any(np.diff(e) <= 0.0):
        raise ValueError("energies must be strictly increasing")
    if np.any(np.diff(c) < 0.0):
        raise ValueError("cdf_values must be non-decreasing")
    c = (c - c[0]) / (c[-1] - c[0])

    def cdf(x):
        x = np.clip(np.asarray(x, dtype=float), e[0], e[-1])
        return np.interp(x, e, c)

    return cdf