ppad-eproc 0.3.0 → 0.4.0
raw patch · 12 files changed
+1291/−12 lines, 12 files
Files
- CHANGELOG +11/−4
- bench/Main.hs +49/−0
- bench/Weight.hs +44/−0
- lib/Numeric/Eproc/Bernoulli.hs +44/−0
- lib/Numeric/Eproc/Bernoulli/TwoSided.hs +35/−0
- lib/Numeric/Eproc/Bounded.hs +44/−0
- lib/Numeric/Eproc/Common.hs +23/−2
- lib/Numeric/Eproc/ConfSeq.hs +327/−0
- lib/Numeric/Eproc/Mixture.hs +297/−0
- lib/Numeric/Eproc/Paired.hs +36/−0
- ppad-eproc.cabal +10/−6
- test/Main.hs +371/−0
CHANGELOG view
@@ -1,13 +1,20 @@ # Changelog +- 0.4.0 (2026-07-03)+ * Adds calibrated evidence accessors to every test module:+ 'log_evalue', 'log_evalue_sup', and the anytime-valid 'p_value'.+ * Adds Numeric.Eproc.Mixture: uniform convex mixtures of+ e-processes, for testing a null against a union of qualitatively+ different alternatives at a single Ville threshold.+ * Adds Numeric.Eproc.ConfSeq: anytime-valid confidence sequences+ for bounded means, via the hedged-capital construction of+ Waudby-Smith & Ramdas (2024).+ * Adds InvalidComponentCount and InvalidGridSize to ConfigError.+ - 0.3.0 (2026-07-02) * Introduces a breaking API change: 'log_wealth' now returns the current log-wealth, whereas the supremum-thus-far statistic is exposed as 'log_wealth_sup'.--- 0.2.2 (2026-07-02)- * Adds a Numeric.Eproc.Bernoulli.TwoSided module for a two-sided- Bernoulli rate test. - 0.2.2 (2026-07-02) * Adds a Numeric.Eproc.Bernoulli.TwoSided module for a two-sided
bench/Main.hs view
@@ -7,6 +7,8 @@ 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.ConfSeq as CS+import qualified Numeric.Eproc.Mixture as Mix import qualified Numeric.Eproc.Paired as P import Criterion.Main @@ -17,6 +19,7 @@ instance NFData P.State where rnf !_ = () instance NFData Bern.State where rnf !_ = () instance NFData BernTS.State where rnf !_ = ()+instance NFData Mix.State where rnf !_ = () instance NFData Bounded.Verdict where rnf !_ = () -- partial helper for benches: configs here are hardcoded valid, so a@@ -35,6 +38,10 @@ , bern_stream , bern_ts_update , bern_ts_stream+ , mix_update+ , mix_stream+ , confseq_update+ , confseq_stream ] update :: Benchmark@@ -138,4 +145,46 @@ bench "fixed" $ nf (run_b cfg_f) xs , bench "adaptive" $ nf (run_b cfg_a) xs , bench "newton" $ nf (run_b cfg_o) xs+ ]++mix_update :: Benchmark+mix_update =+ let !cfg = ok (Mix.config 4 1.0e-3)+ !st = Mix.initial cfg+ !v = force [0.1, -0.2, 0.3, 0.0]+ in bgroup "Mixture.update (one step)" [+ bench "K=4" $ nf (Mix.update cfg st) v+ ]++mix_stream :: Benchmark+mix_stream =+ let !vs = force (take 1000 (cycle+ [[0.1, -0.2, 0.3, 0.0], [-0.3, 0.2, 0.0, 0.1]]))+ !cfg = ok (Mix.config 4 1.0e-3)+ run_x c = foldl' (Mix.update c) (Mix.initial c)+ in bgroup "Mixture.update (1000-step fold)" [+ bench "K=4" $ nf (run_x cfg) vs+ ]++-- ConfSeq.State carries a list of live grid candidates rather than+-- only unboxed fields, but 'initial' and 'update' construct that+-- list fully forced, so WHNF == NF holds here by construction too.+instance NFData CS.State where rnf !_ = ()++confseq_update :: Benchmark+confseq_update =+ let !cfg = ok (CS.config 0.0 1.0 0.05 200)+ !st = CS.initial cfg+ !x = 0.7+ in bgroup "ConfSeq.update (one step, g = 200)" [+ bench "plug-in" $ nf (CS.update cfg st) x+ ]++confseq_stream :: Benchmark+confseq_stream =+ let !xs = force (take 1000 (cycle [0.3, 0.7]))+ !cfg = ok (CS.config 0.0 1.0 0.05 200)+ run_c = foldl' (CS.update cfg) (CS.initial cfg)+ in bgroup "ConfSeq.update (1000-sample fold, g = 200)" [+ bench "plug-in" $ nf run_c xs ]
bench/Weight.hs view
@@ -7,6 +7,8 @@ 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.ConfSeq as CS+import qualified Numeric.Eproc.Mixture as Mix import qualified Numeric.Eproc.Paired as P import Weigh @@ -14,6 +16,7 @@ instance NFData P.State where rnf !_ = () instance NFData Bern.State where rnf !_ = () instance NFData BernTS.State where rnf !_ = ()+instance NFData Mix.State where rnf !_ = () instance NFData Bounded.Verdict where rnf !_ = () -- partial helper for benches: configs here are hardcoded valid.@@ -32,6 +35,10 @@ bern_stream bern_ts_update bern_ts_stream+ mix_update+ mix_stream+ confseq_update+ confseq_stream update :: Weigh () update =@@ -126,3 +133,40 @@ func "fixed" (run_b cfg_f) xs func "adaptive" (run_b cfg_a) xs func "newton" (run_b cfg_o) xs++mix_update :: Weigh ()+mix_update =+ let !cfg = ok (Mix.config 4 1.0e-3)+ !st = Mix.initial cfg+ !v = force [0.1, -0.2, 0.3, 0.0]+ in wgroup "Mixture.update (one step)" $ do+ func "K=4" (Mix.update cfg st) v++mix_stream :: Weigh ()+mix_stream =+ let !vs = force (take 1000 (cycle+ [[0.1, -0.2, 0.3, 0.0], [-0.3, 0.2, 0.0, 0.1]]))+ !cfg = ok (Mix.config 4 1.0e-3)+ run_x c = foldl' (Mix.update c) (Mix.initial c)+ in wgroup "Mixture.update (1000-step fold)" $ do+ func "K=4" (run_x cfg) vs++-- ConfSeq.State carries a list of live grid candidates rather than+-- only unboxed fields, but 'initial' and 'update' construct that+-- list fully forced, so WHNF == NF holds here by construction too.+instance NFData CS.State where rnf !_ = ()++confseq_update :: Weigh ()+confseq_update =+ let !cfg = ok (CS.config 0.0 1.0 0.05 200)+ !st = CS.initial cfg+ in wgroup "ConfSeq.update (one step, g = 200)" $ do+ func "plug-in" (CS.update cfg st) 0.7++confseq_stream :: Weigh ()+confseq_stream =+ let !xs = force (take 1000 (cycle [0.3, 0.7]))+ !cfg = ok (CS.config 0.0 1.0 0.05 200)+ run_c = foldl' (CS.update cfg) (CS.initial cfg)+ in wgroup "ConfSeq.update (1000-sample fold, g = 200)" $ do+ func "plug-in" run_c xs
lib/Numeric/Eproc/Bernoulli.hs view
@@ -73,6 +73,9 @@ -- * Inspection , log_wealth , log_wealth_sup+ , log_evalue+ , log_evalue_sup+ , p_value , samples ) where @@ -254,6 +257,47 @@ log_wealth_sup :: State -> Double log_wealth_sup = st_sup_log_w {-# INLINE log_wealth_sup #-}++-- | The current log e-value. For this one-sided test the single+-- wealth process is itself the e-process (a fresh state already+-- sits at wealth @1@), so this coincides with 'log_wealth'; the+-- accessor exists so that e-values read uniformly across test+-- modules regardless of their internal hedging, e.g. when+-- convex-combining several e-processes. Not monotone; bounded+-- above by 'log_evalue_sup'.+--+-- >>> log_evalue s0+-- 0.0+log_evalue :: State -> Double+log_evalue = st_log_w+{-# INLINE log_evalue #-}++-- | The supremum-so-far of the log e-value; coincides with+-- 'log_wealth_sup' for this one-sided test. Monotone+-- nondecreasing, starting at @0@; 'decide' rejects exactly when+-- it crosses @log(1 \/ alpha)@.+--+-- >>> log_evalue_sup s0+-- 0.0+log_evalue_sup :: State -> Double+log_evalue_sup = st_sup_log_w+{-# INLINE log_evalue_sup #-}++-- | The anytime-valid p-value: the reciprocal of the largest+-- e-value attained so far, @min 1 (exp (negate (log_evalue_sup+-- s)))@.+--+-- Monotone nonincreasing in the sample count, and valid under+-- optional stopping: under @H_0@,+-- @P(exists t: p_t <= alpha) <= alpha@ for every @alpha@+-- simultaneously. 'decide' returns 'Reject' exactly when this+-- value has reached the configured @alpha@ or below.+--+-- >>> p_value s0+-- 1.0+p_value :: State -> Double+p_value s = min 1 (exp (negate (log_evalue_sup s)))+{-# INLINE p_value #-} -- | The number of samples consumed so far. --
lib/Numeric/Eproc/Bernoulli/TwoSided.hs view
@@ -55,6 +55,9 @@ -- * Inspection , log_wealth , log_wealth_sup+ , log_evalue+ , log_evalue_sup+ , p_value , samples ) where @@ -141,6 +144,38 @@ log_wealth_sup :: State -> Double log_wealth_sup (State s) = Bounded.log_wealth_sup s {-# INLINE log_wealth_sup #-}++-- | The current log e-value of the underlying bounded-mean test:+-- 'log_wealth' minus @log 2@, normalized so a fresh state sits at+-- @0@. Not monotone; bounded above by 'log_evalue_sup'.+--+-- >>> log_evalue s0+-- 0.0+log_evalue :: State -> Double+log_evalue (State s) = Bounded.log_evalue s+{-# INLINE log_evalue #-}++-- | The supremum-so-far of the log e-value: 'log_wealth_sup' minus+-- @log 2@. Monotone nondecreasing, starting at @0@; 'decide'+-- rejects exactly when it crosses @log(1 \/ alpha)@.+--+-- >>> log_evalue_sup s0+-- 0.0+log_evalue_sup :: State -> Double+log_evalue_sup (State s) = Bounded.log_evalue_sup s+{-# INLINE log_evalue_sup #-}++-- | The anytime-valid p-value: the reciprocal of the largest+-- e-value attained so far. Monotone nonincreasing; under @H_0@,+-- @P(exists t: p_t <= alpha) <= alpha@ for every @alpha@+-- simultaneously. 'decide' returns 'Reject' exactly when this+-- value has reached the configured @alpha@ or below.+--+-- >>> p_value s0+-- 1.0+p_value :: State -> Double+p_value (State s) = Bounded.p_value s+{-# INLINE p_value #-} -- | The number of samples consumed so far. --
lib/Numeric/Eproc/Bounded.hs view
@@ -84,6 +84,9 @@ -- * Inspection , log_wealth , log_wealth_sup+ , log_evalue+ , log_evalue_sup+ , p_value , samples ) where @@ -310,6 +313,47 @@ log_wealth_sup :: State -> Double log_wealth_sup State{..} = st_sup_log_sum {-# INLINE log_wealth_sup #-}++-- | The current log e-value of the convex-hedge e-process: the log+-- of @(K^+_t + K^-_t) \/ 2@, i.e. 'log_wealth' minus @log 2@.+--+-- Unlike 'log_wealth', this is normalized so that a fresh state+-- sits at @0@ (e-value @1@): it is directly comparable across+-- test modules regardless of their internal hedging, and is the+-- form to use when convex-combining several e-processes. Not+-- monotone; bounded above by 'log_evalue_sup'.+--+-- >>> log_evalue s0+-- 0.0+log_evalue :: State -> Double+log_evalue s = log_wealth s - log2_dbl+{-# INLINE log_evalue #-}++-- | The supremum-so-far of the log e-value: 'log_wealth_sup' minus+-- @log 2@. Monotone nondecreasing, starting at @0@; 'decide'+-- rejects exactly when it crosses @log(1 \/ alpha)@.+--+-- >>> log_evalue_sup s0+-- 0.0+log_evalue_sup :: State -> Double+log_evalue_sup s = log_wealth_sup s - log2_dbl+{-# INLINE log_evalue_sup #-}++-- | The anytime-valid p-value: the reciprocal of the largest+-- e-value attained so far, @min 1 (exp (negate (log_evalue_sup+-- s)))@.+--+-- Monotone nonincreasing in the sample count, and valid under+-- optional stopping: under @H_0@,+-- @P(exists t: p_t <= alpha) <= alpha@ for every @alpha@+-- simultaneously. 'decide' returns 'Reject' exactly when this+-- value has reached the configured @alpha@ or below.+--+-- >>> p_value s0+-- 1.0+p_value :: State -> Double+p_value s = min 1 (exp (negate (log_evalue_sup s)))+{-# INLINE p_value #-} -- | The number of samples consumed so far. --
lib/Numeric/Eproc/Common.hs view
@@ -73,6 +73,21 @@ -- rate @2 \/ (2 - log 3)@. Achieves logarithmic regret against -- the best constant bet in hindsight and is in practice the -- strongest of the three bettors under most signal regimes.+--+-- One deliberate deviation from WSR: Algorithm 2 seeds the+-- squared-gradient accumulator at @1@, which presumes+-- observations scaled to @[0, 1]@. On raw-scale data that+-- constant is dimensionally wrong -- negligible when+-- @z^2 >> 1@, paralysing when @z^2 << 1@ -- so the accumulator+-- here is instead seeded near zero, making the update+-- scale-adaptive. The trade is bold early play: the first+-- nonzero observation typically drives the bet straight to+-- the @lambda_max@ ceiling, annealing back toward the Kelly+-- point as gradients accumulate. Validity is unaffected --+-- predictability and clipping are all it needs -- and regret+-- stays logarithmic with a somewhat larger constant. The+-- visible effect is higher-variance early wealth: a supremum+-- modestly above its floor is expected even under @H_0@. data Bettor = Fixed {-# UNPACK #-} !Double | Adaptive@@ -97,8 +112,10 @@ -- | Reasons that a test-configuration smart constructor can reject -- its inputs. Returned by 'Numeric.Eproc.Bounded.config',--- 'Numeric.Eproc.Bernoulli.config', and--- 'Numeric.Eproc.Paired.config'.+-- 'Numeric.Eproc.Bernoulli.config',+-- 'Numeric.Eproc.Paired.config',+-- 'Numeric.Eproc.Mixture.config', and+-- 'Numeric.Eproc.ConfSeq.config'. data ConfigError = -- | significance level outside @(0, 1)@ InvalidAlpha {-# UNPACK #-} !Double@@ -112,6 +129,10 @@ {-# UNPACK #-} !Double -- hi -- | baseline rate outside @(0, 1)@ | InvalidBaselineRate {-# UNPACK #-} !Double+ -- | component count not positive+ | InvalidComponentCount {-# UNPACK #-} !Int+ -- | grid size below @1@+ | InvalidGridSize {-# UNPACK #-} !Int deriving (Eq, Show) -- | True iff the argument is a finite IEEE-754 double (not NaN, not
+ lib/Numeric/Eproc/ConfSeq.hs view
@@ -0,0 +1,327 @@+{-# OPTIONS_HADDOCK prune #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE RecordWildCards #-}++-- |+-- Module: Numeric.Eproc.ConfSeq+-- Copyright: (c) 2026 Jared Tobin+-- License: MIT+-- Maintainer: Jared Tobin <jared@ppad.tech>+--+-- Anytime-valid confidence sequence for the mean of bounded+-- observations.+--+-- For samples @x_t@ in @[lo, hi]@ with common conditional mean+--+-- @mu = E[x_t | F_{t-1}] for all t@+--+-- (@F_{t-1}@ being the filtration generated by everything observed+-- strictly before time @t@; for i.i.d. samples this is just+-- @E[x]@), the running state yields a confidence interval @C_t@+-- after every observation, with time-uniform coverage:+--+-- @P(for all t, mu in C_t) >= 1 - alpha@+--+-- whenever @C_t@ is reported at all (see 'interval' for the empty+-- case). The guarantee holds uniformly over time, so the user may+-- inspect the interval after every observation and stop at any+-- data-dependent time -- optional stopping does not erode coverage.+--+-- The construction is the /hedged capital/ confidence sequence of+-- Waudby-Smith & Ramdas (2024), Theorem 3, evaluated over a finite+-- grid of candidate means. All arithmetic is carried out in+-- @[0, 1]@ coordinates internally; observations are mapped affinely+-- at the boundary. Each candidate @m@ runs a pair of betting+-- processes: a /positive-direction/ capital @K^+_t(m)@ wagering+-- that the mean exceeds @m@, and a /negative-direction/ capital+-- @K^-_t(m)@ wagering the reverse. The base bet is a single+-- predictable plug-in (their eq. (26)), computed once per update+-- from the running regularized mean and variance of the data and+-- shared by every candidate: it never depends on @m@, and only a+-- final truncation to @c \/ m@ (respectively @c \/ (1 - m)@), with+-- @c = 1\/2@, is candidate-specific. This @m@-freeness is what+-- makes the survivor set provably an interval (Theorem 3);+-- @m@-dependent bets can produce non-interval survivor sets (their+-- Section E.4), which is why this module does not use the library's+-- 'Numeric.Eproc.Common.Bettor' strategies.+--+-- A candidate @m@ is rejected once the max-hedge (@theta = 1\/2@)+-- capital @max(K^+_t(m), K^-_t(m)) \/ 2@ crosses @1 \/ alpha@.+-- Under the truth @m = mu@ each capital process is a nonnegative+-- supermartingale, the max is dominated by the convex combination+-- @(K^+ + K^-) \/ 2@, and Ville's inequality bounds the probability+-- that the truth is ever rejected by @alpha@. No multiplicity+-- correction across grid candidates is needed: coverage concerns+-- only the true mean's own test, and rejection of other candidates+-- merely tightens the interval.+--+-- Grid resolution is an accuracy\/cost knob. Interval endpoints are+-- quantized to the grid -- a @g@-point grid resolves them to within+-- @(hi - lo) \/ (g + 1)@ -- and per-update cost is @O(live+-- candidates)@, shrinking as evidence accumulates and candidates+-- are rejected.+--+-- == Example+--+-- Estimate the mean of a stream in @[0, 1]@ with empirical mean+-- @0.8@, at level @alpha = 0.05@ on a 100-point grid:+--+-- >>> let Right cfg = config 0.0 1.0 0.05 100+-- >>> let xs = concat (replicate 50 [1, 1, 0, 1, 1, 0, 1, 1, 1, 1])+-- >>> interval cfg (foldl' (update cfg) (initial cfg) xs)+-- Just (0.7326732673267327,0.8514851485148515)++module Numeric.Eproc.ConfSeq (+ -- * Confidence-sequence configuration and state+ Config+ , State+ , ConfigError(..)++ -- * Construction+ , config+ , initial++ -- * Streaming+ , update++ -- * Inspection+ , interval+ , samples+ ) where++import GHC.Float (log1p)+import Numeric.Eproc.Common (ConfigError(..), finite)++-- types ----------------------------------------------------------------------++-- | Confidence-sequence configuration. Build with 'config'.+--+-- Carries the sample bounds, the significance level, the grid+-- size, and the precomputed per-candidate rejection threshold+-- @log(2 \/ alpha)@ along with the bet numerator+-- @2 log(2 \/ alpha)@.+data Config = Config {+ cfg_lo :: {-# UNPACK #-} !Double -- ^ sample lower bound+ , cfg_hi :: {-# UNPACK #-} !Double -- ^ sample upper bound+ , cfg_alpha :: {-# UNPACK #-} !Double -- ^ significance level+ , cfg_grid :: {-# UNPACK #-} !Int -- ^ grid size @g@+ , cfg_log_thresh :: {-# UNPACK #-} !Double -- ^ @log(2 \/ alpha)@+ , cfg_bet_num :: {-# UNPACK #-} !Double -- ^ @2 log(2 \/ alpha)@+ }++-- | One live grid candidate: its grid index and the running+-- log-capitals of the positive- and negative-direction bets.+data Point = Point+ {-# UNPACK #-} !Int -- grid index j+ {-# UNPACK #-} !Double -- log K^++ {-# UNPACK #-} !Double -- log K^-++-- | Streaming confidence-sequence state. Construct with 'initial'+-- and fold observations through 'update'.+--+-- Carries the sample count, the shared plug-in bettor statistics+-- (regularized running sums in @[0, 1]@ coordinates), and the+-- live grid candidates. Rejected candidates are dropped+-- permanently, so the reported intervals are nested.+--+-- Invariant: 'initial' and 'update' construct the live list fully+-- forced -- no thunks in the spine or the elements -- so a 'State'+-- in WHNF is already in normal form.+data State = State {+ st_n :: {-# UNPACK #-} !Int -- ^ sample count+ , st_sum_y :: {-# UNPACK #-} !Double -- ^ @sum y_i@+ , st_sum_dev2 :: {-# UNPACK #-} !Double -- ^ @sum (y_i - mu_i)^2@+ , st_live :: ![Point] -- ^ live grid candidates+ }++-- | WSR (2024) truncation level @c = 1\/2@. Bets are capped at+-- @c \/ m@ (positive direction) and @c \/ (1 - m)@ (negative+-- direction), keeping every capital factor at least @1 - c > 0@.+trunc_c :: Double+trunc_c = 0.5+{-# INLINE trunc_c #-}++-- construction ---------------------------------------------------------------++-- | Build a 'Config' for the confidence sequence.+--+-- The candidate means form the interior grid+--+-- @m_j = lo + (j \/ (g + 1)) * (hi - lo), j = 1 .. g@+--+-- (endpoints excluded, so that in @[0, 1]@ coordinates the bet+-- truncations @c \/ m@ and @c \/ (1 - m)@ stay finite). The+-- per-candidate rejection threshold @log(2 \/ alpha)@ and the bet+-- numerator @2 log(2 \/ alpha)@ are precomputed.+--+-- Returns 'Left' with a 'ConfigError' on inputs that would leave+-- the mathematical regime: @alpha@ non-finite or outside+-- @(0, 1)@; @lo@ or @hi@ non-finite, or @lo >= hi@; or a grid+-- size below @1@.+--+-- >>> let Right cfg = config 0.0 1.0 0.05 100+config+ :: Double -- ^ sample lower bound @lo@+ -> Double -- ^ sample upper bound @hi@+ -> Double -- ^ significance level @alpha@+ -> Int -- ^ grid size @g@+ -> Either ConfigError Config+config !lo !hi !alpha !g+ | not (finite alpha && alpha > 0 && alpha < 1) =+ Left (InvalidAlpha alpha)+ | not (finite lo && finite hi && lo < hi) =+ Left (InvalidBounds lo hi)+ | g < 1 =+ Left (InvalidGridSize g)+ | otherwise = Right Config {+ cfg_lo = lo+ , cfg_hi = hi+ , cfg_alpha = alpha+ , cfg_grid = g+ , cfg_log_thresh = log (2 / alpha)+ , cfg_bet_num = 2 * log (2 / alpha)+ }+{-# INLINE config #-}++-- | The initial 'State' for a fresh confidence sequence.+--+-- Every grid candidate starts live with both log-capitals at @0@+-- (i.e., @K^+ = K^- = 1@); the shared bettor statistics start+-- from their regularized priors (@mu_0 = 1\/2@,+-- @sigma^2_0 = 1\/4@ in @[0, 1]@ coordinates).+--+-- >>> let s0 = initial cfg+initial :: Config -> State+initial Config{..} = State {+ st_n = 0+ , st_sum_y = 0+ , st_sum_dev2 = 0+ , st_live = points 1+ }+ where+ -- built eagerly: the tail is forced before consing, so the+ -- whole list is in normal form on construction.+ points !j+ | j > cfg_grid = []+ | otherwise =+ let !p = Point j 0 0+ !rest = points (j + 1)+ in p : rest+{-# INLINE initial #-}++-- streaming ------------------------------------------------------------------++-- | Fold one observation into the running 'State'.+--+-- Maps the observation to @[0, 1]@ coordinates via+-- @y = (x - lo) \/ (hi - lo)@ and computes the shared predictable+-- plug-in bet from the statistics accumulated through the+-- /previous/ step (Waudby-Smith & Ramdas (2024), eq. (26)):+--+-- @lambda_t = min c (sqrt (2 log(2 \/ alpha)+-- \/ (sigma^2_{t-1} * t * log(1 + t))))@+--+-- with @c = 1\/2@. The bet is computed once and shared across all+-- live candidates -- its independence from @m@ is what keeps the+-- survivor set an interval. Each live candidate @m@ then updates+-- its pair of log-capitals with the truncated bets+-- @min lambda_t (c \/ m)@ and @min lambda_t (c \/ (1 - m))@, and+-- is dropped iff @max(log K^+, log K^-)@ has reached+-- @log(2 \/ alpha)@. Finally @y@ is folded into the shared+-- statistics, preserving predictability of the next bet.+--+-- /Precondition/: @x@ must lie in the @[lo, hi]@ interval given+-- to 'config'. The coverage guarantee of the sequence depends on+-- it. Out-of-range observations can drive a capital factor+-- negative, taking the construction out of the supermartingale+-- regime entirely; the function does not check for this.+--+-- >>> let s1 = update cfg s0 0.7+update :: Config -> State -> Double -> State+update Config{..} State{..} !x =+ let !y = (x - cfg_lo) / (cfg_hi - cfg_lo)+ !t = st_n + 1+ !td = fromIntegral t+ !gp1 = fromIntegral (cfg_grid + 1)+ -- sigma^2_{t-1} = (1/4 + sum_{i<=t-1} (y_i - mu_i)^2) / t+ !sig2 = (0.25 + st_sum_dev2) / td+ !lam = min trunc_c+ (sqrt (cfg_bet_num / (sig2 * td * log1p td)))+ -- built eagerly, as in 'initial': the tail is forced before+ -- consing, so the new live list is in normal form on+ -- construction.+ go [] = []+ go (Point j lp ln : ps) =+ let !m = fromIntegral j / gp1+ !d = y - m+ !lp' = lp + log1p (min lam (trunc_c / m) * d)+ !ln' = ln + log1p (negate (min lam (trunc_c / (1 - m)))+ * d)+ !rest = go ps+ in if max lp' ln' >= cfg_log_thresh+ then rest+ else Point j lp' ln' : rest+ !live = go st_live+ -- fold y into the shared statistics only now: the bet above+ -- used statistics through t-1, so predictability holds. the+ -- deviation at step t uses the current-inclusive mean mu_t.+ !sum_y' = st_sum_y + y+ !mu = (0.5 + sum_y') / (td + 1)+ !dev = y - mu+ !dev2' = st_sum_dev2 + dev * dev+ in State t sum_y' dev2' live+{-# INLINE update #-}++-- inspection -----------------------------------------------------------------++-- | The current confidence interval, in the original @[lo, hi]@+-- coordinates.+--+-- The interval spans the surviving grid candidates, widened by+-- one grid step at each end (or clamped to @lo@ \/ @hi@ at the+-- grid's edges). The widening is what makes off-grid true means+-- safe: Theorem 3 guarantees the ideal continuum survivor set is+-- an interval, so its endpoints are bracketed by the nearest+-- /rejected/ grid candidates, and reporting those sentinels+-- yields a superset of the continuum interval. Whenever the+-- result is 'Just', it therefore covers the true mean uniformly+-- over time with probability at least @1 - alpha@ -- no+-- multiplicity correction across candidates is needed, since+-- coverage concerns only the true mean's own test.+--+-- 'Nothing' means every grid candidate has been rejected: the+-- evidence has resolved the mean below the grid's resolution.+-- For a true mean lying exactly on the grid this has probability+-- at most @alpha@ (its own test must have rejected). For an+-- off-grid true mean it additionally occurs once the continuum+-- survivor interval shrinks inside a single grid cell -- a+-- quantization horizon far beyond the point where the reported+-- width is comparable to the grid spacing. Treat 'Nothing' as a+-- signal to rerun with a larger grid, not as an inference.+--+-- >>> interval cfg (initial cfg)+-- Just (0.0,1.0)+interval :: Config -> State -> Maybe (Double, Double)+interval Config{..} State{..} = case st_live of+ [] -> Nothing+ (Point j0 _ _ : ps) ->+ let !jmin = foldl' (\acc (Point j _ _) -> min acc j) j0 ps+ !jmax = foldl' (\acc (Point j _ _) -> max acc j) j0 ps+ !gp1 = fromIntegral (cfg_grid + 1)+ !w = cfg_hi - cfg_lo+ !l | jmin == 1 = cfg_lo+ | otherwise =+ cfg_lo + fromIntegral (jmin - 1) / gp1 * w+ !u | jmax == cfg_grid = cfg_hi+ | otherwise =+ cfg_lo + fromIntegral (jmax + 1) / gp1 * w+ in Just (l, u)+{-# INLINE interval #-}++-- | The number of samples consumed so far.+--+-- >>> samples s0+-- 0+samples :: State -> Int+samples = st_n+{-# INLINE samples #-}
+ lib/Numeric/Eproc/Mixture.hs view
@@ -0,0 +1,297 @@+{-# OPTIONS_HADDOCK prune #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE RecordWildCards #-}++-- |+-- Module: Numeric.Eproc.Mixture+-- Copyright: (c) 2026 Jared Tobin+-- License: MIT+-- Maintainer: Jared Tobin <jared@ppad.tech>+--+-- Uniform convex mixture of e-processes.+--+-- Given @K@ component e-processes @E^1_t, ..., E^K_t@ adapted to a+-- common filtration -- each testing (its facet of) a shared null+-- @H_0@ -- their arithmetic mean+--+-- @M_t = (E^1_t + ... + E^K_t) \/ K@+--+-- is itself an e-process with @M_0 = 1@: convex combinations+-- preserve the nonnegative-supermartingale property. By Ville's+-- inequality @P(sup_t M_t >= 1 \/ alpha) <= alpha@ under @H_0@, so a+-- level-@alpha@ test of the /combined/ null rejects when+-- @sup_t log(E^1_t + ... + E^K_t)@ crosses @log(K \/ alpha)@ -- no+-- Bonferroni correction, and strictly more powerful than one, since+-- the sum dominates the max. Use a mixture when the alternative has+-- several qualitatively different faces (a location shift, a shape+-- change, a rare-outlier channel, ...) and you want a single test+-- with power against their union.+--+-- This module does not own or update the components: they may be+-- heterogeneous (different test modules, different observation+-- transformations), so the caller steps each component itself and+-- feeds 'update' the vector of their current log e-values, as+-- reported by each module's @log_evalue@ accessor, one entry per+-- component in a fixed order.+--+-- Two preconditions are the caller's responsibility, and the+-- type-I guarantee depends on both:+--+-- 1. Each entry must be the current log e-value of a genuine+-- e-process for @H_0@, and all components must be adapted to+-- the same filtration and stepped in lockstep -- 'update' is+-- called exactly once per underlying observation, after all+-- components have absorbed it.+--+-- 2. The vector must have exactly the @K@ entries declared in+-- 'config', always in the same order.+--+-- The rejection latch is kept on the supremum of the /mixture's/+-- log-wealth. Latching (or summing) per-component suprema instead+-- would combine peaks attained at different times -- a quantity+-- that can exceed anything the mixture ever reached, silently+-- inflating the effective alpha. Ville's inequality bounds the+-- mixture's own supremum; that is the only sound latch, and it is+-- the one this module maintains.+--+-- == Example+--+-- Combine a sign test and a magnitude test running against the same+-- stream of differences @d_t@ (the shape used for two-channel+-- symmetry testing):+--+-- >>> import qualified Numeric.Eproc.Bernoulli.TwoSided as Sign+-- >>> import qualified Numeric.Eproc.Bounded as Magn+-- >>> import qualified Numeric.Eproc.Mixture as Mix+-- >>> let Right sc = Sign.config 0.5 1.0e-3 Sign.Newton+-- >>> let Right mc = Magn.config 0.0 (-1.0) 1.0 1.0e-3 Magn.Newton+-- >>> let Right xc = Mix.config 2 1.0e-3+-- >>> :{+-- let step (s, m, x) d =+-- let s' = Sign.update sc s (d > 0)+-- m' = Magn.update mc m d+-- in (s', m', Mix.update xc x+-- [Sign.log_evalue s', Magn.log_evalue m'])+-- :}+-- >>> let ds = take 400 (cycle [0.6, 0.7, -0.2, 0.8])+-- >>> let z0 = (Sign.initial sc, Magn.initial mc, Mix.initial xc)+-- >>> let (_, _, xf) = foldl' step z0 ds+-- >>> Mix.decide xc xf+-- Reject+-- >>> Mix.p_value xc xf+-- 9.482234479673792e-34++module Numeric.Eproc.Mixture (+ -- * Mixture configuration and state+ Config+ , State+ , Verdict(..)+ , ConfigError(..)++ -- * Construction+ , config+ , initial++ -- * Streaming+ , update+ , decide++ -- * Inspection+ , log_wealth+ , log_wealth_sup+ , log_evalue+ , log_evalue_sup+ , p_value+ , samples+ ) where++import Numeric.Eproc.Common (Verdict(..), ConfigError(..), finite)++-- types ----------------------------------------------------------------------++-- | Mixture configuration. Build with 'config'.+--+-- Carries the component count @K@, the significance level, the+-- precomputed rejection threshold @log(K \/ alpha)@, and @log K@+-- (the mixture log-wealth of a fresh state).+data Config = Config {+ -- ^ component count @K@+ cfg_k :: {-# UNPACK #-} !Int+ -- ^ significance level @alpha@+ , cfg_alpha :: {-# UNPACK #-} !Double+ -- ^ rejection threshold @log(K \/ alpha)@+ , cfg_log_thresh :: {-# UNPACK #-} !Double+ -- ^ @log K@+ , cfg_log_k :: {-# UNPACK #-} !Double+ }++-- | Streaming mixture state. Construct with 'initial' and fold+-- per-step component log e-value vectors through 'update'.+--+-- Tracks the current mixture log-wealth @log(sum_i E^i_t)@ and+-- its latched supremum, which is what 'decide' tests against the+-- rejection threshold.+data State = State {+ st_n :: {-# UNPACK #-} !Int -- ^ update count+ , st_log_sum :: {-# UNPACK #-} !Double -- ^ log(sum_i E^i)+ , st_sup_log_sum :: {-# UNPACK #-} !Double -- ^ sup of the above+ }++-- construction ---------------------------------------------------------------++-- | Build a 'Config' for a @K@-component uniform mixture at level+-- @alpha@.+--+-- The rejection threshold is precomputed as @log(K \/ alpha)@:+-- the mixture @M_t = (sum_i E^i_t) \/ K@ crosses @1 \/ alpha@+-- exactly when the sum crosses @K \/ alpha@.+--+-- Returns 'Left' with a 'ConfigError' on inputs outside the+-- mathematical regime: @K < 1@, or @alpha@ non-finite or outside+-- @(0, 1)@.+--+-- >>> let Right cfg = config 4 1.0e-3+config+ :: Int -- ^ component count @K@+ -> Double -- ^ significance level @alpha@+ -> Either ConfigError Config+config !k !alpha+ | k < 1 =+ Left (InvalidComponentCount k)+ | not (finite alpha && alpha > 0 && alpha < 1) =+ Left (InvalidAlpha alpha)+ | otherwise =+ let !kd = fromIntegral k+ in Right Config {+ cfg_k = k+ , cfg_alpha = alpha+ , cfg_log_thresh = log (kd / alpha)+ , cfg_log_k = log kd+ }+{-# INLINE config #-}++-- | The initial 'State' for a fresh mixture.+--+-- Every component starts at e-value @1@, so the mixture log-sum+-- (and its supremum) starts at @log K@.+--+-- >>> let s0 = initial cfg+initial :: Config -> State+initial Config{..} = State {+ st_n = 0+ , st_log_sum = cfg_log_k+ , st_sup_log_sum = cfg_log_k+ }+{-# INLINE initial #-}++-- streaming ------------------------------------------------------------------++-- | Fold one step's component log e-values into the running+-- 'State': computes the current mixture log-sum via a numerically+-- stable log-sum-exp and latches its supremum.+--+-- /Preconditions/ (documented in the module header, unchecked+-- here): the vector holds exactly the @K@ log e-values of+-- components adapted to a common filtration, in a fixed order,+-- with 'update' called once per underlying observation. The+-- degenerate empty vector leaves the state unchanged.+--+-- >>> let s1 = update cfg s0 [0.1, -0.2, 0.0, 0.4]+update :: Config -> State -> [Double] -> State+update _ st@State{..} les = case les of+ [] -> st+ (l : ls) ->+ let !m = foldl' max l ls+ !s = foldl' (\ !acc v -> acc + exp (v - m)) 0 les+ -- all components at e-value zero: the mixture log-sum is+ -- -Infinity, and (m +) would poison it into NaN.+ !cur | isInfinite m && m < 0 = m+ | otherwise = m + log s+ in State {+ st_n = st_n + 1+ , st_log_sum = cur+ , st_sup_log_sum = max st_sup_log_sum cur+ }+{-# INLINE update #-}++-- | Compute the current 'Verdict' from the running 'State'.+--+-- 'Reject' iff the supremum-so-far of @log(sum_i E^i_t)@ has ever+-- crossed @log(K \/ alpha)@ -- equivalently, the mixture+-- e-process @M_t@ has exceeded @1 \/ alpha@ at some point in the+-- stream so far. Under the combined @H_0@, by Ville's inequality,+-- the probability of this ever happening is at most @alpha@,+-- simultaneously over all sample sizes: peek and stop freely.+--+-- >>> decide cfg s0+-- Continue+decide :: Config -> State -> Verdict+decide Config{..} State{..}+ | st_sup_log_sum >= cfg_log_thresh = Reject+ | otherwise = Continue+{-# INLINE decide #-}++-- inspection -----------------------------------------------------------------++-- | The current mixture log-wealth @log(sum_i E^i_t)@, before+-- normalization by @K@. Not monotone; bounded above by+-- 'log_wealth_sup'. Starts at @log K@.+--+-- >>> log_wealth s0+-- 1.3862943611198906+log_wealth :: State -> Double+log_wealth = st_log_sum+{-# INLINE log_wealth #-}++-- | The supremum-so-far of @log(sum_i E^i_t)@. Monotone+-- nondecreasing; 'decide' rejects exactly when it crosses+-- @log(K \/ alpha)@. Starts at @log K@.+--+-- >>> log_wealth_sup s0+-- 1.3862943611198906+log_wealth_sup :: State -> Double+log_wealth_sup = st_sup_log_sum+{-# INLINE log_wealth_sup #-}++-- | The current log e-value of the mixture: the log of+-- @M_t = (sum_i E^i_t) \/ K@, i.e. 'log_wealth' minus @log K@,+-- normalized so a fresh state sits at @0@. This is itself a+-- component-shaped quantity: mixtures nest, so it can in turn be+-- fed to an outer mixture. Not monotone; bounded above by+-- 'log_evalue_sup'.+--+-- >>> log_evalue s0+-- 0.0+log_evalue :: Config -> State -> Double+log_evalue Config{..} State{..} = st_log_sum - cfg_log_k+{-# INLINE log_evalue #-}++-- | The supremum-so-far of the log e-value: 'log_wealth_sup' minus+-- @log K@. Monotone nondecreasing, starting at @0@; 'decide'+-- rejects exactly when it crosses @log(1 \/ alpha)@.+--+-- >>> log_evalue_sup s0+-- 0.0+log_evalue_sup :: Config -> State -> Double+log_evalue_sup Config{..} State{..} = st_sup_log_sum - cfg_log_k+{-# INLINE log_evalue_sup #-}++-- | The anytime-valid p-value: the reciprocal of the largest+-- mixture e-value attained so far. Monotone nonincreasing; under+-- the combined @H_0@, @P(exists t: p_t <= alpha) <= alpha@ for+-- every @alpha@ simultaneously. 'decide' returns 'Reject' exactly+-- when this value has reached the configured @alpha@ or below.+--+-- >>> p_value cfg s0+-- 1.0+p_value :: Config -> State -> Double+p_value cfg s = min 1 (exp (negate (log_evalue_sup cfg s)))+{-# INLINE p_value #-}++-- | The number of 'update' steps consumed so far.+--+-- >>> samples s0+-- 0+samples :: State -> Int+samples = st_n+{-# INLINE samples #-}
lib/Numeric/Eproc/Paired.hs view
@@ -64,6 +64,9 @@ -- * Inspection , log_wealth , log_wealth_sup+ , log_evalue+ , log_evalue_sup+ , p_value , samples ) where @@ -164,6 +167,39 @@ log_wealth_sup :: State -> Double log_wealth_sup (State s) = Bounded.log_wealth_sup s {-# INLINE log_wealth_sup #-}++-- | The current log e-value of the underlying bounded-mean test on+-- the differences: 'log_wealth' minus @log 2@, normalized so a+-- fresh state sits at @0@. Not monotone; bounded above by+-- 'log_evalue_sup'.+--+-- >>> log_evalue s0+-- 0.0+log_evalue :: State -> Double+log_evalue (State s) = Bounded.log_evalue s+{-# INLINE log_evalue #-}++-- | The supremum-so-far of the log e-value: 'log_wealth_sup' minus+-- @log 2@. Monotone nondecreasing, starting at @0@; 'decide'+-- rejects exactly when it crosses @log(1 \/ alpha)@.+--+-- >>> log_evalue_sup s0+-- 0.0+log_evalue_sup :: State -> Double+log_evalue_sup (State s) = Bounded.log_evalue_sup s+{-# INLINE log_evalue_sup #-}++-- | The anytime-valid p-value: the reciprocal of the largest+-- e-value attained so far. Monotone nonincreasing; under @H_0@,+-- @P(exists t: p_t <= alpha) <= alpha@ for every @alpha@+-- simultaneously. 'decide' returns 'Reject' exactly when this+-- value has reached the configured @alpha@ or below.+--+-- >>> p_value s0+-- 1.0+p_value :: State -> Double+p_value (State s) = Bounded.p_value s+{-# INLINE p_value #-} -- | The number of paired observations consumed so far. --
ppad-eproc.cabal view
@@ -1,6 +1,6 @@ cabal-version: 3.0 name: ppad-eproc-version: 0.3.0+version: 0.4.0 synopsis: Anytime-valid sequential testing via e-processes. license: MIT license-file: LICENSE@@ -11,11 +11,13 @@ tested-with: GHC == 9.10.3 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- and- two-sided Bernoulli rate tests with fixed, adaptive (aGRAPA), and- online Newton bettors.+ Anytime-valid sequential hypothesis testing and estimation for+ bounded random variables, via the e-process / betting framework of+ Waudby-Smith 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; anytime-valid p-values+ and e-values; uniform convex mixtures of e-processes; and+ time-uniform confidence sequences for bounded means. flag llvm description: Use GHC's LLVM backend.@@ -38,6 +40,8 @@ Numeric.Eproc.Bernoulli.TwoSided Numeric.Eproc.Bounded Numeric.Eproc.Common+ Numeric.Eproc.ConfSeq+ Numeric.Eproc.Mixture Numeric.Eproc.Paired build-depends: base >= 4.9 && < 5
test/Main.hs view
@@ -8,6 +8,8 @@ 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.ConfSeq as CS+import qualified Numeric.Eproc.Mixture as Mix import qualified Numeric.Eproc.Paired as P import Test.Tasty import Test.Tasty.HUnit@@ -25,6 +27,9 @@ , config_validation_tests , safety_property_tests , two_sided_bernoulli_tests+ , evalue_accessor_tests+ , mixture_tests+ , confseq_tests ] -- partial helper: tests below hardcode valid configs.@@ -624,3 +629,369 @@ vs = map (BernTS.decide cfg) sts in monotone_reject_bern_ts vs ]+++unit_pair :: QC.Gen (Double, Double)+unit_pair = (,) <$> unit_double <*> unit_double++evalue_accessor_tests :: TestTree+evalue_accessor_tests = testGroup "e-value accessors" [+ testCase "fresh states normalize to e-value 1, p-value 1" $ do+ let bcfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton)+ ncfg = ok (Bern.config 0.05 1.0e-3 Bern.Newton)+ tcfg = ok (BernTS.config 0.5 1.0e-3 BernTS.Newton)+ pcfg = ok (P.config 0.0 1.0 1.0e-3 Bounded.Newton)+ Bounded.log_evalue (Bounded.initial bcfg) @?= 0.0+ Bounded.log_evalue_sup (Bounded.initial bcfg) @?= 0.0+ Bounded.p_value (Bounded.initial bcfg) @?= 1.0+ Bern.log_evalue (Bern.initial ncfg) @?= 0.0+ Bern.p_value (Bern.initial ncfg) @?= 1.0+ BernTS.log_evalue (BernTS.initial tcfg) @?= 0.0+ BernTS.p_value (BernTS.initial tcfg) @?= 1.0+ P.log_evalue (P.initial pcfg) @?= 0.0+ P.p_value (P.initial pcfg) @?= 1.0++ , QC.testProperty "Bounded: log_evalue is log_wealth less log 2" $+ QC.forAll arb_bettor $ \b ->+ QC.forAll (QC.listOf unit_double) $ \xs ->+ let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 b)+ st = foldl' (Bounded.update cfg) (Bounded.initial cfg) xs+ in Bounded.log_evalue st == Bounded.log_wealth st - C.log2_dbl++ , QC.testProperty "Bernoulli: log_evalue coincides with log_wealth" $+ QC.forAll arb_bettor $ \b ->+ QC.forAll QC.arbitrary $ \xs ->+ let cfg = ok (Bern.config 0.05 1.0e-3 b)+ st = foldl' (Bern.update cfg) (Bern.initial cfg) (xs :: [Bool])+ in Bern.log_evalue st == Bern.log_wealth st++ , QC.testProperty "Bounded: decide agrees with p_value at alpha" $+ QC.forAll arb_bettor $ \b ->+ QC.forAll (QC.listOf unit_double) $ \xs ->+ let alpha = 0.5+ cfg = ok (Bounded.config 0.5 0.0 1.0 alpha b)+ sts = scanl (Bounded.update cfg) (Bounded.initial cfg) xs+ in all (\s -> (Bounded.decide cfg s == Bounded.Reject)+ == (Bounded.p_value s <= alpha)) sts++ , QC.testProperty "Bernoulli: decide agrees with p_value at alpha" $+ QC.forAll arb_bettor $ \b ->+ QC.forAll QC.arbitrary $ \xs ->+ let alpha = 0.5+ cfg = ok (Bern.config 0.5 alpha b)+ sts = scanl (Bern.update cfg) (Bern.initial cfg)+ (xs :: [Bool])+ in all (\s -> (Bern.decide cfg s == Bern.Reject)+ == (Bern.p_value s <= alpha)) sts++ , QC.testProperty "BernTS: decide agrees with p_value at alpha" $+ QC.forAll arb_bettor $ \b ->+ QC.forAll QC.arbitrary $ \xs ->+ let alpha = 0.5+ cfg = ok (BernTS.config 0.5 alpha b)+ sts = scanl (BernTS.update cfg) (BernTS.initial cfg)+ (xs :: [Bool])+ in all (\s -> (BernTS.decide cfg s == BernTS.Reject)+ == (BernTS.p_value s <= alpha)) sts++ , QC.testProperty "Bounded: p_value monotone nonincreasing" $+ QC.forAll arb_bettor $ \b ->+ QC.forAll (QC.listOf unit_double) $ \xs ->+ let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 b)+ sts = scanl (Bounded.update cfg) (Bounded.initial cfg) xs+ ps = map Bounded.p_value sts+ in and (zipWith (>=) ps (drop 1 ps))++ , QC.testProperty "Paired: p_value monotone nonincreasing" $+ QC.forAll arb_bettor $ \b ->+ QC.forAll (QC.listOf unit_pair) $ \ps ->+ let cfg = ok (P.config 0.0 1.0 1.0e-3 b)+ sts = scanl (P.update cfg) (P.initial cfg) ps+ pv = map P.p_value sts+ in and (zipWith (>=) pv (drop 1 pv))++ , QC.testProperty "Bounded: p_value in [0, 1], evalue below sup" $+ QC.forAll arb_bettor $ \b ->+ QC.forAll (QC.listOf unit_double) $ \xs ->+ let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 b)+ sts = scanl (Bounded.update cfg) (Bounded.initial cfg) xs+ in all (\s -> let p = Bounded.p_value s+ in p >= 0 && p <= 1 &&+ Bounded.log_evalue s+ <= Bounded.log_evalue_sup s) sts++ , QC.testProperty "Bernoulli: p_value in [0, 1], evalue below sup" $+ QC.forAll arb_bettor $ \b ->+ QC.forAll QC.arbitrary $ \xs ->+ let cfg = ok (Bern.config 0.05 1.0e-3 b)+ sts = scanl (Bern.update cfg) (Bern.initial cfg)+ (xs :: [Bool])+ in all (\s -> let p = Bern.p_value s+ in p >= 0 && p <= 1 &&+ Bern.log_evalue s+ <= Bern.log_evalue_sup s) sts+ ]++-- mixture --------------------------------------------------------------------++approx_eq :: Double -> Double -> Bool+approx_eq a b = abs (a - b) <= 1.0e-9 * max 1 (max (abs a) (abs b))++-- step a censor-style two-component hedge (sign + magnitude) over a+-- shared bernoulli stream, feeding the mixture the components'+-- current log e-values, with the early-stopping rule built in.+run_mixture+ :: Mix.Config+ -> BernTS.Config+ -> Bounded.Config+ -> Double -- ^ true bernoulli p+ -> Int -- ^ budget+ -> Gen+ -> (Mix.Verdict, Int)+run_mixture xc sc mc p budget g0 =+ go 0 g0 (BernTS.initial sc) (Bounded.initial mc) (Mix.initial xc)+ where+ go !n !g !s !m !x+ | n >= budget = (Mix.decide xc x, n)+ | otherwise = case Mix.decide xc x of+ Mix.Reject -> (Mix.Reject, n)+ Mix.Continue ->+ let (v, g') = bernoulli p g+ s' = BernTS.update sc s (v == 1.0)+ m' = Bounded.update mc m v+ x' = Mix.update xc x+ [BernTS.log_evalue s', Bounded.log_evalue m']+ in go (n + 1) g' s' m' x'++mixture_rate :: Double -> Double -> Int -> Int -> Word64 -> Double+mixture_rate alpha p budget trials seed =+ let xc = ok (Mix.config 2 alpha)+ sc = ok (BernTS.config 0.5 alpha BernTS.Newton)+ mc = ok (Bounded.config 0.5 0.0 1.0 alpha Bounded.Newton)+ gens = take trials (gen_seq (mk_gen seed))+ rejects = length+ [ () | g <- gens+ , let (v, _) = run_mixture xc sc mc p budget g+ , v == Mix.Reject ]+ in fromIntegral rejects / fromIntegral trials++mixture_tests :: TestTree+mixture_tests = testGroup "mixture" [+ testCase "fresh mixture sits at log K, p-value 1" $ do+ let cfg = ok (Mix.config 4 1.0e-3)+ s0 = Mix.initial cfg+ assertBool "log_wealth is log K" $+ approx_eq (Mix.log_wealth s0) (log 4)+ Mix.log_evalue cfg s0 @?= 0.0+ Mix.log_evalue_sup cfg s0 @?= 0.0+ Mix.p_value cfg s0 @?= 1.0+ Mix.decide cfg s0 @?= Mix.Continue++ , testCase "latch is on the mixture sup, not per-component sups" $ do+ -- two components peak at different times, each attaining log+ -- e-value 1.0. A bogus combination of per-component suprema,+ -- log_sum_exp 1 1 ~ 1.69, crosses the K = 2, alpha = 0.5+ -- threshold log 4 ~ 1.39; the mixture itself never exceeds+ -- ~1.003 and must not reject.+ let cfg = ok (Mix.config 2 0.5)+ s1 = Mix.update cfg (Mix.initial cfg) [1.0, -5.0]+ s2 = Mix.update cfg s1 [-5.0, 1.0]+ Mix.decide cfg s2 @?= Mix.Continue+ assertBool "mixture sup below threshold" $+ Mix.log_wealth_sup s2 < log 4+ assertBool "per-component-sup combination would cross" $+ C.log_sum_exp 1.0 1.0 >= log 4++ , testCase "empty update vector is a no-op" $ do+ let cfg = ok (Mix.config 2 1.0e-3)+ s0 = Mix.initial cfg+ s1 = Mix.update cfg s0 []+ Mix.samples s1 @?= 0+ Mix.log_wealth s1 @?= Mix.log_wealth s0++ , testCase "config validation" $ do+ let assert_left :: Either C.ConfigError Mix.Config -> Assertion+ assert_left e = case e of+ Left _ -> pure ()+ Right _ -> assertFailure "expected Left"+ assert_left (Mix.config 0 0.05)+ assert_left (Mix.config (-3) 0.05)+ assert_left (Mix.config 4 0.0)+ assert_left (Mix.config 4 1.5)+ assert_left (Mix.config 4 (0 / 0))++ , QC.testProperty "K identical components track the component" $+ QC.forAll (QC.choose (1, 6)) $ \k ->+ QC.forAll (QC.listOf unit_double) $ \xs ->+ let bcfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton)+ xcfg = ok (Mix.config k 1.0e-3)+ sts = drop 1 (scanl (Bounded.update bcfg)+ (Bounded.initial bcfg) xs)+ les = map Bounded.log_evalue sts+ mix = foldl'+ (\acc l -> Mix.update xcfg acc (replicate k l))+ (Mix.initial xcfg) les+ cfin = foldl' (Bounded.update bcfg) (Bounded.initial bcfg) xs+ in approx_eq (Mix.log_evalue xcfg mix)+ (Bounded.log_evalue cfin)+ && approx_eq (Mix.log_evalue_sup xcfg mix)+ (Bounded.log_evalue_sup cfin)++ , QC.testProperty "decide agrees with p_value at alpha" $+ QC.forAll (QC.choose (1, 6)) $ \k ->+ QC.forAll (QC.listOf (QC.vectorOf k (QC.choose (-5, 5)))) $ \vs ->+ let alpha = 0.5+ cfg = ok (Mix.config k alpha)+ sts = scanl (Mix.update cfg) (Mix.initial cfg) vs+ in all (\s -> (Mix.decide cfg s == Mix.Reject)+ == (Mix.p_value cfg s <= alpha)) sts++ , QC.testProperty "sup monotone nondecreasing, verdict latched" $+ QC.forAll (QC.choose (1, 6)) $ \k ->+ QC.forAll (QC.listOf (QC.vectorOf k (QC.choose (-5, 5)))) $ \vs ->+ let cfg = ok (Mix.config k 0.5)+ sts = scanl (Mix.update cfg) (Mix.initial cfg) vs+ sups = map Mix.log_wealth_sup sts+ in and (zipWith (<=) sups (drop 1 sups))+ && monotone_reject_bounded (map (Mix.decide cfg) sts)++ , testCase "FPR under H_0 within slack (sign + magnitude hedge)" $ do+ let rate = mixture_rate 0.05 0.5 2000 200 424242+ assertBool ("FPR " ++ show rate ++ " exceeded slack") $+ rate <= 0.08++ , testCase "power against p = 0.7 (sign + magnitude hedge)" $ do+ let rate = mixture_rate 1.0e-3 0.7 5000 100 434343+ assertBool ("power " ++ show rate ++ " too low") $+ rate >= 0.95+ ]+-- confidence sequences -------------------------------------------------------+-- a finite stream of bernoulli(p) samples.+cs_stream :: Double -> Int -> Gen -> [Double]+cs_stream !p n g0 = go n g0+ where+ go 0 _ = []+ go !k !g =+ let (x, g') = bernoulli p g+ in x : go (k - 1) g'++-- do the intervals nest: each contained in its predecessor, with+-- Nothing (empty) absorbing?+cs_nested :: [Maybe (Double, Double)] -> Bool+cs_nested ivs = and (zipWith shrink ivs (drop 1 ivs))+ where+ shrink (Just (l1, u1)) (Just (l2, u2)) = l2 >= l1 && u2 <= u1+ shrink (Just _) Nothing = True+ shrink Nothing Nothing = True+ shrink Nothing (Just _) = False++-- fraction of trials in which the true mean ever escapes the running+-- interval (or the interval goes empty), checked after every+-- observation.+cs_miscoverage_rate+ :: CS.Config+ -> Double -- ^ true mean+ -> Int -- ^ budget per trial+ -> Int -- ^ number of trials+ -> Word64 -- ^ seed+ -> Double+cs_miscoverage_rate cfg p budget trials seed =+ let gens = take trials (gen_seq (mk_gen seed))+ misses = length [ () | g <- gens, cs_trial_missed g ]+ in fromIntegral misses / fromIntegral trials+ where+ cs_trial_missed g0 = go budget g0 (CS.initial cfg)+ where+ go !k !g !st+ | k == 0 = False+ | otherwise =+ let (x, g') = bernoulli p g+ st' = CS.update cfg st x+ in case CS.interval cfg st' of+ Nothing -> True+ Just (l, u)+ | p < l || p > u -> True+ | otherwise -> go (k - 1) g' st'++confseq_tests :: TestTree+confseq_tests = testGroup "confidence sequences" [+ testCase "initial interval is the full range" $ do+ let cfg = ok (CS.config 0.0 1.0 0.05 100)+ CS.interval cfg (CS.initial cfg) @?= Just (0.0, 1.0)+ , testCase "intervals nest along a deterministic stream" $ do+ let cfg = ok (CS.config 0.0 1.0 0.05 50)+ xs = take 500 (cycle [1.0, 1.0, 0.0, 1.0])+ sts = scanl (CS.update cfg) (CS.initial cfg) xs+ ivs = map (CS.interval cfg) sts+ assertBool "nesting violated" (cs_nested ivs)+ -- the stream has empirical mean 0.75; the final interval must+ -- be a strict refinement of the initial one.+ case (ivs, reverse ivs) of+ (iv0 : _, ivn : _) -> assertBool "no shrinkage" (iv0 /= ivn)+ _ -> assertFailure "no intervals"+ , QC.testProperty "intervals nest along any admissible stream" $+ QC.forAll (QC.listOf unit_double) $ \xs ->+ let cfg = ok (CS.config 0.0 1.0 0.05 25)+ sts = scanl (CS.update cfg) (CS.initial cfg) xs+ in cs_nested (map (CS.interval cfg) sts)+ , testCase "coverage: off-grid Bernoulli(0.437) at alpha = 0.05" $ do+ let cfg = ok (CS.config 0.0 1.0 0.05 100)+ rate = cs_miscoverage_rate cfg 0.437 1500 200 991199+ -- expected miscoverage <= 0.05; allow up to 0.08 slack for+ -- sampling variability over 200 trials.+ assertBool ("miscoverage " ++ show rate ++ " exceeded slack") $+ rate <= 0.08+ , testCase "consistency: Bernoulli(0.3) interval shrinks onto mean" $ do+ let cfg = ok (CS.config 0.0 1.0 1.0e-3 200)+ xs = cs_stream 0.3 5000 (mk_gen 424242)+ st = foldl' (CS.update cfg) (CS.initial cfg) xs+ case CS.interval cfg st of+ Nothing -> assertFailure "interval empty"+ Just (l, u) -> do+ assertBool ("interval " ++ show (l, u) ++ " misses mean") $+ l <= 0.3 && 0.3 <= u+ assertBool ("width " ++ show (u - l) ++ " too wide") $+ u - l < 0.2+ , testCase "affine: mean recovered on [-5, 5]" $ do+ -- x = 4 w.p. 0.7, x = -4 w.p. 0.3: true mean 1.6, interior+ -- to the sample bounds and asymmetric about zero.+ let cfg = ok (CS.config (-5.0) 5.0 0.05 100)+ xs = [ if b == 1.0 then 4.0 else (-4.0)+ | b <- cs_stream 0.7 3000 (mk_gen 232323) ]+ st = foldl' (CS.update cfg) (CS.initial cfg) xs+ case CS.interval cfg st of+ Nothing -> assertFailure "interval empty"+ Just (l, u) -> do+ assertBool ("interval " ++ show (l, u) ++ " misses mean") $+ l <= 1.6 && 1.6 <= u+ assertBool ("interval " ++ show (l, u) ++ " not refined") $+ l > -5.0 && u < 5.0+ , testCase "config: grid size 0 rejected" $+ assertLeftCS (CS.config 0.0 1.0 0.05 0)+ , testCase "config: negative grid size rejected" $+ assertLeftCS (CS.config 0.0 1.0 0.05 (-3))+ , testCase "config: alpha out of range rejected" $ do+ assertLeftCS (CS.config 0.0 1.0 0.0 100)+ assertLeftCS (CS.config 0.0 1.0 1.5 100)+ , testCase "config: lo >= hi rejected" $+ assertLeftCS (CS.config 1.0 0.0 0.05 100)+ , testCase "config: non-finite inputs rejected" $ do+ let nan = 0 / 0 :: Double+ pInf = 1 / 0 :: Double+ assertLeftCS (CS.config nan 1.0 0.05 100)+ assertLeftCS (CS.config 0.0 pInf 0.05 100)+ assertLeftCS (CS.config 0.0 1.0 nan 100)+ , QC.testProperty "interval endpoints well-formed on any stream" $+ QC.forAll (QC.listOf unit_double) $ \xs ->+ let cfg = ok (CS.config 0.0 1.0 0.05 25)+ st = foldl' (CS.update cfg) (CS.initial cfg) xs+ in case CS.interval cfg st of+ Nothing -> True+ Just (l, u) ->+ finite l && finite u && 0 <= l && l <= u && u <= 1+ ]+ where+ assertLeftCS :: Either C.ConfigError a -> Assertion+ assertLeftCS e = case e of+ Left _ -> pure ()+ Right _ -> assertFailure "expected Left"