Skip to content

Conformal Test Martingales

Martingales

online_cp.martingale.PluginMartingale

Bases: ConformalTestMartingale

Plugin martingale using a betting strategy for density estimation.

The martingale wraps a BettingStrategy whose bet(p) method provides the predictive density (betting function) at each step. The strategy is updated after the bet — predict-then-learn order — preserving the martingale property.

Protocol per step: 1. Predict: evaluate strategy.bet(p) (uses past data only) 2. Accumulate: logM += log(bet) 3. Learn: call strategy.update(p) 4. Expose: set b_n and B_n for the next step

For cautious behaviour during early steps (small sample), use :class:ExpertAggregationStrategy to mix the primary strategy with a uniform baseline, or wrap in a :class:SleeperStayer.

Parameters:

Name Type Description Default
betting_strategy BettingStrategy or type

An instantiated strategy, or a class to be instantiated with kwargs.

 GaussianKDE
**kwargs Any

Passed to the strategy constructor if a class is given.

 {}

Examples:

>>> from online_cp.betting import FixedStrategy
>>> strat = FixedStrategy(pdf=lambda x: 1.5 if x < 0.5 else 0.5, check_integration=False)
>>> m = PluginMartingale(betting_strategy=strat)
>>> m.update(0.1)
>>> bool(np.isclose(m.M, 1.5))
True
>>> m.update(0.9)
>>> bool(np.isclose(m.M, 0.75))
True
Source code in src/online_cp/martingale/jumpers.py 
 19
 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
 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
class PluginMartingale(ConformalTestMartingale):
    """Plugin martingale using a betting strategy for density estimation.

    The martingale wraps a ``BettingStrategy`` whose ``bet(p)`` method provides
    the predictive density (betting function) at each step. The strategy is
    updated *after* the bet — predict-then-learn order — preserving the
    martingale property.

    Protocol per step:
    1. **Predict**: evaluate ``strategy.bet(p)`` (uses past data only)
    2. **Accumulate**: ``logM += log(bet)``
    3. **Learn**: call ``strategy.update(p)``
    4. **Expose**: set ``b_n`` and ``B_n`` for the *next* step

    For cautious behaviour during early steps (small sample), use
    :class:`ExpertAggregationStrategy` to mix the primary strategy with a
    uniform baseline, or wrap in a :class:`SleeperStayer`.

    Parameters
    ----------
    betting_strategy : BettingStrategy or type
        An instantiated strategy, or a class to be instantiated with kwargs.
    **kwargs
        Passed to the strategy constructor if a class is given.

    Examples
    --------
    >>> from online_cp.betting import FixedStrategy
    >>> strat = FixedStrategy(pdf=lambda x: 1.5 if x < 0.5 else 0.5, check_integration=False)
    >>> m = PluginMartingale(betting_strategy=strat)
    >>> m.update(0.1)
    >>> bool(np.isclose(m.M, 1.5))
    True
    >>> m.update(0.9)
    >>> bool(np.isclose(m.M, 0.75))
    True
    """

    def __init__(
        self,
        betting_strategy: type[BettingStrategy] | BettingStrategy = GaussianKDE,
        store_p_values: bool = True,
        **kwargs: Any,
    ) -> None:
        super().__init__(store_p_values)

        if isinstance(betting_strategy, BettingStrategy):
            self.strategy = betting_strategy
        else:
            betting_kwargs = kwargs if kwargs else {}
            self.strategy = betting_strategy(**betting_kwargs)

        # Expose initial b_n / B_n
        self._update_exposed_functions()

    def _update_exposed_functions(self):
        """Set b_n and B_n for the *next* step (using current strategy state)."""
        strategy_snapshot = deepcopy(self.strategy)

        def b_n(x, _s=strategy_snapshot):
            return _s.bet(x)

        def B_n(x, _s=strategy_snapshot):
            return _s.integrate(x)

        self.b_n = b_n
        self.B_n = B_n

    def update(self, p: float) -> None:
        """Bet on one p-value with the current strategy, then learn from it.

        Evaluates the betting density at ``p``, multiplies it into the
        martingale (guarding against invalid densities), records the value, and
        updates the underlying betting strategy.

        Parameters
        ----------
        p : float
            New p-value in $[0, 1]$.
        """
        # 1. Predict: evaluate current betting function
        b = self.strategy.bet(p)

        # 2. Safeguard: ensure b is a valid positive density
        # Betting functions should be > 0 (probability densities/likelihoods)
        if not np.isfinite(b) or b <= 0:
            warnings.warn(
                f"Betting function returned invalid value b={b} for p={p}. "
                f"Using fallback b=1.0 (uniform betting)",
                RuntimeWarning,
                stacklevel=2
            )
            b = 1.0

        # 3. Accumulate wealth
        self.logM += np.log(b)
        self.log_martingale_values.append(self.logM)

        if self.store_p_values:
            self.p_values.append(p)

        # 4. Learn
        self.strategy.update(p)

        # 5. Mark b_n/B_n stale (lazy recomputation on access)
        self._mark_stale()

update(p: float) -> None

Bet on one p-value with the current strategy, then learn from it.

Evaluates the betting density at p, multiplies it into the martingale (guarding against invalid densities), records the value, and updates the underlying betting strategy.

Parameters:

Name Type Description Default
p float

New p-value in \([0, 1]\).

 required
Source code in src/online_cp/martingale/jumpers.py 
 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
def update(self, p: float) -> None:
    """Bet on one p-value with the current strategy, then learn from it.

    Evaluates the betting density at ``p``, multiplies it into the
    martingale (guarding against invalid densities), records the value, and
    updates the underlying betting strategy.

    Parameters
    ----------
    p : float
        New p-value in $[0, 1]$.
    """
    # 1. Predict: evaluate current betting function
    b = self.strategy.bet(p)

    # 2. Safeguard: ensure b is a valid positive density
    # Betting functions should be > 0 (probability densities/likelihoods)
    if not np.isfinite(b) or b <= 0:
        warnings.warn(
            f"Betting function returned invalid value b={b} for p={p}. "
            f"Using fallback b=1.0 (uniform betting)",
            RuntimeWarning,
            stacklevel=2
        )
        b = 1.0

    # 3. Accumulate wealth
    self.logM += np.log(b)
    self.log_martingale_values.append(self.logM)

    if self.store_p_values:
        self.p_values.append(p)

    # 4. Learn
    self.strategy.update(p)

    # 5. Mark b_n/B_n stale (lazy recomputation on access)
    self._mark_stale()

online_cp.martingale.SimpleMixtureMartingale

Bases: ConformalTestMartingale

Simple Mixture Martingale using the incomplete gamma function.

This is the canonical "parameter-free" test martingale that averages over all power alternatives t^epsilon with epsilon ~ Exp(1). It has a closed-form solution based on the regularized incomplete gamma function.

After each step, b_n and B_n are set analytically: - b_n(p) = M_n(p) / M_{n-1} where M_n(p) is the martingale value if the next observation were p. - B_n(p) = integral of b_n from 0 to p.

Examples:

>>> sm = SimpleMixtureMartingale()
>>> sm.update(0.01)
>>> sm.update(0.01)
>>> bool(sm.M > 1)
True
>>> sm.update(1.0)
>>> bool(sm.M < sm.martingale_values[-2])
True
Source code in src/online_cp/martingale/jumpers.py 
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
class SimpleMixtureMartingale(ConformalTestMartingale):
    """Simple Mixture Martingale using the incomplete gamma function.

    This is the canonical "parameter-free" test martingale that averages over
    all power alternatives t^epsilon with epsilon ~ Exp(1). It has a closed-form
    solution based on the regularized incomplete gamma function.

    After each step, ``b_n`` and ``B_n`` are set analytically:
    - b_n(p) = M_n(p) / M_{n-1} where M_n(p) is the martingale value if the
      next observation were p.
    - B_n(p) = integral of b_n from 0 to p.

    Examples
    --------
    >>> sm = SimpleMixtureMartingale()
    >>> sm.update(0.01)
    >>> sm.update(0.01)
    >>> bool(sm.M > 1)
    True
    >>> sm.update(1.0)
    >>> bool(sm.M < sm.martingale_values[-2])
    True
    """

    def __init__(self, **kwargs):
        super().__init__(**kwargs)
        self.n = 0
        self.sum_log_p = 0.0

    def _compute_logM(self, n, sum_log_p):
        """Compute log-martingale value from sufficient statistics."""
        if n == 0:
            return 0.0
        L = sum_log_p
        if np.isclose(L, 0):
            return -np.log(n + 1)
        arg = -L
        log_gamma_n_plus_1 = gammaln(n + 1)
        val_gammainc = gammainc(n + 1, arg)
        if val_gammainc <= 0:
            log_incomplete_gamma = -700  # underflow floor (~1e-304)
        else:
            log_incomplete_gamma = log_gamma_n_plus_1 + np.log(val_gammainc)
        return -L + log_incomplete_gamma - (n + 1) * np.log(arg)

    def _update_exposed_functions(self):
        """Set analytic b_n and B_n for the next step.

        b_n(p) = M(n+1, sum_log_p + log(p)) / M(n, sum_log_p)
        which simplifies to an analytic ratio of incomplete gamma terms.
        """
        n = self.n
        slp = self.sum_log_p
        current_logM = self.logM

        def b_n(p, _n=n, _slp=slp, _logM=current_logM):
            p_c = np.clip(p, 1e-12, 1.0)
            next_logM = self._compute_logM(_n + 1, _slp + np.log(p_c))
            return np.exp(next_logM - _logM)

        def B_n(p, _n=n, _slp=slp, _logM=current_logM):
            # B_n(p) = integral_0^p b_n(t) dt
            # We compute this via numerical integration (the ratio doesn't simplify to closed-form CDF)
            if p <= 1e-12:
                return 0.0
            if p >= 1.0 - 1e-12:
                return 1.0
            val, _ = quad(b_n, 1e-12, p, limit=50)
            return val

        self.b_n = b_n
        self.B_n = B_n

    def update(self, p: float) -> None:
        """Advance the simple mixture martingale by one p-value.

        Accumulates $\\sum \\log p$ and recomputes the closed-form mixture
        log-martingale from the running count and log-sum.

        Parameters
        ----------
        p : float
            New p-value in $[0, 1]$.
        """
        self.n += 1
        p_clipped = np.clip(p, 1e-12, 1.0)
        self.sum_log_p += np.log(p_clipped)

        self.logM = self._compute_logM(self.n, self.sum_log_p)

        if self.store_p_values:
            self.p_values.append(p)
        self.log_martingale_values.append(self.logM)

        self._mark_stale()

update(p: float) -> None

Advance the simple mixture martingale by one p-value.

Accumulates \(\sum \log p\) and recomputes the closed-form mixture log-martingale from the running count and log-sum.

Parameters:

Name Type Description Default
p float

New p-value in \([0, 1]\).

 required
