eproc

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

commit dcd9754da35f9ded3faeef878fbf234497ed11c2
parent 8b73a8dc0a86c40d9e7845ce275580f9e87d083a
Author: Jared Tobin <jared@jtobin.io>
Date:   Wed,  3 Jun 2026 11:28:29 -0230

rename Mean to Bounded

Module Numeric.Eproc.Mean -- which provides the two-sided bounded-mean
test, not just mean calculation -- becomes Numeric.Eproc.Bounded.
'Bounded' parallels 'Paired' as a noun describing the test family.

Import alias changes from 'M' to the spelled-out 'Bounded' (the
natural single-letter 'B' is already in use by Bettor); call sites in
Paired, tests, bench, and the README updated.

Internal helpers and labels updated to match: local test helper
'run_mean_bernoulli' -> 'run_bounded_bernoulli'; benchmark bgroup
labels 'Mean.update ...' / 'Mean.decide' -> 'Bounded.update ...' /
'Bounded.decide'; README's captured benchmark output relabeled.

No semantic / API changes; tests pass.

Diffstat:
MREADME.md | 28++++++++++++++--------------
Mbench/Main.hs | 44++++++++++++++++++++++----------------------
Mbench/Weight.hs | 44++++++++++++++++++++++----------------------
Mlib/Numeric/Eproc/Bettor.hs | 2+-
Alib/Numeric/Eproc/Bounded.hs | 310+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dlib/Numeric/Eproc/Mean.hs | 310-------------------------------------------------------------------------------
Mlib/Numeric/Eproc/Paired.hs | 24++++++++++++------------
Mppad-eproc.cabal | 2+-
Mtest/Main.hs | 68++++++++++++++++++++++++++++++++++----------------------------------
9 files changed, 416 insertions(+), 416 deletions(-)

