diff --git a/CHANGELOG b/CHANGELOG --- a/CHANGELOG +++ b/CHANGELOG @@ -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 diff --git a/bench/Main.hs b/bench/Main.hs --- a/bench/Main.hs +++ b/bench/Main.hs @@ -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 ] diff --git a/bench/Weight.hs b/bench/Weight.hs --- a/bench/Weight.hs +++ b/bench/Weight.hs @@ -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 diff --git a/lib/Numeric/Eproc/Bernoulli.hs b/lib/Numeric/Eproc/Bernoulli.hs --- a/lib/Numeric/Eproc/Bernoulli.hs +++ b/lib/Numeric/Eproc/Bernoulli.hs @@ -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. -- diff --git a/lib/Numeric/Eproc/Bernoulli/TwoSided.hs b/lib/Numeric/Eproc/Bernoulli/TwoSided.hs --- a/lib/Numeric/Eproc/Bernoulli/TwoSided.hs +++ b/lib/Numeric/Eproc/Bernoulli/TwoSided.hs @@ -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. -- diff --git a/lib/Numeric/Eproc/Bounded.hs b/lib/Numeric/Eproc/Bounded.hs --- a/lib/Numeric/Eproc/Bounded.hs +++ b/lib/Numeric/Eproc/Bounded.hs @@ -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. -- diff --git a/lib/Numeric/Eproc/Common.hs b/lib/Numeric/Eproc/Common.hs --- a/lib/Numeric/Eproc/Common.hs +++ b/lib/Numeric/Eproc/Common.hs @@ -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 diff --git a/lib/Numeric/Eproc/ConfSeq.hs b/lib/Numeric/Eproc/ConfSeq.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Eproc/ConfSeq.hs @@ -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 +-- +-- 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 #-} diff --git a/lib/Numeric/Eproc/Mixture.hs b/lib/Numeric/Eproc/Mixture.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Eproc/Mixture.hs @@ -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 +-- +-- 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 #-} diff --git a/lib/Numeric/Eproc/Paired.hs b/lib/Numeric/Eproc/Paired.hs --- a/lib/Numeric/Eproc/Paired.hs +++ b/lib/Numeric/Eproc/Paired.hs @@ -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. -- diff --git a/ppad-eproc.cabal b/ppad-eproc.cabal --- a/ppad-eproc.cabal +++ b/ppad-eproc.cabal @@ -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 diff --git a/test/Main.hs b/test/Main.hs --- a/test/Main.hs +++ b/test/Main.hs @@ -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"