Source code in src/online_cp/martingale/jumpers.py 
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
def update(self, p: float) -> None:
    """Advance the simple mixture martingale by one p-value.

    Accumulates $\\sum \\log p$ and recomputes the closed-form mixture
    log-martingale from the running count and log-sum.

    Parameters
    ----------
    p : float
        New p-value in $[0, 1]$.
    """
    self.n += 1
    p_clipped = np.clip(p, 1e-12, 1.0)
    self.sum_log_p += np.log(p_clipped)

    self.logM = self._compute_logM(self.n, self.sum_log_p)

    if self.store_p_values:
        self.p_values.append(p)
    self.log_martingale_values.append(self.logM)

    self._mark_stale()

online_cp.martingale.SimpleJumper

Bases: ConformalTestMartingale

Simple Jumper betting martingale (Algorithm 8.1 of ALRW2).

Uses a set of experts indexed by epsilon with betting functions f_epsilon(p) = 1 + epsilon*(p - 0.5). A Markov chain with jump rate J tracks the best expert, enabling adaptation to changing alternatives.

Parameters:

Name Type Description Default
J float

Jump rate (probability of switching expert per step).

 0.01
E list of float or None

Expert grid. Default is [-1, -0.5, 0, 0.5, 1] (Algorithm 8.1).

 None
References

Vovk, Gammerman & Shafer (2022). Algorithmic Learning in a Random World, 2nd edition, Algorithm 8.1. Cambridge University Press.

Examples:

>>> sj = SimpleJumper(J=0.1)
>>> for _ in range(5):
...     sj.update(0.0001)
>>> bool(sj.M > 1.0)
True
Source code in src/online_cp/martingale/jumpers.py 
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
class SimpleJumper(ConformalTestMartingale):
    """Simple Jumper betting martingale (Algorithm 8.1 of ALRW2).

    Uses a set of experts indexed by epsilon with betting functions
    f_epsilon(p) = 1 + epsilon*(p - 0.5). A Markov chain with jump rate J
    tracks the best expert, enabling adaptation to changing alternatives.

    Parameters
    ----------
    J : float
        Jump rate (probability of switching expert per step).
    E : list of float or None
        Expert grid. Default is [-1, -0.5, 0, 0.5, 1] (Algorithm 8.1).

    References
    ----------
    Vovk, Gammerman & Shafer (2022). *Algorithmic Learning in a Random World*,
    2nd edition, Algorithm 8.1. Cambridge University Press.

    Examples
    --------
    >>> sj = SimpleJumper(J=0.1)
    >>> for _ in range(5):
    ...     sj.update(0.0001)
    >>> bool(sj.M > 1.0)
    True
    """

    def __init__(self, J=0.01, E=None, store_p_values=True, **kwargs):
        super().__init__(store_p_values)
        self.J = J
        if E is None:
            self.E = [-1, -0.5, 0, 0.5, 1]
        else:
            self.E = list(E)
        self._n_experts = len(self.E)
        self.log_C_epsilon = {eps: -np.log(self._n_experts) for eps in self.E}
        self.log_C = 0.0

        self.b_epsilon = lambda u, epsilon: 1 + epsilon * (u - 1 / 2)
        self.B_n_inv = lambda x: x

    def update(self, p: float) -> None:
        r"""Advance the Simple Jumper mixture by one p-value.

        Mixes a sleeping family of constant betting functions indexed by
        $\epsilon$ with jump probability ``J``, updating the per-$\epsilon$ and
        pooled log-wealth in log-space.

        Parameters
        ----------
        p : float
            New p-value in $[0, 1]$.
        """
        if self.store_p_values:
            self.p_values.append(p)

        if self.J == 1:
            log_1_minus_J = -np.inf
        else:
            log_1_minus_J = np.log(1 - self.J)

        log_J_div_E = np.log(self.J / self._n_experts)

        new_log_C_epsilon = {}

        for epsilon in self.E:
            term1 = log_1_minus_J + self.log_C_epsilon[epsilon]
            term2 = log_J_div_E + self.log_C
            log_C_mixed = np.logaddexp(term1, term2)
            bet_val = self.b_epsilon(p, epsilon)
            # Safeguard: betting function must be > 0
            if bet_val <= 0 or not np.isfinite(bet_val):
                warnings.warn(
                    f"Betting function returned invalid value for epsilon={epsilon}, "
                    f"p={p}: b={bet_val}. Using fallback b=1.0",
                    RuntimeWarning,
                    stacklevel=2
                )
                bet_val = 1.0
            new_log_C_epsilon[epsilon] = log_C_mixed + np.log(bet_val)

        self.log_C_epsilon = new_log_C_epsilon
        self.log_C = logsumexp(list(self.log_C_epsilon.values()))

        self.logM = self.log_C
        self.log_martingale_values.append(self.logM)

        # Mark b_n/B_n stale (lazy recomputation on access)
        self._mark_stale()

    def _update_exposed_functions(self):
        """Compute wealth-weighted effective betting function for next step."""
        # b_n(u) = sum_eps w_eps * (1 + eps*(u - 0.5))
        #        = 1 + eps_bar * (u - 0.5)
        # where eps_bar = (1-J) * sum_eps eps * exp(log_C_eps - log_C)
        weights = {eps: np.exp(self.log_C_epsilon[eps] - self.log_C) for eps in self.E}
        epsilon_bar = (1 - self.J) * sum(eps * weights[eps] for eps in self.E)

        self.b_n = lambda u: 1 + epsilon_bar * (u - 1 / 2)
        self.B_n = lambda u: (epsilon_bar / 2) * u**2 + (1 - epsilon_bar / 2) * u
        self.B_n_inv = lambda u: (
            (epsilon_bar - 2) / (2 * epsilon_bar)
            + np.sqrt(epsilon_bar * (8 * u + epsilon_bar - 4) + 4) / (2 * epsilon_bar)
            if abs(epsilon_bar) > 1e-9
            else u
        )

update(p: float) -> None

Advance the Simple Jumper mixture by one p-value.

Mixes a sleeping family of constant betting functions indexed by \(\epsilon\) with jump probability J, updating the per-\(\epsilon\) and pooled log-wealth in log-space.

Parameters:

Name Type Description Default
p float

New p-value in \([0, 1]\).

 required
Source code in src/online_cp/martingale/jumpers.py 
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
def update(self, p: float) -> None:
    r"""Advance the Simple Jumper mixture by one p-value.

    Mixes a sleeping family of constant betting functions indexed by
    $\epsilon$ with jump probability ``J``, updating the per-$\epsilon$ and
    pooled log-wealth in log-space.

    Parameters
    ----------
    p : float
        New p-value in $[0, 1]$.
    """
    if self.store_p_values:
        self.p_values.append(p)

    if self.J == 1:
        log_1_minus_J = -np.inf
    else:
        log_1_minus_J = np.log(1 - self.J)

    log_J_div_E = np.log(self.J / self._n_experts)

    new_log_C_epsilon = {}

    for epsilon in self.E:
        term1 = log_1_minus_J + self.log_C_epsilon[epsilon]
        term2 = log_J_div_E + self.log_C
        log_C_mixed = np.logaddexp(term1, term2)
        bet_val = self.b_epsilon(p, epsilon)
        # Safeguard: betting function must be > 0
        if bet_val <= 0 or not np.isfinite(bet_val):
            warnings.warn(
                f"Betting function returned invalid value for epsilon={epsilon}, "
                f"p={p}: b={bet_val}. Using fallback b=1.0",
                RuntimeWarning,
                stacklevel=2
            )
            bet_val = 1.0
        new_log_C_epsilon[epsilon] = log_C_mixed + np.log(bet_val)

    self.log_C_epsilon = new_log_C_epsilon
    self.log_C = logsumexp(list(self.log_C_epsilon.values()))

    self.logM = self.log_C
    self.log_martingale_values.append(self.logM)

    # Mark b_n/B_n stale (lazy recomputation on access)
    self._mark_stale()

online_cp.martingale.CompositeJumper

Bases: ConformalTestMartingale

Composite Jumper that averages over multiple jump rates.

Examples:

>>> cj = CompositeJumper()
>>> for _ in range(5):
...     cj.update(0.001)
>>> bool(cj.M > 1)
True
Source code in src/online_cp/martingale/jumpers.py 
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
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
class CompositeJumper(ConformalTestMartingale):
    """Composite Jumper that averages over multiple jump rates.

    Examples
    --------
    >>> cj = CompositeJumper()
    >>> for _ in range(5):
    ...     cj.update(0.001)
    >>> bool(cj.M > 1)
    True
    """

    def __init__(self, J=None, store_p_values=True):
        super().__init__(store_p_values)
        if J is None:
            self.J = [10 ** (-4), 10 ** (-3), 10 ** (-2), 10 ** (-1), 1]
        else:
            self.J = J

        self.Jumpers = {j: SimpleJumper(J=j, store_p_values=False) for j in self.J}

    def update(self, p: float) -> None:
        """Advance every sub-jumper and recompute the pooled martingale.

        Updates each component jumper on ``p`` and sets the composite
        log-martingale to the (equal-weight) log-mean of the components.

        Parameters
        ----------
        p : float
            New p-value in $[0, 1]$.
        """
        if self.store_p_values:
            self.p_values.append(p)

        for m in self.Jumpers.values():
            m.update(p)

        log_M_values = [m.logM for m in self.Jumpers.values()]
        self.logM = logsumexp(log_M_values) - np.log(len(self.Jumpers))

        self.log_martingale_values.append(self.logM)

        # Mark b_n/B_n stale (lazy recomputation on access)
        self._mark_stale()

    def _update_exposed_functions(self):
        """Compute wealth-weighted betting function from sub-jumpers."""
        log_M_values = [m.logM for m in self.Jumpers.values()]
        log_sum_M = self.logM + np.log(len(self.Jumpers))
        weights = np.exp(np.array(log_M_values) - log_sum_M)

        current_b_ns = [m.b_n for m in self.Jumpers.values()]
        current_B_ns = [m.B_n for m in self.Jumpers.values()]

        self.b_n = lambda u: np.dot(weights, [f(u) for f in current_b_ns])
        self.B_n = lambda u: np.dot(weights, [F(u) for F in current_B_ns])

update(p: float) -> None

Advance every sub-jumper and recompute the pooled martingale.

Updates each component jumper on p and sets the composite log-martingale to the (equal-weight) log-mean of the components.

Parameters:

Name Type Description Default
p float

New p-value in \([0, 1]\).

 required
Source code in src/online_cp/martingale/jumpers.py 
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
def update(self, p: float) -> None:
    """Advance every sub-jumper and recompute the pooled martingale.

    Updates each component jumper on ``p`` and sets the composite
    log-martingale to the (equal-weight) log-mean of the components.

    Parameters
    ----------
    p : float
        New p-value in $[0, 1]$.
    """
    if self.store_p_values:
        self.p_values.append(p)

    for m in self.Jumpers.values():
        m.update(p)

    log_M_values = [m.logM for m in self.Jumpers.values()]
    self.logM = logsumexp(log_M_values) - np.log(len(self.Jumpers))

    self.log_martingale_values.append(self.logM)

    # Mark b_n/B_n stale (lazy recomputation on access)
    self._mark_stale()

online_cp.martingale.SleeperStayer

Bases: ConformalTestMartingale

Sleeper/Stayer conformal test martingale (Algorithm 9.4 of ALRW2).