diff --git a/README.md b/README.md @@ -18,30 +18,30 @@ A sample GHCi session: ``` > -- import qualified > import qualified Numeric.Eproc.Bettor as B - > import qualified Numeric.Eproc.Mean as M + > import qualified Numeric.Eproc.Bounded as Bounded > > -- test H_0: E[X] = 0.5 for samples in [0, 1] at alpha = 1e-3, > -- with the ONS bettor - > let cfg = M.config 0.5 0.0 1.0 1.0e-3 B.Ons + > let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Ons > > -- streaming interface: 'initial' then fold observations through 'update' - > let s0 = M.initial cfg + > let s0 = Bounded.initial cfg > let xs = [1, 1, 0, 1, 1, 0, 1, 1, 1, 1] -- mean 0.8, drifts from H_0 - > let s10 = foldl (M.update cfg) s0 xs + > let s10 = foldl (Bounded.update cfg) s0 xs > > -- inspect wealth and verdict at any point - > M.samples s10 + > Bounded.samples s10 10 - > M.log_wealth s10 + > Bounded.log_wealth s10 0.7182493502552663 - > M.decide cfg s10 + > Bounded.decide cfg s10 Continue > > -- with enough evidence the test rejects - > let s300 = foldl (M.update cfg) s0 (concat (replicate 30 xs)) - > M.log_wealth s300 + > let s300 = foldl (Bounded.update cfg) s0 (concat (replicate 30 xs)) + > Bounded.log_wealth s300 53.092214534054165 - > M.decide cfg s300 + > Bounded.decide cfg s300 Reject ``` @@ -62,25 +62,25 @@ Current benchmark figures on an M4 Silicon MacBook Air look like (use `cabal bench` to run the benchmark suite): ``` - benchmarking Mean.update (one step)/ons + benchmarking Bounded.update (one step)/ons time 13.05 ns (12.95 ns .. 13.17 ns) 1.000 R² (0.999 R² .. 1.000 R²) mean 13.03 ns (12.95 ns .. 13.15 ns) std dev 314.0 ps (248.3 ps .. 422.3 ps) - benchmarking Mean.update (1000-sample fold)/fixed + benchmarking Bounded.update (1000-sample fold)/fixed time 4.840 μs (4.819 μs .. 4.867 μs) 1.000 R² (1.000 R² .. 1.000 R²) mean 4.828 μs (4.817 μs .. 4.847 μs) std dev 44.90 ns (30.94 ns .. 61.54 ns) - benchmarking Mean.update (1000-sample fold)/agrapa + benchmarking Bounded.update (1000-sample fold)/agrapa time 15.67 μs (15.66 μs .. 15.69 μs) 1.000 R² (1.000 R² .. 1.000 R²) mean 15.67 μs (15.65 μs .. 15.69 μs) std dev 63.74 ns (55.65 ns .. 75.07 ns) - benchmarking Mean.update (1000-sample fold)/ons + benchmarking Bounded.update (1000-sample fold)/ons time 14.43 μs (14.42 μs .. 14.44 μs) 1.000 R² (1.000 R² .. 1.000 R²) mean 14.43 μs (14.42 μs .. 14.44 μs) diff --git a/bench/Main.hs b/bench/Main.hs @@ -5,16 +5,16 @@ module Main where import Control.DeepSeq import qualified Numeric.Eproc.Bettor as B -import qualified Numeric.Eproc.Mean as M +import qualified Numeric.Eproc.Bounded as Bounded import qualified Numeric.Eproc.Paired as P import Criterion.Main -- all relevant fields are strict (and UNPACK'd for the doubles), so -- WHNF == NF for these types. orphan instances keep the library API -- untouched. -instance NFData M.State where rnf !_ = () +instance NFData Bounded.State where rnf !_ = () instance NFData P.State where rnf !_ = () -instance NFData M.Verdict where rnf !_ = () +instance NFData Bounded.Verdict where rnf !_ = () main :: IO () main = defaultMain [ @@ -26,35 +26,35 @@ main = defaultMain [ update :: Benchmark update = - let !cfg_f = M.config 0.5 0.0 1.0 1.0e-3 (B.Fixed 0.5) - !cfg_a = M.config 0.5 0.0 1.0 1.0e-3 B.Agrapa - !cfg_o = M.config 0.5 0.0 1.0 1.0e-3 B.Ons - !st_f = M.initial cfg_f - !st_a = M.initial cfg_a - !st_o = M.initial cfg_o + let !cfg_f = Bounded.config 0.5 0.0 1.0 1.0e-3 (B.Fixed 0.5) + !cfg_a = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Agrapa + !cfg_o = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Ons + !st_f = Bounded.initial cfg_f + !st_a = Bounded.initial cfg_a + !st_o = Bounded.initial cfg_o !x = 0.7 - in bgroup "Mean.update (one step)" [ - bench "fixed" $ nf (M.update cfg_f st_f) x - , bench "agrapa" $ nf (M.update cfg_a st_a) x - , bench "ons" $ nf (M.update cfg_o st_o) x + in bgroup "Bounded.update (one step)" [ + bench "fixed" $ nf (Bounded.update cfg_f st_f) x + , bench "agrapa" $ nf (Bounded.update cfg_a st_a) x + , bench "ons" $ nf (Bounded.update cfg_o st_o) x ] decide :: Benchmark decide = - let !cfg = M.config 0.5 0.0 1.0 1.0e-3 B.Ons - !st = M.initial cfg - in bgroup "Mean.decide" [ - bench "initial state" $ nf (M.decide cfg) st + let !cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Ons + !st = Bounded.initial cfg + in bgroup "Bounded.decide" [ + bench "initial state" $ nf (Bounded.decide cfg) st ] stream :: Benchmark stream = let !xs = force (take 1000 (cycle [0.3, 0.7])) - !cfg_f = M.config 0.5 0.0 1.0 1.0e-3 (B.Fixed 0.5) - !cfg_a = M.config 0.5 0.0 1.0 1.0e-3 B.Agrapa - !cfg_o = M.config 0.5 0.0 1.0 1.0e-3 B.Ons - run_m cfg = foldl' (M.update cfg) (M.initial cfg) - in bgroup "Mean.update (1000-sample fold)" [ + !cfg_f = Bounded.config 0.5 0.0 1.0 1.0e-3 (B.Fixed 0.5) + !cfg_a = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Agrapa + !cfg_o = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Ons + run_m cfg = foldl' (Bounded.update cfg) (Bounded.initial cfg) + in bgroup "Bounded.update (1000-sample fold)" [ bench "fixed" $ nf (run_m cfg_f) xs , bench "agrapa" $ nf (run_m cfg_a) xs , bench "ons" $ nf (run_m cfg_o) xs diff --git a/bench/Weight.hs b/bench/Weight.hs @@ -5,13 +5,13 @@ module Main where import Control.DeepSeq import qualified Numeric.Eproc.Bettor as B -import qualified Numeric.Eproc.Mean as M +import qualified Numeric.Eproc.Bounded as Bounded import qualified Numeric.Eproc.Paired as P import Weigh -instance NFData M.State where rnf !_ = () +instance NFData Bounded.State where rnf !_ = () instance NFData P.State where rnf !_ = () -instance NFData M.Verdict where rnf !_ = () +instance NFData Bounded.Verdict where rnf !_ = () -- note that 'weigh' doesn't work properly in a repl main :: IO () @@ -23,32 +23,32 @@ main = mainWith $ do update :: Weigh () update = - let !cfg_f = M.config 0.5 0.0 1.0 1.0e-3 (B.Fixed 0.5) - !cfg_a = M.config 0.5 0.0 1.0 1.0e-3 B.Agrapa - !cfg_o = M.config 0.5 0.0 1.0 1.0e-3 B.Ons - !st_f = M.initial cfg_f - !st_a = M.initial cfg_a - !st_o = M.initial cfg_o - in wgroup "Mean.update (one step)" $ do - func "fixed" (M.update cfg_f st_f) 0.7 - func "agrapa" (M.update cfg_a st_a) 0.7 - func "ons" (M.update cfg_o st_o) 0.7 + let !cfg_f = Bounded.config 0.5 0.0 1.0 1.0e-3 (B.Fixed 0.5) + !cfg_a = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Agrapa + !cfg_o = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Ons + !st_f = Bounded.initial cfg_f + !st_a = Bounded.initial cfg_a + !st_o = Bounded.initial cfg_o + in wgroup "Bounded.update (one step)" $ do + func "fixed" (Bounded.update cfg_f st_f) 0.7 + func "agrapa" (Bounded.update cfg_a st_a) 0.7 + func "ons" (Bounded.update cfg_o st_o) 0.7 decide :: Weigh () decide = - let !cfg = M.config 0.5 0.0 1.0 1.0e-3 B.Ons - !st = M.initial cfg - in wgroup "Mean.decide" $ do - func "initial state" (M.decide cfg) st + let !cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Ons + !st = Bounded.initial cfg + in wgroup "Bounded.decide" $ do + func "initial state" (Bounded.decide cfg) st stream :: Weigh () stream = let !xs = force (take 1000 (cycle [0.3, 0.7])) - !cfg_f = M.config 0.5 0.0 1.0 1.0e-3 (B.Fixed 0.5) - !cfg_a = M.config 0.5 0.0 1.0 1.0e-3 B.Agrapa - !cfg_o = M.config 0.5 0.0 1.0 1.0e-3 B.Ons - run_m cfg = foldl' (M.update cfg) (M.initial cfg) - in wgroup "Mean.update (1000-sample fold)" $ do + !cfg_f = Bounded.config 0.5 0.0 1.0 1.0e-3 (B.Fixed 0.5) + !cfg_a = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Agrapa + !cfg_o = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Ons + run_m cfg = foldl' (Bounded.update cfg) (Bounded.initial cfg) + in wgroup "Bounded.update (1000-sample fold)" $ do func "fixed" (run_m cfg_f) xs func "agrapa" (run_m cfg_a) xs func "ons" (run_m cfg_o) xs diff --git a/lib/Numeric/Eproc/Bettor.hs b/lib/Numeric/Eproc/Bettor.hs @@ -34,7 +34,7 @@ module Numeric.Eproc.Bettor ( -- For 'Agrapa' and 'Ons', a per-direction safe-bet ceiling -- @lambda_max@ is derived from the sample bounds supplied to the -- surrounding test configuration (e.g. --- 'Numeric.Eproc.Mean.config') -- bets get clipped to +-- 'Numeric.Eproc.Bounded.config') -- bets get clipped to -- @[0, lambda_max]@ so that the wealth factor @1 + lambda * z@ -- stays nonnegative for every admissible observation. -- diff --git a/lib/Numeric/Eproc/Bounded.hs b/lib/Numeric/Eproc/Bounded.hs @@ -0,0 +1,310 @@ +{-# OPTIONS_HADDOCK prune #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE RecordWildCards #-} + +-- | +-- Module: Numeric.Eproc.Bounded +-- Copyright: (c) 2026 Jared Tobin +-- License: MIT +-- Maintainer: Jared Tobin <jared@ppad.tech> +-- +-- Two-sided bounded-mean anytime-valid test. +-- +-- For samples @x_t@ in @[lo, hi]@, tests @H_0: E[x] = m@ against +-- @H_1: E[x] /= m@. +-- +-- Internally two one-sided e-processes are run in parallel: a +-- /positive-direction/ process betting against the alternative +-- @E[x] > m@ (using centred observations @z = x - m@), and a +-- /negative-direction/ process betting against @E[x] < m@ (using +-- @-z@). Each maintains its own log-wealth and bettor state. The +-- test rejects when either side's wealth crosses @2 \/ alpha@; the +-- factor of 2 is the Bonferroni adjustment for the two-sided union. +-- +-- The test is /anytime-valid/: under @H_0@ the wealth process is a +-- nonnegative supermartingale, so by Ville's inequality the +-- probability of ever crossing the threshold is at most @alpha@, +-- regardless of when the user decides to stop streaming samples. + +module Numeric.Eproc.Bounded ( + -- * Test configuration and state + Config + , State + , Verdict(..) + + -- * Construction + , config + , initial + + -- * Streaming + , update + , decide + + -- * Inspection + , log_wealth + , samples + ) where + +import GHC.Exts (Double(D#)) +import Numeric.Eproc.Bettor + +-- types ---------------------------------------------------------------------- + +-- | Test outcome at the current sample count. +-- +-- 'Reject' means the wealth process has crossed the Bonferroni +-- threshold, so @H_0@ is rejected at level @alpha@. 'Continue' +-- means there is not yet enough evidence; collect more samples (or +-- stop and report no rejection -- the type-I error guarantee holds +-- for /any/ stopping rule). +data Verdict = + Reject + | Continue + deriving (Eq, Show) + +-- per-direction bettor state. one constructor per 'Bettor' alternative; +-- the constructor used in a given 'State' matches the 'Bettor' chosen +-- in the enclosing 'Config'. +data BetState = + SFixed + | SAgrapa + {-# UNPACK #-} !Double -- sum of z (centred observation) + {-# UNPACK #-} !Double -- sum of z^2 (for online variance) + {-# UNPACK #-} !Int -- count + | SOns + {-# UNPACK #-} !Double -- current bet lambda + {-# UNPACK #-} !Double -- running sum of per-step squared gradients + +-- | Bounded-mean test configuration. Build with 'config'. +-- +-- Carries the bettor strategy, the null mean, the significance +-- level, the precomputed Bonferroni-adjusted log-wealth threshold, +-- and the per-direction safe-bet ceilings (see 'config' for how +-- the latter are derived from the sample bounds). +data Config = Config { + cfg_bettor :: !Bettor -- ^ bettor strategy + , cfg_lam_max_pos :: {-# UNPACK #-} !Double -- ^ pos-direction safe-bet ceiling + , cfg_lam_max_neg :: {-# UNPACK #-} !Double -- ^ neg-direction safe-bet ceiling + , cfg_null_mean :: {-# UNPACK #-} !Double -- ^ null mean @m@ + , cfg_alpha :: {-# UNPACK #-} !Double -- ^ significance level @alpha@ + , cfg_log_thresh :: {-# UNPACK #-} !Double -- ^ rejection threshold @log(2 \/ alpha)@ + } + +-- | Streaming test state. Construct with 'initial' and fold +-- observations through 'update'. +-- +-- The two log-wealth fields track the running log-wealth of the +-- positive- and negative-direction e-processes separately; +-- 'decide' compares each to the threshold and 'log_wealth' returns +-- the larger of the two. The per-direction bettor states carry +-- whatever the chosen 'Bettor' needs (running sums, current bet, +-- etc.). +data State = State { + st_n :: {-# UNPACK #-} !Int -- ^ sample count + , st_log_w_pos :: {-# UNPACK #-} !Double -- ^ log-wealth, pos-direction process + , st_log_w_neg :: {-# UNPACK #-} !Double -- ^ log-wealth, neg-direction process + , st_bet_pos :: !BetState -- ^ bettor state, pos-direction + , st_bet_neg :: !BetState -- ^ bettor state, neg-direction + } + +-- internal ------------------------------------------------------------------- + +-- floor for the wealth factor before taking a log; keeps the running +-- log-wealth finite when a step pushes the factor to (or below) zero. +-- NB. written via MagicHash because the fractional literal '1.0e-300' +-- compiles as 'fromRational (1.0e-300 :: Rational)', and GHC does +-- not constant-fold the conversion -- leaving a per-step +-- '$wrationalToDouble' call in the worker. +tiny :: Double +tiny = D# 1.0e-300## +{-# INLINE tiny #-} + +-- per-bettor initial state. +init_bet :: Bettor -> BetState +init_bet b = case b of + Fixed _ -> SFixed + Agrapa -> SAgrapa 0 0 0 + Ons -> SOns 0 1.0e-6 -- small acc seed avoids div-by-zero on first step +{-# INLINE init_bet #-} + +-- compute the next bet 'lambda' from the bettor and its current +-- state; 'lam_max' is the direction-specific safety bound. for +-- Agrapa we form a Kelly-style plug-in from the running sample mean +-- and variance; for Ons the bet is just the last lambda chosen by the +-- Newton step (updated during 'step_bet'). +bet_lambda :: Bettor -> Double -> BetState -> Double +bet_lambda b !lam_max !s = case b of + Fixed lam -> lam + Agrapa -> case s of + SAgrapa !sm !sm2 !n + | n == 0 -> 0 + | otherwise -> + let !nd = fromIntegral n + !mu = sm / nd + !mu2 = mu * mu + !var = max 0 (sm2 / nd - mu2) + !den = var + mu2 + !raw = if den == 0 then 0 else mu / den + in max 0 (min lam_max raw) + _ -> 0 + Ons -> case s of + SOns !lam _ -> lam + _ -> 0 +{-# INLINE bet_lambda #-} + +-- update bettor state with newly observed centred value 'z'. for +-- Agrapa this is just accumulating sums; for Ons we take one Newton +-- step on the per-step log-wealth loss '-log(1 + lambda * z)', +-- accumulating squared gradients for adaptive scaling. +step_bet :: Bettor -> Double -> BetState -> Double -> BetState +step_bet b !lam_max !s !z = case b of + Fixed _ -> SFixed + Agrapa -> case s of + SAgrapa !sm !sm2 !n -> SAgrapa (sm + z) (sm2 + z * z) (n + 1) + _ -> SAgrapa z (z * z) 1 + Ons -> case s of + SOns !lam !acc -> + let !denom = 1 + lam * z + !g = if denom == 0 then 0 else negate z / denom + !acc' = acc + g * g + !lam' = lam - g / acc' + !clp = max 0 (min lam_max lam') + in SOns clp acc' + _ -> SOns 0 1.0e-6 +{-# INLINE step_bet #-} + +-- construction --------------------------------------------------------------- + +-- | Build a 'Config' for the bounded-mean test. +-- +-- Each per-direction safe-bet ceiling @lambda_max@ is set so that +-- the wealth factor stays nonnegative for every admissible +-- observation: +-- +-- * The positive-direction factor is @1 + lambda_p * (x - m)@. +-- Since @x@ can dip to @lo@, @x - m@ can reach @lo - m@ (the +-- most negative value), so we need +-- @lambda_p <= 1 \/ (m - lo)@. The ceiling stored is half this +-- to leave numerical margin -- the WSR safety recommendation. +-- +-- * The negative-direction factor is @1 - lambda_n * (x - m)@. +-- Since @x@ can rise to @hi@, @x - m@ can reach @hi - m@, so we +-- need @lambda_n <= 1 \/ (hi - m)@; again the ceiling is set to +-- half this. +-- +-- The log-wealth rejection threshold is precomputed as +-- @log(2 \/ alpha)@; the 2 is the Bonferroni union-bound +-- adjustment for the two one-sided e-processes. +-- +-- >>> import qualified Numeric.Eproc.Bettor as B +-- >>> let cfg = config 0.5 0.0 1.0 1.0e-3 B.Ons +config + :: Double -- ^ null mean @m@ + -> Double -- ^ sample lower bound @lo@ + -> Double -- ^ sample upper bound @hi@ + -> Double -- ^ significance level @alpha@ + -> Bettor -- ^ bettor strategy + -> Config +config !m !lo !hi !alpha !b = Config { + cfg_bettor = b + , cfg_lam_max_pos = 0.5 / (m - lo) + , cfg_lam_max_neg = 0.5 / (hi - m) + , cfg_null_mean = m + , cfg_alpha = alpha + , cfg_log_thresh = log (2 / alpha) + } +{-# INLINE config #-} + +-- | The initial 'State' for a fresh streaming test. +-- +-- Both directional log-wealths start at @0@ (i.e., wealth @1@) and +-- both bettors start in the per-strategy initial state appropriate +-- for the 'Bettor' chosen in the 'Config'. +-- +-- >>> let s0 = initial cfg +initial :: Config -> State +initial Config{..} = + let !s0 = init_bet cfg_bettor + in State { + st_n = 0 + , st_log_w_pos = 0 + , st_log_w_neg = 0 + , st_bet_pos = s0 + , st_bet_neg = s0 + } +{-# INLINE initial #-} + +-- streaming ------------------------------------------------------------------ + +-- | Fold one observation into the running 'State'. +-- +-- Computes the centred observation @z = x - m@, queries the two +-- directional bettors for their predictable bets, accumulates +-- per-direction log-wealth via +-- +-- @log_w' = log_w + log (1 + lambda * z)@ +-- +-- (with the symmetric @-lambda@ for the negative direction), and +-- then steps the bettor states given the newly observed @z@. The +-- per-step wealth factor is floored at a tiny positive value to +-- keep the log finite when a marginal bet drives the factor to (or +-- below) zero. +-- +-- >>> let s1 = update cfg s0 0.7 +update :: Config -> State -> Double -> State +update Config{..} State{..} !x = + let !z = x - cfg_null_mean + !lam_p = bet_lambda cfg_bettor cfg_lam_max_pos st_bet_pos + !lam_n = bet_lambda cfg_bettor cfg_lam_max_neg st_bet_neg + !fac_p = 1 + lam_p * z + !fac_n = 1 - lam_n * z + !logw_p = st_log_w_pos + log (max tiny fac_p) + !logw_n = st_log_w_neg + log (max tiny fac_n) + !sp = step_bet cfg_bettor cfg_lam_max_pos st_bet_pos z + !sn = step_bet cfg_bettor cfg_lam_max_neg st_bet_neg (negate z) + in State (st_n + 1) logw_p logw_n sp sn +{-# INLINE update #-} + +-- | Compute the current 'Verdict' from the running 'State'. +-- +-- 'Reject' iff either directional log-wealth has crossed the +-- Bonferroni-adjusted threshold @log(2 \/ alpha)@; equivalently, +-- the wealth process on either side has exceeded @2 \/ alpha@. +-- Under @H_0@, by Ville's inequality, the probability of this ever +-- happening is at most @alpha@ -- and crucially this bound holds +-- at /every/ sample size simultaneously, so the user is free to +-- peek at the verdict as often as they like and stop on the first +-- 'Reject'. +-- +-- >>> decide cfg s0 +-- Continue +decide :: Config -> State -> Verdict +decide Config{..} State{..} + | st_log_w_pos >= cfg_log_thresh = Reject + | st_log_w_neg >= cfg_log_thresh = Reject + | otherwise = Continue +{-# INLINE decide #-} + +-- inspection ----------------------------------------------------------------- + +-- | The current log-wealth, taken as the maximum of the two +-- directional processes. +-- +-- This is the natural \"test statistic\": it is monotone in the +-- evidence against @H_0@ accumulated so far, and the test rejects +-- exactly when it crosses @log(2 \/ alpha)@. +-- +-- >>> log_wealth s0 +-- 0.0 +log_wealth :: State -> Double +log_wealth State{..} = max st_log_w_pos st_log_w_neg +{-# INLINE log_wealth #-} + +-- | The number of samples consumed so far. +-- +-- >>> samples s0 +-- 0 +samples :: State -> Int +samples = st_n +{-# INLINE samples #-} diff --git a/lib/Numeric/Eproc/Mean.hs b/lib/Numeric/Eproc/Mean.hs @@ -1,310 +0,0 @@ -{-# OPTIONS_HADDOCK prune #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE MagicHash #-} -{-# LANGUAGE RecordWildCards #-} - --- | --- Module: Numeric.Eproc.Mean --- Copyright: (c) 2026 Jared Tobin --- License: MIT --- Maintainer: Jared Tobin <jared@ppad.tech> --- --- Two-sided bounded-mean anytime-valid test. --- --- For samples @x_t@ in @[lo, hi]@, tests @H_0: E[x] = m@ against --- @H_1: E[x] /= m@. --- --- Internally two one-sided e-processes are run in parallel: a --- /positive-direction/ process betting against the alternative --- @E[x] > m@ (using centred observations @z = x - m@), and a --- /negative-direction/ process betting against @E[x] < m@ (using --- @-z@). Each maintains its own log-wealth and bettor state. The --- test rejects when either side's wealth crosses @2 \/ alpha@; the --- factor of 2 is the Bonferroni adjustment for the two-sided union. --- --- The test is /anytime-valid/: under @H_0@ the wealth process is a --- nonnegative supermartingale, so by Ville's inequality the --- probability of ever crossing the threshold is at most @alpha@, --- regardless of when the user decides to stop streaming samples. - -module Numeric.Eproc.Mean ( - -- * Test configuration and state - Config - , State - , Verdict(..) - - -- * Construction - , config - , initial - - -- * Streaming - , update - , decide - - -- * Inspection - , log_wealth - , samples - ) where - -import GHC.Exts (Double(D#)) -import Numeric.Eproc.Bettor - --- types ---------------------------------------------------------------------- - --- | Test outcome at the current sample count. --- --- 'Reject' means the wealth process has crossed the Bonferroni --- threshold, so @H_0@ is rejected at level @alpha@. 'Continue' --- means there is not yet enough evidence; collect more samples (or --- stop and report no rejection -- the type-I error guarantee holds --- for /any/ stopping rule). -data Verdict = - Reject - | Continue - deriving (Eq, Show) - --- per-direction bettor state. one constructor per 'Bettor' alternative; --- the constructor used in a given 'State' matches the 'Bettor' chosen --- in the enclosing 'Config'. -data BetState = - SFixed - | SAgrapa - {-# UNPACK #-} !Double -- sum of z (centred observation) - {-# UNPACK #-} !Double -- sum of z^2 (for online variance) - {-# UNPACK #-} !Int -- count - | SOns - {-# UNPACK #-} !Double -- current bet lambda - {-# UNPACK #-} !Double -- running sum of per-step squared gradients - --- | Bounded-mean test configuration. Build with 'config'. --- --- Carries the bettor strategy, the null mean, the significance --- level, the precomputed Bonferroni-adjusted log-wealth threshold, --- and the per-direction safe-bet ceilings (see 'config' for how --- the latter are derived from the sample bounds). -data Config = Config { - cfg_bettor :: !Bettor -- ^ bettor strategy - , cfg_lam_max_pos :: {-# UNPACK #-} !Double -- ^ pos-direction safe-bet ceiling - , cfg_lam_max_neg :: {-# UNPACK #-} !Double -- ^ neg-direction safe-bet ceiling - , cfg_null_mean :: {-# UNPACK #-} !Double -- ^ null mean @m@ - , cfg_alpha :: {-# UNPACK #-} !Double -- ^ significance level @alpha@ - , cfg_log_thresh :: {-# UNPACK #-} !Double -- ^ rejection threshold @log(2 \/ alpha)@ - } - --- | Streaming test state. Construct with 'initial' and fold --- observations through 'update'. --- --- The two log-wealth fields track the running log-wealth of the --- positive- and negative-direction e-processes separately; --- 'decide' compares each to the threshold and 'log_wealth' returns --- the larger of the two. The per-direction bettor states carry --- whatever the chosen 'Bettor' needs (running sums, current bet, --- etc.). -data State = State { - st_n :: {-# UNPACK #-} !Int -- ^ sample count - , st_log_w_pos :: {-# UNPACK #-} !Double -- ^ log-wealth, pos-direction process - , st_log_w_neg :: {-# UNPACK #-} !Double -- ^ log-wealth, neg-direction process - , st_bet_pos :: !BetState -- ^ bettor state, pos-direction - , st_bet_neg :: !BetState -- ^ bettor state, neg-direction - } - --- internal ------------------------------------------------------------------- - --- floor for the wealth factor before taking a log; keeps the running --- log-wealth finite when a step pushes the factor to (or below) zero. --- NB. written via MagicHash because the fractional literal '1.0e-300' --- compiles as 'fromRational (1.0e-300 :: Rational)', and GHC does --- not constant-fold the conversion -- leaving a per-step --- '$wrationalToDouble' call in the worker. -tiny :: Double -tiny = D# 1.0e-300## -{-# INLINE tiny #-} - --- per-bettor initial state. -init_bet :: Bettor -> BetState -init_bet b = case b of - Fixed _ -> SFixed - Agrapa -> SAgrapa 0 0 0 - Ons -> SOns 0 1.0e-6 -- small acc seed avoids div-by-zero on first step -{-# INLINE init_bet #-} - --- compute the next bet 'lambda' from the bettor and its current --- state; 'lam_max' is the direction-specific safety bound. for --- Agrapa we form a Kelly-style plug-in from the running sample mean --- and variance; for Ons the bet is just the last lambda chosen by the --- Newton step (updated during 'step_bet'). -bet_lambda :: Bettor -> Double -> BetState -> Double -bet_lambda b !lam_max !s = case b of - Fixed lam -> lam - Agrapa -> case s of - SAgrapa !sm !sm2 !n - | n == 0 -> 0 - | otherwise -> - let !nd = fromIntegral n - !mu = sm / nd - !mu2 = mu * mu - !var = max 0 (sm2 / nd - mu2) - !den = var + mu2 - !raw = if den == 0 then 0 else mu / den - in max 0 (min lam_max raw) - _ -> 0 - Ons -> case s of - SOns !lam _ -> lam - _ -> 0 -{-# INLINE bet_lambda #-} - --- update bettor state with newly observed centred value 'z'. for --- Agrapa this is just accumulating sums; for Ons we take one Newton --- step on the per-step log-wealth loss '-log(1 + lambda * z)', --- accumulating squared gradients for adaptive scaling. -step_bet :: Bettor -> Double -> BetState -> Double -> BetState -step_bet b !lam_max !s !z = case b of - Fixed _ -> SFixed - Agrapa -> case s of - SAgrapa !sm !sm2 !n -> SAgrapa (sm + z) (sm2 + z * z) (n + 1) - _ -> SAgrapa z (z * z) 1 - Ons -> case s of - SOns !lam !acc -> - let !denom = 1 + lam * z - !g = if denom == 0 then 0 else negate z / denom - !acc' = acc + g * g - !lam' = lam - g / acc' - !clp = max 0 (min lam_max lam') - in SOns clp acc' - _ -> SOns 0 1.0e-6 -{-# INLINE step_bet #-} - --- construction --------------------------------------------------------------- - --- | Build a 'Config' for the bounded-mean test. --- --- Each per-direction safe-bet ceiling @lambda_max@ is set so that --- the wealth factor stays nonnegative for every admissible --- observation: --- --- * The positive-direction factor is @1 + lambda_p * (x - m)@. --- Since @x@ can dip to @lo@, @x - m@ can reach @lo - m@ (the --- most negative value), so we need --- @lambda_p <= 1 \/ (m - lo)@. The ceiling stored is half this --- to leave numerical margin -- the WSR safety recommendation. --- --- * The negative-direction factor is @1 - lambda_n * (x - m)@. --- Since @x@ can rise to @hi@, @x - m@ can reach @hi - m@, so we --- need @lambda_n <= 1 \/ (hi - m)@; again the ceiling is set to --- half this. --- --- The log-wealth rejection threshold is precomputed as --- @log(2 \/ alpha)@; the 2 is the Bonferroni union-bound --- adjustment for the two one-sided e-processes. --- --- >>> import qualified Numeric.Eproc.Bettor as B --- >>> let cfg = config 0.5 0.0 1.0 1.0e-3 B.Ons -config - :: Double -- ^ null mean @m@ - -> Double -- ^ sample lower bound @lo@ - -> Double -- ^ sample upper bound @hi@ - -> Double -- ^ significance level @alpha@ - -> Bettor -- ^ bettor strategy - -> Config -config !m !lo !hi !alpha !b = Config { - cfg_bettor = b - , cfg_lam_max_pos = 0.5 / (m - lo) - , cfg_lam_max_neg = 0.5 / (hi - m) - , cfg_null_mean = m - , cfg_alpha = alpha - , cfg_log_thresh = log (2 / alpha) - } -{-# INLINE config #-} - --- | The initial 'State' for a fresh streaming test. --- --- Both directional log-wealths start at @0@ (i.e., wealth @1@) and --- both bettors start in the per-strategy initial state appropriate --- for the 'Bettor' chosen in the 'Config'. --- --- >>> let s0 = initial cfg -initial :: Config -> State -initial Config{..} = - let !s0 = init_bet cfg_bettor - in State { - st_n = 0 - , st_log_w_pos = 0 - , st_log_w_neg = 0 - , st_bet_pos = s0 - , st_bet_neg = s0 - } -{-# INLINE initial #-} - --- streaming ------------------------------------------------------------------ - --- | Fold one observation into the running 'State'. --- --- Computes the centred observation @z = x - m@, queries the two --- directional bettors for their predictable bets, accumulates --- per-direction log-wealth via --- --- @log_w' = log_w + log (1 + lambda * z)@ --- --- (with the symmetric @-lambda@ for the negative direction), and --- then steps the bettor states given the newly observed @z@. The --- per-step wealth factor is floored at a tiny positive value to --- keep the log finite when a marginal bet drives the factor to (or --- below) zero. --- --- >>> let s1 = update cfg s0 0.7 -update :: Config -> State -> Double -> State -update Config{..} State{..} !x = - let !z = x - cfg_null_mean - !lam_p = bet_lambda cfg_bettor cfg_lam_max_pos st_bet_pos - !lam_n = bet_lambda cfg_bettor cfg_lam_max_neg st_bet_neg - !fac_p = 1 + lam_p * z - !fac_n = 1 - lam_n * z - !logw_p = st_log_w_pos + log (max tiny fac_p) - !logw_n = st_log_w_neg + log (max tiny fac_n) - !sp = step_bet cfg_bettor cfg_lam_max_pos st_bet_pos z - !sn = step_bet cfg_bettor cfg_lam_max_neg st_bet_neg (negate z) - in State (st_n + 1) logw_p logw_n sp sn -{-# INLINE update #-} - --- | Compute the current 'Verdict' from the running 'State'. --- --- 'Reject' iff either directional log-wealth has crossed the --- Bonferroni-adjusted threshold @log(2 \/ alpha)@; equivalently, --- the wealth process on either side has exceeded @2 \/ alpha@. --- Under @H_0@, by Ville's inequality, the probability of this ever --- happening is at most @alpha@ -- and crucially this bound holds --- at /every/ sample size simultaneously, so the user is free to --- peek at the verdict as often as they like and stop on the first --- 'Reject'. --- --- >>> decide cfg s0 --- Continue -decide :: Config -> State -> Verdict -decide Config{..} State{..} - | st_log_w_pos >= cfg_log_thresh = Reject - | st_log_w_neg >= cfg_log_thresh = Reject - | otherwise = Continue -{-# INLINE decide #-} - --- inspection ----------------------------------------------------------------- - --- | The current log-wealth, taken as the maximum of the two --- directional processes. --- --- This is the natural \"test statistic\": it is monotone in the --- evidence against @H_0@ accumulated so far, and the test rejects --- exactly when it crosses @log(2 \/ alpha)@. --- --- >>> log_wealth s0 --- 0.0 -log_wealth :: State -> Double -log_wealth State{..} = max st_log_w_pos st_log_w_neg -{-# INLINE log_wealth #-} - --- | The number of samples consumed so far. --- --- >>> samples s0 --- 0 -samples :: State -> Int -samples = st_n -{-# INLINE samples #-} diff --git a/lib/Numeric/Eproc/Paired.hs b/lib/Numeric/Eproc/Paired.hs @@ -16,7 +16,7 @@ -- The reduction is straightforward: under the null, the differences -- @d_t = a_t - b_t@ have mean zero, and differences of @[lo, hi]@ -- values lie in @[lo - hi, hi - lo]@. So the paired test is just --- the bounded-mean test ("Numeric.Eproc.Mean") on @d_t@ with +-- the bounded-mean test ("Numeric.Eproc.Bounded") on @d_t@ with -- null mean @0@ and sample bounds @[lo - hi, hi - lo]@. -- -- Pairing is required: independent two-sample testing without @@ -43,20 +43,20 @@ module Numeric.Eproc.Paired ( , samples ) where -import qualified Numeric.Eproc.Mean as M -import Numeric.Eproc.Mean (Verdict(..)) +import qualified Numeric.Eproc.Bounded as Bounded +import Numeric.Eproc.Bounded (Verdict(..)) import Numeric.Eproc.Bettor (Bettor) -- types ---------------------------------------------------------------------- -- | Paired two-sample test configuration. Build with 'config'. Wraps --- a 'Numeric.Eproc.Mean.Config' for the underlying +-- a 'Numeric.Eproc.Bounded.Config' for the underlying -- difference test. -newtype Config = Config M.Config +newtype Config = Config Bounded.Config -- | Streaming paired two-sample test state. Construct with 'initial' -- and fold paired observations through 'update'. -newtype State = State M.State +newtype State = State Bounded.State -- construction --------------------------------------------------------------- @@ -77,14 +77,14 @@ config -> Config config !lo !hi !alpha b = let !d = hi - lo - in Config (M.config 0 (negate d) d alpha b) + in Config (Bounded.config 0 (negate d) d alpha b) {-# INLINE config #-} -- | The initial 'State' for a fresh streaming test. -- -- >>> let s0 = initial cfg initial :: Config -> State -initial (Config c) = State (M.initial c) +initial (Config c) = State (Bounded.initial c) {-# INLINE initial #-} -- streaming ------------------------------------------------------------------ @@ -97,7 +97,7 @@ initial (Config c) = State (M.initial c) -- >>> let s1 = update cfg s0 (0.3, 0.7) update :: Config -> State -> (Double, Double) -> State update (Config c) (State s) (!a, !b) = - State (M.update c s (a - b)) + State (Bounded.update c s (a - b)) {-# INLINE update #-} -- | Compute the current 'Verdict' from the running 'State'. @@ -109,7 +109,7 @@ update (Config c) (State s) (!a, !b) = -- >>> decide cfg s0 -- Continue decide :: Config -> State -> Verdict -decide (Config c) (State s) = M.decide c s +decide (Config c) (State s) = Bounded.decide c s {-# INLINE decide #-} -- inspection ----------------------------------------------------------------- @@ -120,7 +120,7 @@ decide (Config c) (State s) = M.decide c s -- >>> log_wealth s0 -- 0.0 log_wealth :: State -> Double -log_wealth (State s) = M.log_wealth s +log_wealth (State s) = Bounded.log_wealth s {-# INLINE log_wealth #-} -- | The number of paired observations consumed so far. @@ -128,5 +128,5 @@ log_wealth (State s) = M.log_wealth s -- >>> samples s0 -- 0 samples :: State -> Int -samples (State s) = M.samples s +samples (State s) = Bounded.samples s {-# INLINE samples #-} diff --git a/ppad-eproc.cabal b/ppad-eproc.cabal @@ -35,7 +35,7 @@ library ghc-options: -fllvm -O2 exposed-modules: Numeric.Eproc.Bettor - Numeric.Eproc.Mean + Numeric.Eproc.Bounded Numeric.Eproc.Paired build-depends: base >= 4.9 && < 5 diff --git a/test/Main.hs b/test/Main.hs @@ -5,7 +5,7 @@ module Main where import Data.Bits import Data.Word import qualified Numeric.Eproc.Bettor as B -import qualified Numeric.Eproc.Mean as M +import qualified Numeric.Eproc.Bounded as Bounded import qualified Numeric.Eproc.Paired as P import Test.Tasty import Test.Tasty.HUnit @@ -53,26 +53,26 @@ gen_seq g = let (_, g') = step_gen g in g : gen_seq g' -- run a sequential mean test on a stream of n bernoulli(p) samples, -- with the early-stopping rule built in. returns (verdict, samples -- consumed). -run_mean_bernoulli - :: M.Config +run_bounded_bernoulli + :: Bounded.Config -> Double -- ^ p -> Int -- ^ budget -> Gen - -> (M.Verdict, Int) -run_mean_bernoulli cfg p budget g0 = go 0 g0 (M.initial cfg) + -> (Bounded.Verdict, Int) +run_bounded_bernoulli cfg p budget g0 = go 0 g0 (Bounded.initial cfg) where go !n !g !st - | n >= budget = (M.decide cfg st, n) - | otherwise = case M.decide cfg st of - M.Reject -> (M.Reject, n) - M.Continue -> + | n >= budget = (Bounded.decide cfg st, n) + | otherwise = case Bounded.decide cfg st of + Bounded.Reject -> (Bounded.Reject, n) + Bounded.Continue -> let (x, g') = bernoulli p g - st' = M.update cfg st x + st' = Bounded.update cfg st x in go (n + 1) g' st' -- fraction of trials that rejected. rejection_rate - :: M.Config + :: Bounded.Config -> Double -- ^ true bernoulli p -> Int -- ^ budget per trial -> Int -- ^ number of trials @@ -82,8 +82,8 @@ rejection_rate cfg p budget trials seed = let gens = take trials (gen_seq (mk_gen seed)) rejects = length [ () | g <- gens - , let (v, _) = run_mean_bernoulli cfg p budget g - , v == M.Reject ] + , let (v, _) = run_bounded_bernoulli cfg p budget g + , v == Bounded.Reject ] in fromIntegral rejects / fromIntegral trials run_paired @@ -98,8 +98,8 @@ run_paired cfg pa pb budget g0 = go 0 g0 (P.initial cfg) go !n !g !st | n >= budget = (P.decide cfg st, n) | otherwise = case P.decide cfg st of - M.Reject -> (M.Reject, n) - M.Continue -> + Bounded.Reject -> (Bounded.Reject, n) + Bounded.Continue -> let (a, g1) = bernoulli pa g (b, g2) = bernoulli pb g1 st' = P.update cfg st (a, b) @@ -118,7 +118,7 @@ paired_avg_rate cfg pa pb budget trials seed = rejects = length [ () | g <- gens , let (v, _) = run_paired cfg pa pb budget g - , v == M.Reject ] + , v == Bounded.Reject ] in fromIntegral rejects / fromIntegral trials -- sanity --------------------------------------------------------------------- @@ -127,13 +127,13 @@ paired_avg_rate cfg pa pb budget trials seed = sanity_tests :: TestTree sanity_tests = testGroup "sanity" [ testCase "degenerate input never rejects" $ do - let cfg = M.config 0.5 0.0 1.0 1.0e-6 B.Ons + let cfg = Bounded.config 0.5 0.0 1.0 1.0e-6 B.Ons xs = replicate 5000 0.5 - st = foldl' (M.update cfg) (M.initial cfg) xs - M.decide cfg st @?= M.Continue + st = foldl' (Bounded.update cfg) (Bounded.initial cfg) xs + Bounded.decide cfg st @?= Bounded.Continue , testCase "two-sided thresholds applied symmetrically" $ do - let cfg = M.config 0.5 0.0 1.0 1.0e-6 B.Ons - M.decide cfg (M.initial cfg) @?= M.Continue + let cfg = Bounded.config 0.5 0.0 1.0 1.0e-6 B.Ons + Bounded.decide cfg (Bounded.initial cfg) @?= Bounded.Continue ] -- null calibration ----------------------------------------------------------- @@ -144,14 +144,14 @@ sanity_tests = testGroup "sanity" [ calibration_tests :: TestTree calibration_tests = testGroup "null calibration" [ testCase "ONS, Bernoulli(0.5), m=0.5, alpha=0.05" $ do - let cfg = M.config 0.5 0.0 1.0 0.05 B.Ons + let cfg = Bounded.config 0.5 0.0 1.0 0.05 B.Ons rate = rejection_rate cfg 0.5 2000 200 12345 -- expected rate <= 0.05; allow up to 0.10 slack for sampling -- variability over 200 trials. assertBool ("FPR " ++ show rate ++ " exceeded slack") $ rate <= 0.10 , testCase "aGRAPA, Bernoulli(0.5), m=0.5, alpha=0.05" $ do - let cfg = M.config 0.5 0.0 1.0 0.05 B.Agrapa + let cfg = Bounded.config 0.5 0.0 1.0 0.05 B.Agrapa rate = rejection_rate cfg 0.5 2000 200 67890 assertBool ("FPR " ++ show rate ++ " exceeded slack") $ rate <= 0.10 @@ -163,12 +163,12 @@ calibration_tests = testGroup "null calibration" [ power_tests :: TestTree power_tests = testGroup "power" [ testCase "ONS detects Bernoulli(0.7) vs m=0.5" $ do - let cfg = M.config 0.5 0.0 1.0 1.0e-3 B.Ons + let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Ons rate = rejection_rate cfg 0.7 5000 100 11111 assertBool ("power " ++ show rate ++ " too low") $ rate >= 0.95 , testCase "aGRAPA detects Bernoulli(0.7) vs m=0.5" $ do - let cfg = M.config 0.5 0.0 1.0 1.0e-3 B.Agrapa + let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Agrapa rate = rejection_rate cfg 0.7 5000 100 22222 assertBool ("power " ++ show rate ++ " too low") $ rate >= 0.95 @@ -195,18 +195,18 @@ two_sample_tests = testGroup "two-sample" [ bettor_smoke_tests :: TestTree bettor_smoke_tests = testGroup "bettor smoke" [ testCase "fixed bettor runs without error" $ do - let cfg = M.config 0.5 0.0 1.0 1.0e-3 (B.Fixed 0.5) + let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 (B.Fixed 0.5) xs = take 100 (cycle [0.0, 1.0]) - st = foldl' (M.update cfg) (M.initial cfg) xs - assertBool "samples advanced" (M.samples st == 100) + st = foldl' (Bounded.update cfg) (Bounded.initial cfg) xs + assertBool "samples advanced" (Bounded.samples st == 100) , testCase "ONS bettor runs without error" $ do - let cfg = M.config 0.5 0.0 1.0 1.0e-3 B.Ons + let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Ons xs = take 100 (cycle [0.0, 1.0]) - st = foldl' (M.update cfg) (M.initial cfg) xs - assertBool "samples advanced" (M.samples st == 100) + st = foldl' (Bounded.update cfg) (Bounded.initial cfg) xs + assertBool "samples advanced" (Bounded.samples st == 100) , testCase "aGRAPA bettor runs without error" $ do - let cfg = M.config 0.5 0.0 1.0 1.0e-3 B.Agrapa + let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 B.Agrapa xs = take 100 (cycle [0.0, 1.0]) - st = foldl' (M.update cfg) (M.initial cfg) xs - assertBool "samples advanced" (M.samples st == 100) + st = foldl' (Bounded.update cfg) (Bounded.initial cfg) xs + assertBool "samples advanced" (Bounded.samples st == 100) ]