eproc

Unnamed repository; edit this file 'description' to name the repository.
Log | Files | Refs | README | LICENSE

commit 3d8e7f8c1a33691e3d9b46e483a280f496a96817
parent ffcd36e3be7a982e7e551e1173a2a69530280292
Author: Jared Tobin <jared@jtobin.io>
Date:   Thu,  2 Jul 2026 16:18:01 -0230

Bernoulli: add two-sided variant via convex hedge

Adds 'Numeric.Eproc.Bernoulli.TwoSided', a sibling module to
'Numeric.Eproc.Bernoulli' testing H_0: p = p_0 against p != p_0.
The primary use case is the sign test at p_0 = 1/2, and more
generally any comparison of Bernoulli rates where the shift
direction isn't known in advance (which the existing one-sided
variant can't handle).

Construction is the convex hedge of Waudby-Smith & Ramdas
(2024) §4, mirroring what 'Bounded' does: two per-direction
Bernoulli capital processes K^+ (betting against p > p_0) and
K^- (betting against p < p_0) are combined into the hedged
e-process (K^+ + K^-) / 2 with E[K_0] = 1, and the test rejects
when the supremum of log(K^+ + K^-) has ever crossed
log(2/alpha). State carries per-direction log-wealths plus a
single max-log-sum field, updated each step via log-sum-exp
(moved from 'Bounded' into 'Common' so both consumers can use
it). Rejection is latched.

The sibling module keeps the same 'config' / 'initial' /
'update' / 'decide' / 'log_wealth' / 'samples' shape as the
one-sided version, so downstream code just qualifies both and
picks:

  import qualified Numeric.Eproc.Bernoulli as Bern
  import qualified Numeric.Eproc.Bernoulli.TwoSided as BernTS

Seven new test cases cover upward/downward shift detection,
null calibration, latched rejection, and config validation.

Diffstat:
MCHANGELOG | 8++++++++
Mlib/Numeric/Eproc/Bernoulli.hs | 10++++++----
Alib/Numeric/Eproc/Bernoulli/TwoSided.hs | 226+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Mlib/Numeric/Eproc/Bounded.hs | 10+---------
Mlib/Numeric/Eproc/Common.hs | 13+++++++++++++
Mppad-eproc.cabal | 5+++--
Mtest/Main.hs | 84+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
7 files changed, 341 insertions(+), 15 deletions(-)