Maintains a grid of piecewise-constant betting experts indexed by (a, b) together with a sleeping capital account. At each step, a fraction R of the sleeping capital is redistributed equally to all active experts.

Each expert uses the betting function f_{(a,b)}(p) = b/a if p <= a, else (1-b)/(1-a). This targets change-points where the conformal p-values shift from Uniform to having mass b below threshold a.

Parameters:

Name Type Description Default
R float

Wake-up rate: fraction of sleeping capital redistributed per step.

 0.001
G int

Grid resolution. The grid is {1/G, 2/G, ..., (G-1)/G}^2.

 10
References

Vovk, Gammerman & Shafer (2022). Algorithmic Learning in a Random World, 2nd edition, Algorithm 9.4 (Sleeper). Cambridge University Press.

Examples:

>>> ss = SleeperStayer(R=0.01, G=5)
>>> for _ in range(50):
...     ss.update(0.05)
>>> bool(ss.M > 1.0)
True
Source code in src/online_cp/martingale/sleepers.py 
 12
 13
 14
 15
 16
 17
 18
 19
 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
 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
class SleeperStayer(ConformalTestMartingale):
    """Sleeper/Stayer conformal test martingale (Algorithm 9.4 of ALRW2).

    Maintains a grid of piecewise-constant betting experts indexed by (a, b)
    together with a sleeping capital account. At each step, a fraction R of the
    sleeping capital is redistributed equally to all active experts.

    Each expert uses the betting function f_{(a,b)}(p) = b/a if p <= a, else
    (1-b)/(1-a). This targets change-points where the conformal p-values shift
    from Uniform to having mass b below threshold a.

    Parameters
    ----------
    R : float
        Wake-up rate: fraction of sleeping capital redistributed per step.
    G : int
        Grid resolution. The grid is {1/G, 2/G, ..., (G-1)/G}^2.

    References
    ----------
    Vovk, Gammerman & Shafer (2022). *Algorithmic Learning in a Random World*,
    2nd edition, Algorithm 9.4 (Sleeper). Cambridge University Press.

    Examples
    --------
    >>> ss = SleeperStayer(R=0.01, G=5)
    >>> for _ in range(50):
    ...     ss.update(0.05)
    >>> bool(ss.M > 1.0)
    True
    """

    def __init__(self, R=0.001, G=10, store_p_values=True):
        super().__init__(store_p_values)
        self.R = R
        self.G = G

        # Build the grid: (a, b) pairs with a, b in {1/G, ..., (G-1)/G}
        grid_vals = np.arange(1, G) / G
        self._grid = [(a, b) for a in grid_vals for b in grid_vals]
        self._n_experts = len(self._grid)

        # Precompute betting values for each expert in log-space
        self._left = np.array([b / a for a, b in self._grid])
        self._right = np.array([(1 - b) / (1 - a) for a, b in self._grid])
        self._log_left = np.log(self._left)
        self._log_right = np.log(self._right)
        self._thresholds = np.array([a for a, _ in self._grid])

        # Capital in log-space: all starts in sleeping account
        self._log_S_active = np.full(self._n_experts, -np.inf)  # no capital yet
        self._log_S_sleep = 0.0  # log(1) = 0
        self._log_1_minus_R = math.log(1.0 - R)
        self._log_R_div_n = math.log(R / self._n_experts)
        self._n = 0

    def update(self, p: float) -> None:
        """Advance the Sleeper/Stayer martingale by one p-value.

        Each active threshold expert bets on ``p`` (left vs right of its
        threshold) while a sleeping reserve preserves capital; the martingale is
        the total capital across the sleeping and active experts.

        Parameters
        ----------
        p : float
            New p-value in $[0, 1]$.
        """
        if self.store_p_values:
            self.p_values.append(p)

        self._n += 1

        # Step 1: Bet — add log(bet) to each active expert's log-capital
        log_bets = np.where(p <= self._thresholds, self._log_left, self._log_right)
        self._log_S_active += log_bets

        # Step 2: Output — total capital = exp(log_S_sleep) + sum(exp(log_S_active))
        # Compute in log-space via logsumexp
        log_active_sum = logsumexp(self._log_S_active)
        self.logM = np.logaddexp(self._log_S_sleep, log_active_sum)
        self.log_martingale_values.append(self.logM)

        # Step 3: Redistribute — move fraction R of sleeping capital to active experts
        # log(transfer) = log(R / n_experts) + log_S_sleep
        log_transfer = self._log_R_div_n + self._log_S_sleep
        self._log_S_active = np.logaddexp(self._log_S_active, log_transfer)
        self._log_S_sleep += self._log_1_minus_R

        # Mark b_n/B_n stale (lazy recomputation on access)
        self._mark_stale()

    def _update_exposed_functions(self):
        """Set b_n/B_n as wealth-weighted combination of active experts."""
        log_active_sum = logsumexp(self._log_S_active)
        if np.isfinite(log_active_sum):
            log_weights = self._log_S_active - log_active_sum
            weights = np.exp(log_weights)
            left_vals = self._left
            right_vals = self._right
            thresholds = self._thresholds

            def _b_n(u, w=weights, l=left_vals, r=right_vals, t=thresholds):
                bets = np.where(u <= t, l, r)
                return float(np.dot(w, bets))

            def _B_n(u, w=weights, l=left_vals, r=right_vals, t=thresholds,
                     grid=self._grid):
                # CDF of weighted mixture
                cdfs = np.where(
                    u <= t,
                    l * u,
                    np.array([b + r_i * (u - a) for (a, b), r_i in zip(grid, r)])
                )
                return float(np.dot(w, cdfs))

            self.b_n = _b_n
            self.B_n = _B_n

update(p: float) -> None

Advance the Sleeper/Stayer martingale by one p-value.

Each active threshold expert bets on p (left vs right of its threshold) while a sleeping reserve preserves capital; the martingale is the total capital across the sleeping and active experts.

Parameters:

Name Type Description Default
p float

New p-value in \([0, 1]\).

 required
Source code in src/online_cp/martingale/sleepers.py 
 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
def update(self, p: float) -> None:
    """Advance the Sleeper/Stayer martingale by one p-value.

    Each active threshold expert bets on ``p`` (left vs right of its
    threshold) while a sleeping reserve preserves capital; the martingale is
    the total capital across the sleeping and active experts.

    Parameters
    ----------
    p : float
        New p-value in $[0, 1]$.
    """
    if self.store_p_values:
        self.p_values.append(p)

    self._n += 1

    # Step 1: Bet — add log(bet) to each active expert's log-capital
    log_bets = np.where(p <= self._thresholds, self._log_left, self._log_right)
    self._log_S_active += log_bets

    # Step 2: Output — total capital = exp(log_S_sleep) + sum(exp(log_S_active))
    # Compute in log-space via logsumexp
    log_active_sum = logsumexp(self._log_S_active)
    self.logM = np.logaddexp(self._log_S_sleep, log_active_sum)
    self.log_martingale_values.append(self.logM)

    # Step 3: Redistribute — move fraction R of sleeping capital to active experts
    # log(transfer) = log(R / n_experts) + log_S_sleep
    log_transfer = self._log_R_div_n + self._log_S_sleep
    self._log_S_active = np.logaddexp(self._log_S_active, log_transfer)
    self._log_S_sleep += self._log_1_minus_R

    # Mark b_n/B_n stale (lazy recomputation on access)
    self._mark_stale()

online_cp.martingale.SleeperDrifter

Bases: ConformalTestMartingale

Sleeper/Drifter conformal test martingale (Algorithm 9.5 of ALRW2).

Extension of the Sleeper/Stayer that wakes experts in batches every M steps and uses a drifting threshold that interpolates between the initial guess a and the target b over time.

The drifting threshold for expert (i, a, b) at step n is: a' = (iM/n)a + (1 - iM/n)b

This makes the martingale more sensitive to gradual distribution shifts.

Parameters:

Name Type Description Default
R float

Wake-up rate per batch: fraction of sleeping capital allocated when a new batch wakes up.

 0.001
G int

Grid resolution. The grid is {1/G, 2/G, ..., (G-1)/G}^2.

 10
M int

Batch interval: new experts wake up every M steps.

 100
References

Vovk, Gammerman & Shafer (2022). Algorithmic Learning in a Random World, 2nd edition, Algorithm 9.5 (Drifter). Cambridge University Press.

Examples:

>>> sd = SleeperDrifter(R=0.01, G=5, M=10)
>>> for _ in range(50):
...     sd.update(0.05)
>>> bool(sd.M > 1.0)
True
Source code in src/online_cp/martingale/sleepers.py 
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
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
class SleeperDrifter(ConformalTestMartingale):
    """Sleeper/Drifter conformal test martingale (Algorithm 9.5 of ALRW2).

    Extension of the Sleeper/Stayer that wakes experts in batches every M steps
    and uses a drifting threshold that interpolates between the initial guess a
    and the target b over time.

    The drifting threshold for expert (i, a, b) at step n is:
        a' = (i*M/n)*a + (1 - i*M/n)*b

    This makes the martingale more sensitive to gradual distribution shifts.

    Parameters
    ----------
    R : float
        Wake-up rate per batch: fraction of sleeping capital allocated when
        a new batch wakes up.
    G : int
        Grid resolution. The grid is {1/G, 2/G, ..., (G-1)/G}^2.
    M : int
        Batch interval: new experts wake up every M steps.

    References
    ----------
    Vovk, Gammerman & Shafer (2022). *Algorithmic Learning in a Random World*,
    2nd edition, Algorithm 9.5 (Drifter). Cambridge University Press.

    Examples
    --------
    >>> sd = SleeperDrifter(R=0.01, G=5, M=10)
    >>> for _ in range(50):
    ...     sd.update(0.05)
    >>> bool(sd.M > 1.0)
    True
    """

    def __init__(self, R=0.001, G=10, M=100, store_p_values=True):
        super().__init__(store_p_values)
        self.R = R
        self.G = G
        self._batch_interval = M

        # Grid of (a, b) pairs
        grid_vals = np.arange(1, G) / G
        self._grid = [(a, b) for a in grid_vals for b in grid_vals]
        self._n_grid = len(self._grid)

        # Active experts: dict of (batch_index, grid_index) -> log-capital
        self._experts = {}  # (i, j) -> log_capital
        self._log_S_sleep = 0.0  # log(1) = 0
        self._n = 0
        self._log_prune_threshold = math.log(1e-15)

    def _get_drifted_threshold(self, batch_idx, grid_idx):
        """Compute a' = (i*M/n)*a + (1 - i*M/n)*b for expert (i, a, b)."""
        a, b = self._grid[grid_idx]
        ratio = (batch_idx * self._batch_interval) / self._n
        ratio = min(ratio, 1.0)  # Clamp
        return ratio * a + (1 - ratio) * b

    def update(self, p: float) -> None:
        """Advance the Sleeper/Drifter martingale by one p-value.

        Like :class:`SleeperStayer`, but each expert's threshold *drifts* over
        time, so the bet uses the current drifted threshold of every active
        expert before pooling their capital.

        Parameters
        ----------
        p : float
            New p-value in $[0, 1]$.
        """
        if self.store_p_values:
            self.p_values.append(p)

        self._n += 1

        # Step 1: Bet — update all active experts in log-space
        keys_to_remove = []
        for key, log_capital in self._experts.items():
            batch_idx, grid_idx = key
            a_prime = self._get_drifted_threshold(batch_idx, grid_idx)
            _, b = self._grid[grid_idx]

            # Bet using f_{(a', b)}(p)
            if a_prime <= 0 or a_prime >= 1:
                log_bet_val = 0.0  # log(1) = no bet
            else:
                bet_val = b / a_prime if p <= a_prime else (1 - b) / (1 - a_prime)
                log_bet_val = math.log(bet_val)

            new_log_capital = log_capital + log_bet_val
            if new_log_capital < self._log_prune_threshold:
                keys_to_remove.append(key)
            else:
                self._experts[key] = new_log_capital

        for key in keys_to_remove:
            del self._experts[key]

        # Step 2: Output — total capital in log-space
        if self._experts:
            log_active_sum = logsumexp(list(self._experts.values()))
            self.logM = np.logaddexp(self._log_S_sleep, log_active_sum)
        else:
            self.logM = self._log_S_sleep
        self.log_martingale_values.append(self.logM)

        # Step 3: Wake new batch (if n is divisible by M) — prepares for next step
        if self._n % self._batch_interval == 0:
            batch_idx = self._n // self._batch_interval
            # log(transfer_per_expert) = log(R * M) + log_S_sleep - log(n_grid)
            log_transfer = (math.log(self.R * self._batch_interval)
                           + self._log_S_sleep
                           - math.log(self._n_grid))
            # S_sleep *= (1 - R*M), clamped to avoid negative
            rm = self.R * self._batch_interval
            if rm >= 1.0:
                self._log_S_sleep = -np.inf
            else:
                self._log_S_sleep += math.log(1.0 - rm)
            for j in range(self._n_grid):
                key = (batch_idx, j)
                if key in self._experts:
                    self._experts[key] = np.logaddexp(self._experts[key], log_transfer)
                else:
                    self._experts[key] = log_transfer

        # Mark b_n/B_n stale (lazy recomputation on access)
        self._mark_stale()

    def _update_exposed_functions(self):
        """Set b_n/B_n as wealth-weighted function for next step."""
        if self._experts:
            log_vals = list(self._experts.values())
            log_total_active = logsumexp(log_vals)
            expert_items = list(self._experts.items())
            log_total = log_total_active
            n_next = self._n + 1

            def _b_n(u, _items=expert_items, _log_total=log_total, _n=n_next,
                     _grid=self._grid, _M=self._batch_interval):
                val = 0.0
                for (bi, gi), log_cap in _items:
                    a, b = _grid[gi]
                    ratio = min((bi * _M) / _n, 1.0)
                    a_prime = ratio * a + (1 - ratio) * b
                    if 0 < a_prime < 1:
                        f = b / a_prime if u <= a_prime else (1 - b) / (1 - a_prime)
                    else:
                        f = 1.0
                    val += math.exp(log_cap - _log_total) * f
                return val

            def _B_n(u, _items=expert_items, _log_total=log_total, _n=n_next,
                     _grid=self._grid, _M=self._batch_interval):
                if u <= 0.0:
                    return 0.0
                if u >= 1.0:
                    return 1.0
                val = 0.0
                for (bi, gi), log_cap in _items:
                    a, b = _grid[gi]
                    ratio = min((bi * _M) / _n, 1.0)
                    a_prime = ratio * a + (1 - ratio) * b
                    w = math.exp(log_cap - _log_total)
                    if 0 < a_prime < 1:
                        if u <= a_prime:
                            val += w * (b / a_prime) * u
                        else:
                            val += w * (b + (1 - b) / (1 - a_prime) * (u - a_prime))
                    else:
                        val += w * u
                return val

            self.b_n = _b_n
            self.B_n = _B_n

update(p: float) -> None

Advance the Sleeper/Drifter martingale by one p-value.

Like :class:SleeperStayer, but each expert's threshold drifts over time, so the bet uses the current drifted threshold of every active expert before pooling their capital.

Parameters:

Name Type Description Default
p float

New p-value in \([0, 1]\).

 required
Source code in src/online_cp/martingale/sleepers.py 
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
def update(self, p: float) -> None:
    """Advance the Sleeper/Drifter martingale by one p-value.

    Like :class:`SleeperStayer`, but each expert's threshold *drifts* over
    time, so the bet uses the current drifted threshold of every active
    expert before pooling their capital.

    Parameters
    ----------
    p : float
        New p-value in $[0, 1]$.
    """
    if self.store_p_values:
        self.p_values.append(p)

    self._n += 1

    # Step 1: Bet — update all active experts in log-space
    keys_to_remove = []
    for key, log_capital in self._experts.items():
        batch_idx, grid_idx = key
        a_prime = self._get_drifted_threshold(batch_idx, grid_idx)
        _, b = self._grid[grid_idx]

        # Bet using f_{(a', b)}(p)
        if a_prime <= 0 or a_prime >= 1:
            log_bet_val = 0.0  # log(1) = no bet
        else:
            bet_val = b / a_prime if p <= a_prime else (1 - b) / (1 - a_prime)
            log_bet_val = math.log(bet_val)

        new_log_capital = log_capital + log_bet_val
        if new_log_capital < self._log_prune_threshold:
            keys_to_remove.append(key)
        else:
            self._experts[key] = new_log_capital

    for key in keys_to_remove:
        del self._experts[key]

    # Step 2: Output — total capital in log-space
    if self._experts:
        log_active_sum = logsumexp(list(self._experts.values()))
        self.logM = np.logaddexp(self._log_S_sleep, log_active_sum)
    else:
        self.logM = self._log_S_sleep
    self.log_martingale_values.append(self.logM)

    # Step 3: Wake new batch (if n is divisible by M) — prepares for next step
    if self._n % self._batch_interval == 0:
        batch_idx = self._n // self._batch_interval
        # log(transfer_per_expert) = log(R * M) + log_S_sleep - log(n_grid)
        log_transfer = (math.log(self.R * self._batch_interval)
                       + self._log_S_sleep
                       - math.log(self._n_grid))
        # S_sleep *= (1 - R*M), clamped to avoid negative
        rm = self.R * self._batch_interval
        if rm >= 1.0:
            self._log_S_sleep = -np.inf
        else:
            self._log_S_sleep += math.log(1.0 - rm)
        for j in range(self._n_grid):
            key = (batch_idx, j)
            if key in self._experts:
                self._experts[key] = np.logaddexp(self._experts[key], log_transfer)
            else:
                self._experts[key] = log_transfer

    # Mark b_n/B_n stale (lazy recomputation on access)
    self._mark_stale()

Legendre Jumper Martingales

online_cp.martingale.SimpleLegendreJumper

Bases: ConformalTestMartingale

Simple Legendre Jumper betting martingale (Algorithm 2).

Uses the betting function f_eps^(k)(p) = 1 + eps * P_k(2p-1) with a Markov chain over the state space E and jump rate J.

Parameters:

Name Type Description Default
order int

Degree k >= 1 of the shifted Legendre polynomial.

 1
J float

Jump rate in (0, 1].

 0.01
epsilon_grid tuple of float

State space E. Default is the paper's 5-point grid.

 STANDARD_GRID

Examples:

>>> slj = SimpleLegendreJumper(order=1, J=0.01)
>>> for _ in range(10):
...     slj.update(0.01)
>>> slj.M > 1.0
True
Source code in src/online_cp/martingale/legendre.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
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
class SimpleLegendreJumper(ConformalTestMartingale):
    """Simple Legendre Jumper betting martingale (Algorithm 2).

    Uses the betting function f_eps^(k)(p) = 1 + eps * P_k(2p-1)
    with a Markov chain over the state space E and jump rate J.

    Parameters
    ----------
    order : int
        Degree k >= 1 of the shifted Legendre polynomial.
    J : float
        Jump rate in (0, 1].
    epsilon_grid : tuple of float
        State space E. Default is the paper's 5-point grid.

    Examples
    --------
    >>> slj = SimpleLegendreJumper(order=1, J=0.01)
    >>> for _ in range(10):
    ...     slj.update(0.01)
    >>> slj.M > 1.0
    True
    """

    def __init__(self, order: int = 1, J: float = 0.01, epsilon_grid: tuple[float, ...] = STANDARD_GRID,
                 store_p_values: bool = True) -> None:
        super().__init__(store_p_values)
        if order < 1:
            raise ValueError("order must be >= 1")
        if not (0 < J <= 1):
            raise ValueError("J must be in (0, 1]")
        self.order = order
        self.J = J
        self.epsilon_grid = tuple(epsilon_grid)
        self.N = len(self.epsilon_grid)

        # Log-space state: log(C_eps) for each expert
        self._log_C = np.full(self.N, -np.log(self.N))  # uniform: 1/|E|
        self._log_total = 0.0  # log(C) = log(sum C_eps) = 0 initially

    def _P(self, u):
        """Evaluate shifted Legendre polynomial P_k(2u-1)."""
        return eval_legendre(self.order, 2.0 * u - 1.0)

    def update(self, p: float) -> None:
        r"""Advance the simple Legendre jumper by one p-value.

        Performs the jump (mixing over the $\epsilon$ grid) then bets with the
        shifted-Legendre density of the configured order, all in log-space
        ([LegendreJumper, preprint]).

        Parameters
        ----------
        p : float
            New p-value in $[0, 1]$.
        """
        if self.store_p_values:
            self.p_values.append(p)

        # --- Step 1: Jump ---
        # C_eps <- (1-J) * C_eps + (J / |E|) * C
        # In log-space: log((1-J)*C_eps + (J/N)*C)
        #             = logaddexp(log(1-J) + log_C_eps, log(J/N) + log_total)
        if self.J == 1.0:
            log_1_minus_J = -np.inf
        else:
            log_1_minus_J = np.log(1.0 - self.J)
        log_J_div_N = np.log(self.J / self.N)

        term1 = log_1_minus_J + self._log_C
        term2 = log_J_div_N + self._log_total
        self._log_C = np.logaddexp(term1, term2)

        # --- Step 2: Bet ---
        # C_eps <- C_eps * f_eps^(k)(p)
        # f_eps(p) = 1 + eps * P_k(p)
        P_p = self._P(p)
        bet_values = np.array([1.0 + eps * P_p for eps in self.epsilon_grid])
        self._log_C += np.log(bet_values)

        # --- Step 3: Accumulate ---
        # C = sum(C_eps)
        self._log_total = logsumexp(self._log_C)
        self.logM = self._log_total
        self.log_martingale_values.append(self.logM)

        # --- Mark b_n/B_n stale (lazy recomputation on access) ---
        self._mark_stale()

    def _update_exposed_functions(self):
        """Set b_n and B_n as the wealth-weighted effective betting function."""
        # Post-jump weights for next step
        if self.J == 1.0:
            log_1_minus_J = -np.inf
        else:
            log_1_minus_J = np.log(1.0 - self.J)
        log_J_div_N = np.log(self.J / self.N)

        term1 = log_1_minus_J + self._log_C
        term2 = log_J_div_N + self._log_total
        log_weights_next = np.logaddexp(term1, term2)
        log_sum = logsumexp(log_weights_next)
        weights = np.exp(log_weights_next - log_sum)

        # Effective epsilon: eps_bar = sum_i w_i * eps_i
        eps_bar = np.dot(weights, self.epsilon_grid)
        order = self.order

        def _b_n(u, _eps=eps_bar, _k=order):
            return 1.0 + _eps * eval_legendre(_k, 2.0 * u - 1.0)

        def _B_n(u, _eps=eps_bar, _k=order):
            if u <= 0:
                return 0.0
            if u >= 1:
                return 1.0
            # B_n(u) = u + eps * Q_k(u) where Q_k = integral_0^u P_k(t) dt
            P_kp1 = eval_legendre(_k + 1, 2.0 * u - 1.0)
            P_km1 = eval_legendre(_k - 1, 2.0 * u - 1.0)
            Q = (P_kp1 - P_km1) / (2.0 * (2 * _k + 1))
            return u + _eps * Q

        self.b_n = _b_n
        self.B_n = _B_n

        # Closed-form inverse for k=1 (quadratic); otherwise numerical fallback
        if order == 1:
            a = eps_bar

            def _B_n_inv(x, _a=a):
                if x <= 0.0:
                    return 0.0
                if x >= 1.0:
                    return 1.0
                if abs(_a) < 1e-12:
                    return x
                return (-(1 - _a) + np.sqrt((1 - _a) ** 2 + 4 * _a * x)) / (2 * _a)

            self.B_n_inv = _B_n_inv