diff --git a/CHANGELOG b/CHANGELOG @@ -1,5 +1,13 @@ # Changelog +- unreleased + * New module 'Numeric.Eproc.Bernoulli.TwoSided' providing a + two-sided Bernoulli rate test (H_0: p = p_0) via the same + convex-hedge construction 'Bounded' uses. Canonical use is the + sign test at p_0 = 1/2. Same 'config' / 'initial' / 'update' / + 'decide' / 'log_wealth' / 'samples' shape as the sibling + one-sided 'Numeric.Eproc.Bernoulli'. + - 0.2.1 (2026-07-02) * Two-sided bounded-mean tests now reject faster, or at least never later. diff --git a/lib/Numeric/Eproc/Bernoulli.hs b/lib/Numeric/Eproc/Bernoulli.hs @@ -8,7 +8,9 @@ -- License: MIT -- Maintainer: Jared Tobin <jared@ppad.tech> -- --- One-sided Bernoulli rate anytime-valid test. +-- One-sided Bernoulli rate anytime-valid test. See +-- "Numeric.Eproc.Bernoulli.TwoSided" for the two-sided companion +-- (used for the sign test at @p_0 = 1\/2@, among other things). -- -- For samples @x_t@ in @{0, 1}@, tests -- @@ -36,9 +38,9 @@ -- 'Reject' even if subsequent observations drive the current -- wealth back below threshold. -- --- Unlike "Numeric.Eproc.Bounded", the alternative here is one-sided, --- so a single wealth process suffices and no Bonferroni adjustment --- is needed -- the rejection threshold is @log(1 \/ alpha)@. +-- The alternative here is one-sided, so a single wealth process +-- suffices and no Bonferroni or hedge adjustment is needed -- the +-- rejection threshold is @log(1 \/ alpha)@. -- -- == Example -- diff --git a/lib/Numeric/Eproc/Bernoulli/TwoSided.hs b/lib/Numeric/Eproc/Bernoulli/TwoSided.hs @@ -0,0 +1,226 @@ +{-# OPTIONS_HADDOCK prune #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE RecordWildCards #-} + +-- | +-- Module: Numeric.Eproc.Bernoulli.TwoSided +-- Copyright: (c) 2026 Jared Tobin +-- License: MIT +-- Maintainer: Jared Tobin <jared@ppad.tech> +-- +-- Two-sided Bernoulli rate anytime-valid test. Companion to +-- "Numeric.Eproc.Bernoulli", which handles the one-sided case; +-- reach for this module when you want to test +-- +-- @H_0: E[x_t | F_{t-1}] = p_0 for all t@ +-- +-- against the negation. The canonical case is the sign test at +-- @p_0 = 1\/2@. +-- +-- The construction is the convex hedge of Waudby-Smith & Ramdas +-- (2024) §4: two per-direction Bernoulli capital processes +-- @K^+_t@ (betting against @p > p_0@ via @z = x - p_0@) and +-- @K^-_t@ (betting against @p < p_0@ via @-z@) are combined into +-- the hedged e-process @K_t = (K^+_t + K^-_t) \/ 2@ with +-- @E[K_0] = 1@. By Ville's inequality +-- @P(sup_t K_t >= 1 \/ alpha) <= alpha@, so the test rejects when +-- the supremum of @K^+_t + K^-_t@ has ever crossed @2 \/ alpha@; +-- the threshold is @log(2 \/ alpha)@. This is the same construction +-- "Numeric.Eproc.Bounded" uses to combine its two directional +-- processes. +-- +-- The test is /anytime-valid/ and rejection is /latched/ in the +-- running state. +-- +-- == Example +-- +-- Sign test at @p_0 = 1\/2@ with a downward shift: +-- +-- >>> let Right cfg = config 0.5 1.0e-3 Newton +-- >>> let s0 = initial cfg +-- >>> let xs = take 500 (cycle [False, False, False, True]) +-- >>> decide cfg (foldl' (update cfg) s0 xs) +-- Reject + +module Numeric.Eproc.Bernoulli.TwoSided ( + -- * Test configuration and state + Config + , State + , Verdict(..) + , ConfigError(..) + + -- * Bettor strategies + , Bettor(..) + + -- * Construction + , config + , initial + + -- * Streaming + , update + , decide + + -- * Inspection + , log_wealth + , samples + ) where + +import Numeric.Eproc.Common ( + Bettor(..), Verdict(..), ConfigError(..) + , BetState, init_bet, bet_lambda, step_bet + , finite, log_sum_exp + ) + +-- types ---------------------------------------------------------------------- + +-- | Two-sided Bernoulli rate test configuration. Build with 'config'. +-- +-- Carries the bettor strategy, the baseline rate, the significance +-- level, the precomputed convex-hedge log-wealth threshold +-- @log(2 \/ alpha)@, and the per-direction safe-bet ceilings. +data Config = Config { + cfg_bettor :: !Bettor + , cfg_lam_max_pos :: {-# UNPACK #-} !Double -- 0.5 / p0 + , cfg_lam_max_neg :: {-# UNPACK #-} !Double -- 0.5 / (1 - p0) + , cfg_p0 :: {-# UNPACK #-} !Double + , cfg_alpha :: {-# UNPACK #-} !Double + , cfg_log_thresh :: {-# UNPACK #-} !Double -- log(2/alpha) + } + +-- | Streaming test state. Construct with 'initial' and fold +-- observations through 'update'. +-- +-- The two log-wealth fields track the running log-wealth of the +-- positive- and negative-direction Bernoulli e-processes +-- separately; the /max log-sum/ field latches the supremum so +-- far of @log(K^+_t + K^-_t)@, which is the statistic the +-- convex-hedge construction actually monitors. +data State = State { + st_n :: {-# UNPACK #-} !Int + , st_log_w_pos :: {-# UNPACK #-} !Double + , st_log_w_neg :: {-# UNPACK #-} !Double + , st_max_log_sum :: {-# UNPACK #-} !Double + , st_bet_pos :: !BetState + , st_bet_neg :: !BetState + } + +-- construction --------------------------------------------------------------- + +-- | Build a 'Config' for the two-sided Bernoulli rate test. +-- +-- Per-direction safe-bet ceilings are @0.5 \/ p_0@ (positive) and +-- @0.5 \/ (1 - p_0)@ (negative), chosen so that each wealth factor +-- stays nonnegative for both admissible observations. The +-- threshold is @log(2 \/ alpha)@; the 2 reflects that the +-- convex-hedge test monitors the sum @K^+ + K^-@, whose initial +-- value is @2@ (each side starts at @K = 1@). +-- +-- Returns 'Left' with a 'ConfigError' on inputs that would leave +-- the mathematical regime: either of @p_0@ or @alpha@ non-finite +-- (NaN or infinite); @p_0@ outside @(0, 1)@; or @alpha@ outside +-- @(0, 1)@. +-- +-- >>> let Right cfg = config 0.5 1.0e-3 Newton +config + :: Double -- ^ baseline rate @p_0@, in @(0, 1)@ + -> Double -- ^ significance level @alpha@, in @(0, 1)@ + -> Bettor -- ^ bettor strategy + -> Either ConfigError Config +config !p0 !alpha !b + | not (finite p0 && p0 > 0 && p0 < 1) = + Left (InvalidBaselineRate p0) + | not (finite alpha && alpha > 0 && alpha < 1) = + Left (InvalidAlpha alpha) + | otherwise = Right Config { + cfg_bettor = b + , cfg_lam_max_pos = 0.5 / p0 + , cfg_lam_max_neg = 0.5 / (1 - p0) + , cfg_p0 = p0 + , cfg_alpha = alpha + , cfg_log_thresh = log (2 / alpha) + } +{-# INLINE config #-} + +-- | The initial 'State' for a fresh streaming test. +-- +-- Both per-direction log-wealths start at @0@ (i.e., @K = 1@); +-- the max-log-sum starts at @log 2@ (since @K^+_0 + K^-_0 = 2@); +-- both bettors start in the per-strategy initial state +-- appropriate for the 'Bettor' chosen in the 'Config'. +-- +-- >>> let s0 = initial cfg +initial :: Config -> State +initial Config{..} = + let !s0 = init_bet cfg_bettor + in State { + st_n = 0 + , st_log_w_pos = 0 + , st_log_w_neg = 0 + , st_max_log_sum = log 2 + , st_bet_pos = s0 + , st_bet_neg = s0 + } +{-# INLINE initial #-} + +-- streaming ------------------------------------------------------------------ + +-- | Fold one observation into the running 'State'. +-- +-- Computes the centred observation @z = x - p_0@, queries the two +-- directional bettors, accumulates per-direction log-wealth, then +-- updates the running supremum of @log(K^+ + K^-)@ via +-- log-sum-exp and steps the bettor states. +-- +-- >>> let s1 = update cfg s0 True +update :: Config -> State -> Bool -> State +update Config{..} State{..} !x = + let !xd = if x then 1 else 0 + !z = xd - cfg_p0 + !lam_p = bet_lambda cfg_bettor cfg_lam_max_pos st_bet_pos + !lam_n = bet_lambda cfg_bettor cfg_lam_max_neg st_bet_neg + !fac_p = 1 + lam_p * z + !fac_n = 1 - lam_n * z + !logw_p = st_log_w_pos + log fac_p + !logw_n = st_log_w_neg + log fac_n + !log_sum = log_sum_exp logw_p logw_n + !max_sum = max st_max_log_sum log_sum + !sp = step_bet cfg_bettor cfg_lam_max_pos st_bet_pos z + !sn = step_bet cfg_bettor cfg_lam_max_neg st_bet_neg (negate z) + in State (st_n + 1) logw_p logw_n max_sum sp sn +{-# INLINE update #-} + +-- | Compute the current 'Verdict' from the running 'State'. +-- +-- 'Reject' iff the supremum-so-far of @log(K^+_t + K^-_t)@ has +-- crossed @log(2 \/ alpha)@ at some point; equivalently the +-- convex-hedge e-process @(K^+ + K^-) \/ 2@ has exceeded +-- @1 \/ alpha@. Under @H_0@, Ville's inequality bounds the +-- probability of this ever happening by @alpha@, uniformly +-- across sample counts. +-- +-- >>> decide cfg s0 +-- Continue +decide :: Config -> State -> Verdict +decide Config{..} State{..} + | st_max_log_sum >= cfg_log_thresh = Reject + | otherwise = Continue +{-# INLINE decide #-} + +-- inspection ----------------------------------------------------------------- + +-- | The supremum-so-far of @log(K^+_t + K^-_t)@. Monotone +-- nondecreasing; starts at @log 2@ (since @K^+_0 + K^-_0 = 2@). +-- +-- >>> log_wealth s0 +-- 0.6931471805599453 +log_wealth :: State -> Double +log_wealth = st_max_log_sum +{-# INLINE log_wealth #-} + +-- | The number of samples consumed so far. +-- +-- >>> samples s0 +-- 0 +samples :: State -> Int +samples = st_n +{-# INLINE samples #-} diff --git a/lib/Numeric/Eproc/Bounded.hs b/lib/Numeric/Eproc/Bounded.hs @@ -86,11 +86,10 @@ module Numeric.Eproc.Bounded ( , samples ) where -import GHC.Float (log1p) import Numeric.Eproc.Common ( Bettor(..), Verdict(..), ConfigError(..) , BetState, init_bet, bet_lambda, step_bet - , finite + , finite, log_sum_exp ) -- types ---------------------------------------------------------------------- @@ -254,13 +253,6 @@ update Config{..} State{..} !x = in State (st_n + 1) logw_p logw_n max_sum sp sn {-# INLINE update #-} --- | @log(exp a + exp b)@, computed without intermediate overflow. -log_sum_exp :: Double -> Double -> Double -log_sum_exp !a !b - | a >= b = a + log1p (exp (b - a)) - | otherwise = b + log1p (exp (a - b)) -{-# INLINE log_sum_exp #-} - -- | Compute the current 'Verdict' from the running 'State'. -- -- 'Reject' iff the supremum-so-far of @log(K^+_t + K^-_t)@ has diff --git a/lib/Numeric/Eproc/Common.hs b/lib/Numeric/Eproc/Common.hs @@ -32,8 +32,11 @@ module Numeric.Eproc.Common ( -- * Internal: helpers , finite + , log_sum_exp ) where +import GHC.Float (log1p) + -- | A predictable bettor. -- -- A bettor describes how, given the history of centred @@ -117,6 +120,16 @@ finite :: Double -> Bool finite x = not (isNaN x) && not (isInfinite x) {-# INLINE finite #-} +-- | @log(exp a + exp b)@, computed without intermediate overflow. +-- Used by the convex-hedge two-sided combinations to update the +-- running @log(K^+ + K^-)@ statistic from the two per-direction +-- log-wealths. +log_sum_exp :: Double -> Double -> Double +log_sum_exp !a !b + | a >= b = a + log1p (exp (b - a)) + | otherwise = b + log1p (exp (a - b)) +{-# INLINE log_sum_exp #-} + -- | Per-bettor state. One constructor per 'Bettor' alternative; the -- constructor used in any given state matches the 'Bettor' chosen -- in the enclosing 'Config'. diff --git a/ppad-eproc.cabal b/ppad-eproc.cabal @@ -13,8 +13,8 @@ extra-doc-files: CHANGELOG description: Anytime-valid sequential hypothesis testing for bounded random variables, via the e-process / betting framework of Waudby-Smith and - Ramdas (2024). Provides bounded-mean, paired two-sample, and - one-sided Bernoulli rate tests with fixed, adaptive (aGRAPA), and + Ramdas (2024). Provides bounded-mean, paired two-sample, and one- and + two-sided Bernoulli rate tests with fixed, adaptive (aGRAPA), and online Newton bettors. flag llvm @@ -35,6 +35,7 @@ library ghc-options: -fllvm -O2 exposed-modules: Numeric.Eproc.Bernoulli + Numeric.Eproc.Bernoulli.TwoSided Numeric.Eproc.Bounded Numeric.Eproc.Common Numeric.Eproc.Paired diff --git a/test/Main.hs b/test/Main.hs @@ -5,6 +5,7 @@ module Main where import Data.Bits import Data.Word import qualified Numeric.Eproc.Bernoulli as Bern +import qualified Numeric.Eproc.Bernoulli.TwoSided as BernTS import qualified Numeric.Eproc.Bounded as Bounded import qualified Numeric.Eproc.Common as C import qualified Numeric.Eproc.Paired as P @@ -23,6 +24,7 @@ main = defaultMain $ testGroup "ppad-eproc" [ , latched_rejection_tests , config_validation_tests , safety_property_tests + , two_sided_bernoulli_tests ] -- partial helper: tests below hardcode valid configs. @@ -391,6 +393,88 @@ config_validation_tests = testGroup "config validation" [ Left _ -> pure () Right _ -> assertFailure "expected Left" +-- two-sided bernoulli -------------------------------------------------------- + +run_ts_bernoulli + :: BernTS.Config + -> Double -- ^ true rate p + -> Int -- ^ budget + -> Gen + -> (BernTS.Verdict, Int) +run_ts_bernoulli cfg p budget g0 = + go 0 g0 (BernTS.initial cfg) + where + go !n !g !st + | n >= budget = (BernTS.decide cfg st, n) + | otherwise = case BernTS.decide cfg st of + BernTS.Reject -> (BernTS.Reject, n) + BernTS.Continue -> + let (u, g') = next_double g + !x = u < p + st' = BernTS.update cfg st x + in go (n + 1) g' st' + +ts_bernoulli_rate + :: BernTS.Config + -> Double + -> Int + -> Int + -> Word64 + -> Double +ts_bernoulli_rate cfg p budget trials seed = + let gens = take trials (gen_seq (mk_gen seed)) + rejects = length + [ () | g <- gens + , let (v, _) = run_ts_bernoulli cfg p budget g + , v == BernTS.Reject ] + in fromIntegral rejects / fromIntegral trials + +two_sided_bernoulli_tests :: TestTree +two_sided_bernoulli_tests = testGroup "two-sided bernoulli" [ + testCase "constant at p_0 doesn't reject" $ do + -- Bernoulli(0.5) with p_0 = 0.5 is under the null. + let cfg = ok (BernTS.config 0.5 1.0e-6 BernTS.Newton) + -- alternating True/False keeps the empirical rate at 0.5. + xs = take 5000 (cycle [True, False]) + st = foldl' (BernTS.update cfg) (BernTS.initial cfg) xs + BernTS.decide cfg st @?= BernTS.Continue + , testCase "detects upward shift (p = 0.7 vs p_0 = 0.5)" $ do + let cfg = ok (BernTS.config 0.5 1.0e-3 BernTS.Newton) + rate = ts_bernoulli_rate cfg 0.7 5000 100 111222 + assertBool ("power " ++ show rate ++ " too low") $ + rate >= 0.95 + , testCase "detects downward shift (p = 0.3 vs p_0 = 0.5)" $ do + let cfg = ok (BernTS.config 0.5 1.0e-3 BernTS.Newton) + rate = ts_bernoulli_rate cfg 0.3 5000 100 333444 + assertBool ("power " ++ show rate ++ " too low") $ + rate >= 0.95 + , testCase "FPR at p = p_0 = 0.5 within slack" $ do + let cfg = ok (BernTS.config 0.5 0.05 BernTS.Newton) + rate = ts_bernoulli_rate cfg 0.5 2000 200 555666 + assertBool ("FPR " ++ show rate ++ " exceeded slack") $ + rate <= 0.08 + , testCase "latched: cross then drown stays rejected" $ do + let cfg = ok (BernTS.config 0.5 0.5 (BernTS.Fixed 1.0)) + -- ten 1s push the positive side well past threshold. + xs1 = replicate 10 True + -- then two hundred 0s drop the current wealth, but the + -- latch must hold. + xs2 = replicate 200 False + st1 = foldl' (BernTS.update cfg) (BernTS.initial cfg) xs1 + st2 = foldl' (BernTS.update cfg) st1 xs2 + BernTS.decide cfg st1 @?= BernTS.Reject + BernTS.decide cfg st2 @?= BernTS.Reject + , testCase "config: NaN p0 rejected" $ do + let nan = 0/0 :: Double + case BernTS.config nan 0.05 BernTS.Newton of + Left _ -> pure () + Right _ -> assertFailure "expected Left" + , testCase "config: alpha out of range rejected" $ + case BernTS.config 0.5 1.5 BernTS.Newton of + Left _ -> pure () + Right _ -> assertFailure "expected Left" + ] + -- safety properties ---------------------------------------------------------- unit_double :: QC.Gen Double