update(p: float) -> None

Advance the simple Legendre jumper by one p-value.

Performs the jump (mixing over the \(\epsilon\) grid) then bets with the shifted-Legendre density of the configured order, all in log-space ([LegendreJumper, preprint]).

Parameters:

Name Type Description Default
p float

New p-value in \([0, 1]\).

 required
Source code in src/online_cp/martingale/legendre.py 
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
def update(self, p: float) -> None:
    r"""Advance the simple Legendre jumper by one p-value.

    Performs the jump (mixing over the $\epsilon$ grid) then bets with the
    shifted-Legendre density of the configured order, all in log-space
    ([LegendreJumper, preprint]).

    Parameters
    ----------
    p : float
        New p-value in $[0, 1]$.
    """
    if self.store_p_values:
        self.p_values.append(p)

    # --- Step 1: Jump ---
    # C_eps <- (1-J) * C_eps + (J / |E|) * C
    # In log-space: log((1-J)*C_eps + (J/N)*C)
    #             = logaddexp(log(1-J) + log_C_eps, log(J/N) + log_total)
    if self.J == 1.0:
        log_1_minus_J = -np.inf
    else:
        log_1_minus_J = np.log(1.0 - self.J)
    log_J_div_N = np.log(self.J / self.N)

    term1 = log_1_minus_J + self._log_C
    term2 = log_J_div_N + self._log_total
    self._log_C = np.logaddexp(term1, term2)

    # --- Step 2: Bet ---
    # C_eps <- C_eps * f_eps^(k)(p)
    # f_eps(p) = 1 + eps * P_k(p)
    P_p = self._P(p)
    bet_values = np.array([1.0 + eps * P_p for eps in self.epsilon_grid])
    self._log_C += np.log(bet_values)

    # --- Step 3: Accumulate ---
    # C = sum(C_eps)
    self._log_total = logsumexp(self._log_C)
    self.logM = self._log_total
    self.log_martingale_values.append(self.logM)

    # --- Mark b_n/B_n stale (lazy recomputation on access) ---
    self._mark_stale()

online_cp.martingale.ProductLegendreJumper

Bases: ConformalTestMartingale

Product Legendre Jumper betting martingale (Algorithm 3).

Maintains a single Markov chain over the full Cartesian product state space E = E_1 x E_2 x ... x E_K, with betting function:

f_eps^K(p) = prod_k (1 + eps_k * P_k(p)) / Z(eps)

Parameters:

Name Type Description Default
orders list of int

Set K of Legendre polynomial degrees (each >= 1).

 None
J float

Jump rate in (0, 1].

 0.01
epsilon_grid tuple of float

Per-order state space E_k. Default is the paper's 5-point grid.

 STANDARD_GRID

Examples:

>>> plj = ProductLegendreJumper(orders=[1, 2], J=0.01)
>>> for _ in range(10):
...     plj.update(0.01)
>>> plj.M > 1.0
True
Source code in src/online_cp/martingale/legendre.py 
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
class ProductLegendreJumper(ConformalTestMartingale):
    """Product Legendre Jumper betting martingale (Algorithm 3).

    Maintains a single Markov chain over the full Cartesian product state
    space E = E_1 x E_2 x ... x E_K, with betting function:

        f_eps^K(p) = prod_k (1 + eps_k * P_k(p)) / Z(eps)

    Parameters
    ----------
    orders : list of int
        Set K of Legendre polynomial degrees (each >= 1).
    J : float
        Jump rate in (0, 1].
    epsilon_grid : tuple of float
        Per-order state space E_k. Default is the paper's 5-point grid.

    Examples
    --------
    >>> plj = ProductLegendreJumper(orders=[1, 2], J=0.01)
    >>> for _ in range(10):
    ...     plj.update(0.01)
    >>> plj.M > 1.0
    True
    """

    def __init__(self, orders: list[int] | None = None, J: float = 0.01, epsilon_grid: tuple[float, ...] = STANDARD_GRID,
                 store_p_values: bool = True) -> None:
        super().__init__(store_p_values)
        if orders is None:
            orders = [1, 2, 3]
        if not orders:
            raise ValueError("At least one order is required.")
        if any(k < 1 for k in orders):
            raise ValueError("All orders must be >= 1.")
        if not (0 < J <= 1):
            raise ValueError("J must be in (0, 1]")
        self.orders = list(orders)
        self.K = len(self.orders)
        self.J = J
        self.epsilon_grid = tuple(epsilon_grid)

        # Build full Cartesian product state space
        self.states = list(itertools.product(self.epsilon_grid, repeat=self.K))
        self.N = len(self.states)
        self._states_array = np.array(self.states)  # shape (N, K) for vectorized _eval_bets

        if self.N > 500:
            import warnings as _w
            _w.warn(
                f"Product state space has {self.N} experts "
                f"(orders={self.orders}, grid size={len(self.epsilon_grid)}). "
                f"Consider using fewer orders or a coarser grid.",
                stacklevel=2,
            )

        # Pre-compute Z(eps) for every state (done once at init)
        self._log_Z = np.zeros(self.N)
        for i, eps_vec in enumerate(self.states):
            Z = compute_normalization_Z(self.orders, eps_vec)
            self._log_Z[i] = np.log(Z)

        # Log-space state: uniform prior over all experts
        self._log_C = np.full(self.N, -np.log(self.N))
        self._log_total = 0.0

    def _eval_bets(self, p):
        """Evaluate f_eps^K(p) for all states. Returns array of shape (N,)."""
        # Compute P_k(p) for each order once
        Pk_vals = np.array([eval_legendre(k, 2.0 * p - 1.0) for k in self.orders])

        # Vectorized: (N, K) * (K,) -> (N, K), then sum over K axis
        return np.log1p(self._states_array * Pk_vals).sum(axis=1) - self._log_Z

    def update(self, p: float) -> None:
        r"""Advance the product Legendre jumper by one p-value.

        Bets with a *product* of shifted-Legendre densities over several orders,
        after the jump-mixing step over the $\epsilon$ grid
        ([LegendreJumper, preprint]).

        Parameters
        ----------
        p : float
            New p-value in $[0, 1]$.
        """
        if self.store_p_values:
            self.p_values.append(p)

        # --- Step 1: Jump ---
        # C_eps <- (1-J) * C_eps + (J / |E|) * C
        if self.J == 1.0:
            log_1_minus_J = -np.inf
        else:
            log_1_minus_J = np.log(1.0 - self.J)
        log_J_div_N = np.log(self.J / self.N)

        term1 = log_1_minus_J + self._log_C
        term2 = log_J_div_N + self._log_total
        self._log_C = np.logaddexp(term1, term2)

        # --- Step 2: Bet ---
        # C_eps <- C_eps * f_eps^K(p)
        log_bets = self._eval_bets(p)
        self._log_C += log_bets

        # --- Step 3: Accumulate ---
        self._log_total = logsumexp(self._log_C)
        self.logM = self._log_total
        self.log_martingale_values.append(self.logM)

        # --- Mark b_n/B_n stale (lazy recomputation on access) ---
        self._mark_stale()

    def _update_exposed_functions(self):
        """Set b_n and B_n as the wealth-weighted mixture of expert betting functions."""
        # Post-jump weights for next step
        if self.J == 1.0:
            log_1_minus_J = -np.inf
        else:
            log_1_minus_J = np.log(1.0 - self.J)
        log_J_div_N = np.log(self.J / self.N)

        term1 = log_1_minus_J + self._log_C
        term2 = log_J_div_N + self._log_total
        log_weights_next = np.logaddexp(term1, term2)
        log_sum = logsumexp(log_weights_next)
        weights = np.exp(log_weights_next - log_sum)

        # Capture state for closures
        orders = self.orders
        states = self.states
        log_Z = self._log_Z

        def _b_n(u, _w=weights, _orders=orders, _states=states, _log_Z=log_Z):
            Pk_vals = [eval_legendre(k, 2.0 * u - 1.0) for k in _orders]
            total = 0.0
            for i, eps_vec in enumerate(_states):
                prod = 1.0
                for j, eps in enumerate(eps_vec):
                    prod *= (1.0 + eps * Pk_vals[j])
                total += _w[i] * prod / np.exp(_log_Z[i])
            return total

        def _B_n(u, _b_n_func=None):
            if u <= 0:
                return 0.0
            if u >= 1:
                return 1.0
            # Numerical integration (the mixture doesn't simplify to a closed-form CDF)
            from scipy.integrate import quad
            val, _ = quad(_b_n, 1e-12, u, limit=50)
            return val

        self.b_n = _b_n
        self.B_n = _B_n

update(p: float) -> None

Advance the product Legendre jumper by one p-value.

Bets with a product of shifted-Legendre densities over several orders, after the jump-mixing step over the \(\epsilon\) grid ([LegendreJumper, preprint]).

Parameters:

Name Type Description Default
p float

New p-value in \([0, 1]\).

 required
Source code in src/online_cp/martingale/legendre.py 
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
def update(self, p: float) -> None:
    r"""Advance the product Legendre jumper by one p-value.

    Bets with a *product* of shifted-Legendre densities over several orders,
    after the jump-mixing step over the $\epsilon$ grid
    ([LegendreJumper, preprint]).

    Parameters
    ----------
    p : float
        New p-value in $[0, 1]$.
    """
    if self.store_p_values:
        self.p_values.append(p)

    # --- Step 1: Jump ---
    # C_eps <- (1-J) * C_eps + (J / |E|) * C
    if self.J == 1.0:
        log_1_minus_J = -np.inf
    else:
        log_1_minus_J = np.log(1.0 - self.J)
    log_J_div_N = np.log(self.J / self.N)

    term1 = log_1_minus_J + self._log_C
    term2 = log_J_div_N + self._log_total
    self._log_C = np.logaddexp(term1, term2)

    # --- Step 2: Bet ---
    # C_eps <- C_eps * f_eps^K(p)
    log_bets = self._eval_bets(p)
    self._log_C += log_bets

    # --- Step 3: Accumulate ---
    self._log_total = logsumexp(self._log_C)
    self.logM = self._log_total
    self.log_martingale_values.append(self.logM)

    # --- Mark b_n/B_n stale (lazy recomputation on access) ---
    self._mark_stale()

online_cp.martingale.VariationalLegendreJumper

Bases: ConformalTestMartingale

Variational Legendre Jumper betting martingale (Algorithm 4).

Runs |K| independent sub-jumpers (each an instance of the SLJ logic), one per polynomial degree. At each step, consensus parameters are computed as the wealth-weighted mean epsilon from each sub-jumper, and the global martingale bets with the Z-normalised product betting function evaluated at these consensus parameters.

Computational cost: O(|K| * g) per step (linear, not exponential).

Parameters:

Name Type Description Default
orders list of int

Set K of Legendre polynomial degrees (each >= 1).

 None
J float

Jump rate in (0, 1].

 0.01
epsilon_grid tuple of float

Per-order state space E_k. Default is the paper's 5-point grid.

 STANDARD_GRID

Examples:

>>> vlj = VariationalLegendreJumper(orders=[1, 2], J=0.01)
>>> for _ in range(10):
...     vlj.update(0.01)
>>> vlj.M > 1.0
True
Source code in src/online_cp/martingale/legendre.py 
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
class VariationalLegendreJumper(ConformalTestMartingale):
    """Variational Legendre Jumper betting martingale (Algorithm 4).

    Runs |K| independent sub-jumpers (each an instance of the SLJ logic),
    one per polynomial degree. At each step, consensus parameters are
    computed as the wealth-weighted mean epsilon from each sub-jumper,
    and the global martingale bets with the Z-normalised product betting
    function evaluated at these consensus parameters.

    Computational cost: O(|K| * g) per step (linear, not exponential).

    Parameters
    ----------
    orders : list of int
        Set K of Legendre polynomial degrees (each >= 1).
    J : float
        Jump rate in (0, 1].
    epsilon_grid : tuple of float
        Per-order state space E_k. Default is the paper's 5-point grid.

    Examples
    --------
    >>> vlj = VariationalLegendreJumper(orders=[1, 2], J=0.01)
    >>> for _ in range(10):
    ...     vlj.update(0.01)
    >>> vlj.M > 1.0
    True
    """

    def __init__(self, orders: list[int] | None = None, J: float = 0.01, epsilon_grid: tuple[float, ...] = STANDARD_GRID,
                 store_p_values: bool = True) -> None:
        super().__init__(store_p_values)
        if orders is None:
            orders = [1, 2, 3]
        if not orders:
            raise ValueError("At least one order is required.")
        if any(k < 1 for k in orders):
            raise ValueError("All orders must be >= 1.")
        if not (0 < J <= 1):
            raise ValueError("J must be in (0, 1]")
        self.orders = list(orders)
        self.K = len(self.orders)
        self.J = J
        self.epsilon_grid = tuple(epsilon_grid)
        self.g = len(self.epsilon_grid)
        self._eps_array = np.array(self.epsilon_grid)

        # Log-space sub-jumper state: log(C_{k, eps})
        # Shape: (K, g). Initialised to log(1/g) = -log(g).
        self._log_C = np.full((self.K, self.g), -np.log(self.g))
        # log(C_k) = log(sum_eps C_{k, eps}) = 0 initially
        self._log_C_k = np.zeros(self.K)

        # Global martingale value in log-space
        self._log_S = 0.0

        # Pre-compute log constants for the jump step
        if self.J == 1.0:
            self._log_1_minus_J = -np.inf
        else:
            self._log_1_minus_J = np.log(1.0 - self.J)
        self._log_J_div_g = np.log(self.J / self.g)

        # Pre-compute Gaunt coefficients for fast Z evaluation.
        # Z(eps) = 1 + sum_{|S|>=3} (prod_{k in S} eps_k) * G_S
        # where G_S = integral_0^1 prod_{k in S} P_k(p) dp.
        # For |K| <= 2, Z = 1 always (by orthogonality).
        self._gaunt_terms = []  # list of (index_tuple, G_S)
        if self.K >= 3:
            for size in range(3, self.K + 1):
                for indices in itertools.combinations(range(self.K), size):
                    poly = Polynomial([1.0])
                    for idx in indices:
                        poly = poly * shifted_legendre_poly(self.orders[idx])
                    antideriv = poly.integ()
                    G_S = float(antideriv(1.0) - antideriv(0.0))
                    if abs(G_S) > 1e-15:
                        self._gaunt_terms.append((indices, G_S))

    def _fast_Z(self, eps_bar):
        """Compute Z(eps_bar) using precomputed Gaunt coefficients. O(1) for typical |K|."""
        if not self._gaunt_terms:
            return 1.0
        Z = 1.0
        for indices, G_S in self._gaunt_terms:
            prod = 1.0
            for idx in indices:
                prod *= eps_bar[idx]
            Z += prod * G_S
        return Z

    def update(self, p: float) -> None:
        r"""Advance the variational Legendre jumper by one p-value.

        Maintains per-sub-jumper consensus parameters $\bar\epsilon_k$ via a
        variational update before betting ([LegendreJumper, preprint]).

        Parameters
        ----------
        p : float
            New p-value in $[0, 1]$.
        """
        if self.store_p_values:
            self.p_values.append(p)

        # --- Step 1: Jump (per sub-jumper) ---
        # C_{k,eps} <- (1-J) * C_{k,eps} + (J/g) * C_k
        # In log-space: logaddexp(log(1-J) + log_C, log(J/g) + log_C_k)
        term1 = self._log_1_minus_J + self._log_C
        term2 = self._log_J_div_g + self._log_C_k[:, np.newaxis]
        self._log_C = np.logaddexp(term1, term2)

        # --- Step 2: Consensus parameters ---
        # eps_bar_k = sum_eps eps * (C_{k,eps} / C_k)
        # weights = exp(log_C - log_C_k) are normalized probabilities
        weights = np.exp(self._log_C - self._log_C_k[:, np.newaxis])
        eps_bar = (weights * self._eps_array[np.newaxis, :]).sum(axis=1)

        # --- Step 3: Global bet ---
        # S_n = S_{n-1} * prod_k(1 + eps_bar_k * P_k(p)) / Z(eps_bar)
        Pk_vals = np.array([eval_legendre(k, 2.0 * p - 1.0) for k in self.orders])
        Z = self._fast_Z(eps_bar)
        log_global_bet = np.sum(np.log(1.0 + eps_bar * Pk_vals)) - np.log(Z)
        self._log_S += log_global_bet

        # --- Step 4: Update sub-jumpers ---
        # C_{k,eps} <- C_{k,eps} * (1 + eps * P_k(p))
        marginal_bets = 1.0 + self._eps_array[np.newaxis, :] * Pk_vals[:, np.newaxis]
        self._log_C += np.log(marginal_bets)

        # --- Step 5: Recompute log_C_k ---
        self._log_C_k = logsumexp(self._log_C, axis=1)

        # --- Update martingale ---
        self.logM = self._log_S
        self.log_martingale_values.append(self.logM)

        # --- Mark b_n/B_n stale (lazy recomputation on access) ---
        self._mark_stale()

    def _update_exposed_functions(self):
        """Set b_n and B_n using consensus parameters for the next step."""
        # Simulate post-jump state (what would happen next step)
        term1 = self._log_1_minus_J + self._log_C
        term2 = self._log_J_div_g + self._log_C_k[:, np.newaxis]
        log_C_next = np.logaddexp(term1, term2)
        log_C_k_next = logsumexp(log_C_next, axis=1)

        weights = np.exp(log_C_next - log_C_k_next[:, np.newaxis])
        eps_bar = (weights * self._eps_array[np.newaxis, :]).sum(axis=1)

        orders = self.orders
        gaunt_terms = self._gaunt_terms

        def _b_n(u, _eps_bar=eps_bar.copy(), _orders=orders, _gaunt=gaunt_terms):  # noqa: B008
            Pk_vals = [eval_legendre(k, 2.0 * u - 1.0) for k in _orders]
            # Fast Z from precomputed Gaunt coefficients
            Z = 1.0
            for indices, G_S in _gaunt:
                prod = 1.0
                for idx in indices:
                    prod *= _eps_bar[idx]
                Z += prod * G_S
            val = 1.0
            for j, eps in enumerate(_eps_bar):
                val *= (1.0 + eps * Pk_vals[j])
            return val / Z

        def _B_n(u):
            if u <= 0:
                return 0.0
            if u >= 1:
                return 1.0
            from scipy.integrate import quad
            val, _ = quad(_b_n, 1e-12, u, limit=50)
            return val

        self.b_n = _b_n
        self.B_n = _B_n

update(p: float) -> None

Advance the variational Legendre jumper by one p-value.

Maintains per-sub-jumper consensus parameters \(\bar\epsilon_k\) via a variational update before betting ([LegendreJumper, preprint]).

Parameters:

Name Type Description Default
p float

New p-value in \([0, 1]\).

 required
Source code in src/online_cp/martingale/legendre.py 
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
def update(self, p: float) -> None:
    r"""Advance the variational Legendre jumper by one p-value.

    Maintains per-sub-jumper consensus parameters $\bar\epsilon_k$ via a
    variational update before betting ([LegendreJumper, preprint]).

    Parameters
    ----------
    p : float
        New p-value in $[0, 1]$.
    """
    if self.store_p_values:
        self.p_values.append(p)

    # --- Step 1: Jump (per sub-jumper) ---
    # C_{k,eps} <- (1-J) * C_{k,eps} + (J/g) * C_k
    # In log-space: logaddexp(log(1-J) + log_C, log(J/g) + log_C_k)
    term1 = self._log_1_minus_J + self._log_C
    term2 = self._log_J_div_g + self._log_C_k[:, np.newaxis]
    self._log_C = np.logaddexp(term1, term2)

    # --- Step 2: Consensus parameters ---
    # eps_bar_k = sum_eps eps * (C_{k,eps} / C_k)
    # weights = exp(log_C - log_C_k) are normalized probabilities
    weights = np.exp(self._log_C - self._log_C_k[:, np.newaxis])
    eps_bar = (weights * self._eps_array[np.newaxis, :]).sum(axis=1)

    # --- Step 3: Global bet ---
    # S_n = S_{n-1} * prod_k(1 + eps_bar_k * P_k(p)) / Z(eps_bar)
    Pk_vals = np.array([eval_legendre(k, 2.0 * p - 1.0) for k in self.orders])
    Z = self._fast_Z(eps_bar)
    log_global_bet = np.sum(np.log(1.0 + eps_bar * Pk_vals)) - np.log(Z)
    self._log_S += log_global_bet

    # --- Step 4: Update sub-jumpers ---
    # C_{k,eps} <- C_{k,eps} * (1 + eps * P_k(p))
    marginal_bets = 1.0 + self._eps_array[np.newaxis, :] * Pk_vals[:, np.newaxis]
    self._log_C += np.log(marginal_bets)

    # --- Step 5: Recompute log_C_k ---
    self._log_C_k = logsumexp(self._log_C, axis=1)

    # --- Update martingale ---
    self.logM = self._log_S
    self.log_martingale_values.append(self.logM)

    # --- Mark b_n/B_n stale (lazy recomputation on access) ---
    self._mark_stale()

online_cp.martingale.CompositeLegendreJumper

Bases: ConformalTestMartingale

Composite Legendre Jumper that averages over multiple jump rates.

Creates multiple instances of a base Legendre Jumper class (one per jumping rate) and computes the martingale as their arithmetic mean. This is the direct analogue of the Composite Jumper for Legendre martingales.

Parameters:

Name Type Description Default
base_class class

The base martingale class to instantiate. Must accept a J parameter. Default is SimpleLegendreJumper.

 None
J list of float or None

List of jumping rates. Default is [1e-4, 1e-3, 1e-2, 1e-1, 1.0].

 None
**kwargs Any

Additional keyword arguments forwarded to base_class (e.g. order, orders, epsilon_grid).

 {}

Examples:

>>> clj = CompositeLegendreJumper()
>>> for _ in range(5):
...     clj.update(0.01)
>>> clj.M > 1.0
True

>>> from online_cp.martingale import VariationalLegendreJumper
>>> clj = CompositeLegendreJumper(
...     base_class=VariationalLegendreJumper, orders=[1, 2]
... )
>>> for _ in range(5):
...     clj.update(0.01)
>>> clj.M > 1.0
True
Source code in src/online_cp/martingale/legendre.py 
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
class CompositeLegendreJumper(ConformalTestMartingale):
    """Composite Legendre Jumper that averages over multiple jump rates.

    Creates multiple instances of a base Legendre Jumper class (one per
    jumping rate) and computes the martingale as their arithmetic mean.
    This is the direct analogue of the Composite Jumper for Legendre
    martingales.

    Parameters
    ----------
    base_class : class
        The base martingale class to instantiate. Must accept a ``J``
        parameter. Default is ``SimpleLegendreJumper``.
    J : list of float or None
        List of jumping rates. Default is [1e-4, 1e-3, 1e-2, 1e-1, 1.0].
    **kwargs
        Additional keyword arguments forwarded to ``base_class`` (e.g.
        ``order``, ``orders``, ``epsilon_grid``).

    Examples
    --------
    >>> clj = CompositeLegendreJumper()
    >>> for _ in range(5):
    ...     clj.update(0.01)
    >>> clj.M > 1.0
    True

    >>> from online_cp.martingale import VariationalLegendreJumper
    >>> clj = CompositeLegendreJumper(
    ...     base_class=VariationalLegendreJumper, orders=[1, 2]
    ... )
    >>> for _ in range(5):
    ...     clj.update(0.01)
    >>> clj.M > 1.0
    True
    """

    def __init__(self, base_class: type[ConformalTestMartingale] | None = None, J: list[float] | None = None,
                 store_p_values: bool = True, **kwargs: Any) -> None:
        super().__init__(store_p_values)
        if base_class is None:
            base_class = SimpleLegendreJumper
        if J is None:
            J = [1e-4, 1e-3, 1e-2, 1e-1, 1.0]

        self.J = list(J)
        self.base_class = base_class
        self._jumpers = [
            base_class(J=j, store_p_values=False, **kwargs)
            for j in self.J
        ]
        self._n_jumpers = len(self._jumpers)
        self._log_n = np.log(self._n_jumpers)

    def update(self, p: float) -> None:
        """Advance every component Legendre jumper and pool them.

        Updates each sub-jumper on ``p`` and sets the composite log-martingale to
        the equal-weight log-mean of the components ([LegendreJumper, preprint]).

        Parameters
        ----------
        p : float
            New p-value in $[0, 1]$.
        """
        if self.store_p_values:
            self.p_values.append(p)

        for m in self._jumpers:
            m.update(p)

        log_Ms = np.array([m.logM for m in self._jumpers])
        self.logM = logsumexp(log_Ms) - self._log_n

        self.log_martingale_values.append(self.logM)

        self._mark_stale()

update(p: float) -> None

Advance every component Legendre jumper and pool them.

Updates each sub-jumper on p and sets the composite log-martingale to the equal-weight log-mean of the components ([LegendreJumper, preprint]).

Parameters:

Name Type Description Default
p float

New p-value in \([0, 1]\).

 required
Source code in src/online_cp/martingale/legendre.py 
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
def update(self, p: float) -> None:
    """Advance every component Legendre jumper and pool them.

    Updates each sub-jumper on ``p`` and sets the composite log-martingale to
    the equal-weight log-mean of the components ([LegendreJumper, preprint]).

    Parameters
    ----------
    p : float
        New p-value in $[0, 1]$.
    """
    if self.store_p_values:
        self.p_values.append(p)

    for m in self._jumpers:
        m.update(p)

    log_Ms = np.array([m.logM for m in self._jumpers])
    self.logM = logsumexp(log_Ms) - self._log_n

    self.log_martingale_values.append(self.logM)

    self._mark_stale()

Change-Point Detection Wrappers

online_cp.martingale.VilleWrapper

Ville's inequality procedure for change-point detection.

The simplest test based on a conformal test martingale: reject the exchangeability hypothesis when the running maximum of the martingale exceeds a threshold c. By Ville's inequality:

P(∃n : S_n >= c) <= 1/c

So threshold c = 20 gives a 5% significance level, c = 100 gives 1%, etc.

Parameters:

Name Type Description Default
martingale ConformalTestMartingale

The underlying martingale to wrap.

 required
threshold float

Default alarm threshold (default 20, i.e. 5% significance).

 20
References

Vovk, Gammerman & Shafer (2022). Algorithmic Learning in a Random World, 2nd edition, §8.4.1 (The Ville Procedure). Cambridge University Press.

Examples:

>>> from online_cp.martingale import SimpleJumper, VilleWrapper
>>> sj = SimpleJumper(J=0.1)
>>> ville = VilleWrapper(sj, threshold=20)
>>> for _ in range(10):
...     ville.update(0.5)
>>> bool(ville.rejected)
False
Source code in src/online_cp/martingale/wrappers.py 
 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
class VilleWrapper:
    """Ville's inequality procedure for change-point detection.

    The simplest test based on a conformal test martingale: reject the
    exchangeability hypothesis when the running maximum of the martingale
    exceeds a threshold c. By Ville's inequality:

        P(∃n : S_n >= c) <= 1/c

    So threshold c = 20 gives a 5% significance level, c = 100 gives 1%, etc.

    Parameters
    ----------
    martingale : ConformalTestMartingale
        The underlying martingale to wrap.
    threshold : float
        Default alarm threshold (default 20, i.e. 5% significance).

    References
    ----------
    Vovk, Gammerman & Shafer (2022). *Algorithmic Learning in a Random World*,
    2nd edition, §8.4.1 (The Ville Procedure). Cambridge University Press.

    Examples
    --------
    >>> from online_cp.martingale import SimpleJumper, VilleWrapper
    >>> sj = SimpleJumper(J=0.1)
    >>> ville = VilleWrapper(sj, threshold=20)
    >>> for _ in range(10):
    ...     ville.update(0.5)
    >>> bool(ville.rejected)
    False
    """

    def __init__(self, martingale, threshold=20):
        self.martingale = martingale
        self.threshold = threshold
        self._log_max = 0.0
        self._n = 0
        self._rejection_time = None

    @property
    def log_max(self):
        """Log of the running maximum of the martingale."""
        return self._log_max

    @property
    def max(self):
        """Running maximum of the martingale."""
        return np.exp(self._log_max)

    @property
    def rejected(self):
        """Whether the exchangeability hypothesis has been rejected."""
        return self._log_max >= np.log(self.threshold)

    @property
    def rejection_time(self):
        """Step at which the hypothesis was first rejected, or None."""
        return self._rejection_time

    def update(self, p: float) -> None:
        """Update the inner martingale and track the running maximum."""
        self.martingale.update(p)
        self._n += 1
        if self.martingale.logM > self._log_max:
            self._log_max = self.martingale.logM
        if self._rejection_time is None and self._log_max >= np.log(self.threshold):
            self._rejection_time = self._n

    def alarm(self, threshold=None):
        """Check whether max(S_n) exceeds the threshold.

        Parameters
        ----------
        threshold : float or None
            Override threshold. If None, uses the threshold set at construction.

        Returns
        -------
        bool
            True if the running maximum exceeds the threshold.
        """
        t = threshold if threshold is not None else self.threshold
        return self._log_max >= np.log(t)

log_max property

Log of the running maximum of the martingale.

max property

Running maximum of the martingale.

rejected property

Whether the exchangeability hypothesis has been rejected.

rejection_time property

Step at which the hypothesis was first rejected, or None.

update(p: float) -> None

Update the inner martingale and track the running maximum.

Source code in src/online_cp/martingale/wrappers.py 
87
88
89
90
91
92
93
94
def update(self, p: float) -> None:
    """Update the inner martingale and track the running maximum."""
    self.martingale.update(p)
    self._n += 1
    if self.martingale.logM > self._log_max:
        self._log_max = self.martingale.logM
    if self._rejection_time is None and self._log_max >= np.log(self.threshold):
        self._rejection_time = self._n

alarm(threshold=None)

Check whether max(S_n) exceeds the threshold.

Parameters:

Name Type Description Default
threshold float or None

Override threshold. If None, uses the threshold set at construction.

 None

Returns:

Type Description
bool

True if the running maximum exceeds the threshold.
Source code in src/online_cp/martingale/wrappers.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
def alarm(self, threshold=None):
    """Check whether max(S_n) exceeds the threshold.

    Parameters
    ----------
    threshold : float or None
        Override threshold. If None, uses the threshold set at construction.

    Returns
    -------
    bool
        True if the running maximum exceeds the threshold.
    """
    t = threshold if threshold is not None else self.threshold
    return self._log_max >= np.log(t)

online_cp.martingale.CUSUMWrapper

CUSUM change-detection wrapper for any conformal test martingale.

Computes the Page CUSUM statistic as the ratio of the current martingale value to its running minimum:

gamma_n = S_n / min_{i <= n} S_i

In log-space: log(gamma_n) = logM_n - min_{i <= n} logM_i

This removes any accumulated "debt" from an initial in-control period, giving faster detection after the change-point. Optionally accepts a linear barrier for controlling the false alarm rate over long horizons.

Parameters:

Name Type Description Default
martingale ConformalTestMartingale

The underlying martingale to wrap.

 required
barrier_slope float or None

If not None, the alarm threshold grows linearly as barrier_slope * n.

 None
References

Vovk, Gammerman & Shafer (2022). Algorithmic Learning in a Random World, 2nd edition, §8.3. Cambridge University Press.

Examples:

>>> from online_cp.martingale import SimpleJumper, CUSUMWrapper
>>> sj = SimpleJumper(J=0.01)
>>> cusum = CUSUMWrapper(sj)
>>> for _ in range(10):
...     cusum.update(0.5)
>>> bool(cusum.gamma >= 1.0)  # gamma is always >= 1 (since S_n >= min S_i is not guaranteed)
True
Source code in src/online_cp/martingale/wrappers.py 
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
class CUSUMWrapper:
    """CUSUM change-detection wrapper for any conformal test martingale.

    Computes the Page CUSUM statistic as the ratio of the current martingale
    value to its running minimum:

        gamma_n = S_n / min_{i <= n} S_i

    In log-space: log(gamma_n) = logM_n - min_{i <= n} logM_i

    This removes any accumulated "debt" from an initial in-control period,
    giving faster detection after the change-point. Optionally accepts a linear
    barrier for controlling the false alarm rate over long horizons.

    Parameters
    ----------
    martingale : ConformalTestMartingale
        The underlying martingale to wrap.
    barrier_slope : float or None
        If not None, the alarm threshold grows linearly as barrier_slope * n.

    References
    ----------
    Vovk, Gammerman & Shafer (2022). *Algorithmic Learning in a Random World*,
    2nd edition, §8.3. Cambridge University Press.

    Examples
    --------
    >>> from online_cp.martingale import SimpleJumper, CUSUMWrapper
    >>> sj = SimpleJumper(J=0.01)
    >>> cusum = CUSUMWrapper(sj)
    >>> for _ in range(10):
    ...     cusum.update(0.5)
    >>> bool(cusum.gamma >= 1.0)  # gamma is always >= 1 (since S_n >= min S_i is not guaranteed)
    True
    """

    def __init__(self, martingale, barrier_slope=None):
        self.martingale = martingale
        self.barrier_slope = barrier_slope
        self._log_min = 0.0  # log(S_0) = 0
        self._n = 0
        self._log_gamma_values = [0.0]

    @property
    def gamma(self):
        """Current CUSUM statistic."""
        return np.exp(self._log_gamma_values[-1])

    @property
    def log_gamma(self):
        """Current log CUSUM statistic."""
        return self._log_gamma_values[-1]

    @property
    def cusum_values(self):
        """All CUSUM statistic values."""
        return np.exp(self._log_gamma_values)

    @property
    def log_cusum_values(self):
        """All log CUSUM statistic values."""
        return self._log_gamma_values

    def update(self, p: float) -> None:
        """Update the inner martingale and recompute the CUSUM statistic."""
        self.martingale.update(p)
        self._n += 1

        logM = self.martingale.logM
        if logM < self._log_min:
            self._log_min = logM

        log_gamma = logM - self._log_min
        self._log_gamma_values.append(log_gamma)

    def alarm(self, threshold):
        """Check whether gamma_n exceeds the threshold (optionally with barrier).

        Parameters
        ----------
        threshold : float
            The alarm threshold. If barrier_slope is set, the effective
            threshold at step n is threshold + barrier_slope * n.

        Returns
        -------
        bool
            True if gamma_n exceeds the (possibly time-varying) threshold.
        """
        effective = threshold
        if self.barrier_slope is not None:
            effective = threshold + self.barrier_slope * self._n
        return self.gamma > effective

gamma property

Current CUSUM statistic.

log_gamma property

Current log CUSUM statistic.

cusum_values property

All CUSUM statistic values.

log_cusum_values property

All log CUSUM statistic values.

update(p: float) -> None

Update the inner martingale and recompute the CUSUM statistic.

Source code in src/online_cp/martingale/wrappers.py 
177
178
179
180
181
182
183
184
185
186
187
def update(self, p: float) -> None:
    """Update the inner martingale and recompute the CUSUM statistic."""
    self.martingale.update(p)
    self._n += 1

    logM = self.martingale.logM
    if logM < self._log_min:
        self._log_min = logM

    log_gamma = logM - self._log_min
    self._log_gamma_values.append(log_gamma)

alarm(threshold)

Check whether gamma_n exceeds the threshold (optionally with barrier).

Parameters:

Name Type Description Default
threshold float

The alarm threshold. If barrier_slope is set, the effective threshold at step n is threshold + barrier_slope * n.

 required

Returns:

Type Description
bool

True if gamma_n exceeds the (possibly time-varying) threshold.
Source code in src/online_cp/martingale/wrappers.py 
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
def alarm(self, threshold):
    """Check whether gamma_n exceeds the threshold (optionally with barrier).

    Parameters
    ----------
    threshold : float
        The alarm threshold. If barrier_slope is set, the effective
        threshold at step n is threshold + barrier_slope * n.

    Returns
    -------
    bool
        True if gamma_n exceeds the (possibly time-varying) threshold.
    """
    effective = threshold
    if self.barrier_slope is not None:
        effective = threshold + self.barrier_slope * self._n
    return self.gamma > effective

online_cp.martingale.ShiryaevRobertsWrapper

Shiryaev-Roberts change-detection wrapper for any conformal test martingale.

Computes the Shiryaev-Roberts statistic as:

R_n = sum_{i=1}^{n} S_n / S_i

In log-space: R_n = sum_{i=1}^{n} exp(logM_n - logM_{i-1})

This is always >= the CUSUM statistic (sum >= max), giving a slightly different power/false-alarm trade-off.

Parameters:

Name Type Description Default
martingale ConformalTestMartingale

The underlying martingale to wrap.

 required
References

Vovk, Gammerman & Shafer (2022). Algorithmic Learning in a Random World, 2nd edition, §8.3. Cambridge University Press.

Examples:

>>> from online_cp.martingale import SimpleJumper, ShiryaevRobertsWrapper
>>> sj = SimpleJumper(J=0.01)
>>> sr = ShiryaevRobertsWrapper(sj)
>>> for _ in range(10):
...     sr.update(0.5)
>>> sr.R >= 0
True
Source code in src/online_cp/martingale/wrappers.py 
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
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
class ShiryaevRobertsWrapper:
    """Shiryaev-Roberts change-detection wrapper for any conformal test martingale.

    Computes the Shiryaev-Roberts statistic as:

        R_n = sum_{i=1}^{n} S_n / S_i

    In log-space: R_n = sum_{i=1}^{n} exp(logM_n - logM_{i-1})

    This is always >= the CUSUM statistic (sum >= max), giving a slightly
    different power/false-alarm trade-off.

    Parameters
    ----------
    martingale : ConformalTestMartingale
        The underlying martingale to wrap.

    References
    ----------
    Vovk, Gammerman & Shafer (2022). *Algorithmic Learning in a Random World*,
    2nd edition, §8.3. Cambridge University Press.

    Examples
    --------
    >>> from online_cp.martingale import SimpleJumper, ShiryaevRobertsWrapper
    >>> sj = SimpleJumper(J=0.01)
    >>> sr = ShiryaevRobertsWrapper(sj)
    >>> for _ in range(10):
    ...     sr.update(0.5)
    >>> sr.R >= 0
    True
    """

    def __init__(self, martingale):
        """Wrap a martingale with the Shiryaev–Roberts statistic.

        Parameters
        ----------
        martingale : ConformalTestMartingale
            The underlying conformal test martingale whose increments drive the
            Shiryaev–Roberts statistic.
        """
        self.martingale = martingale
        self._n = 0
        self._sr_values = [0.0]  # R_0 = 0 (no terms in the sum)
        self._prev_logM = 0.0  # logM_{n-1}, starts at 0

    @property
    def R(self):
        """Current Shiryaev-Roberts statistic."""
        return self._sr_values[-1]

    @property
    def sr_values(self):
        """All Shiryaev-Roberts statistic values."""
        return self._sr_values

    def update(self, p: float) -> None:
        """Update the inner martingale and recompute the SR statistic.

        Uses the O(1) recursive formula (eq. 8.18 of ALRW2):
            R_n = (S_n / S_{n-1}) * (R_{n-1} + 1)
        """
        self.martingale.update(p)
        self._n += 1

        logM_n = self.martingale.logM
        # S_n / S_{n-1} = exp(logM_n - logM_{n-1})
        ratio = np.exp(logM_n - self._prev_logM)
        R_n = ratio * (self._sr_values[-1] + 1)
        self._sr_values.append(float(R_n))
        self._prev_logM = logM_n

    def alarm(self, threshold):
        """Check whether R_n exceeds the threshold.

        Parameters
        ----------
        threshold : float
            The alarm threshold.

        Returns
        -------
        bool
            True if R_n exceeds the threshold.
        """
        return self.R > threshold

R property

Current Shiryaev-Roberts statistic.

sr_values property

All Shiryaev-Roberts statistic values.

__init__(martingale)

Wrap a martingale with the Shiryaev–Roberts statistic.

Parameters:

Name Type Description Default
martingale ConformalTestMartingale

The underlying conformal test martingale whose increments drive the Shiryaev–Roberts statistic.

 required
Source code in src/online_cp/martingale/wrappers.py 
242
243
244
245
246
247
248
249
250
251
252
253
254
def __init__(self, martingale):
    """Wrap a martingale with the Shiryaev–Roberts statistic.

    Parameters
    ----------
    martingale : ConformalTestMartingale
        The underlying conformal test martingale whose increments drive the
        Shiryaev–Roberts statistic.
    """
    self.martingale = martingale
    self._n = 0
    self._sr_values = [0.0]  # R_0 = 0 (no terms in the sum)
    self._prev_logM = 0.0  # logM_{n-1}, starts at 0

update(p: float) -> None

Update the inner martingale and recompute the SR statistic.

Uses the O(1) recursive formula (eq. 8.18 of ALRW2): R_n = (S_n / S_{n-1}) * (R_{n-1} + 1)

Source code in src/online_cp/martingale/wrappers.py 
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
def update(self, p: float) -> None:
    """Update the inner martingale and recompute the SR statistic.

    Uses the O(1) recursive formula (eq. 8.18 of ALRW2):
        R_n = (S_n / S_{n-1}) * (R_{n-1} + 1)
    """
    self.martingale.update(p)
    self._n += 1

    logM_n = self.martingale.logM
    # S_n / S_{n-1} = exp(logM_n - logM_{n-1})
    ratio = np.exp(logM_n - self._prev_logM)
    R_n = ratio * (self._sr_values[-1] + 1)
    self._sr_values.append(float(R_n))
    self._prev_logM = logM_n

alarm(threshold)

Check whether R_n exceeds the threshold.

Parameters:

Name Type Description Default
threshold float

The alarm threshold.

 required

Returns:

Type Description
bool

True if R_n exceeds the threshold.
Source code in src/online_cp/martingale/wrappers.py 
282
283
284
285
286
287
288
289
290
291
292
293
294
295
def alarm(self, threshold):
    """Check whether R_n exceeds the threshold.

    Parameters
    ----------
    threshold : float
        The alarm threshold.

    Returns
    -------
    bool
        True if R_n exceeds the threshold.
    """
    return self.R > threshold

For betting strategies (density estimators used by martingales), see Betting Strategies.