ppad-eproc (empty) → 0.1.0
raw patch · 10 files changed
+1362/−0 lines, 10 filesdep +basedep +criteriondep +deepseq
Dependencies added: base, criterion, deepseq, ppad-eproc, tasty, tasty-hunit, tasty-quickcheck, weigh
Files
- CHANGELOG +6/−0
- LICENSE +20/−0
- bench/Main.hs +73/−0
- bench/Weight.hs +65/−0
- lib/Numeric/Eproc/Bernoulli.hs +283/−0
- lib/Numeric/Eproc/Bounded.hs +321/−0
- lib/Numeric/Eproc/Common.hs +70/−0
- lib/Numeric/Eproc/Paired.hs +144/−0
- ppad-eproc.cabal +92/−0
- test/Main.hs +288/−0
+ CHANGELOG view
@@ -0,0 +1,6 @@+# Changelog++- 0.1.0 (2026-06-03)+ * Initial release, supporting anytime-valid sequential testing via+ e-processes: bounded-mean, Bernoulli, and paired two-sample tests,+ with fixed-lambda, aGRAPA, and ONS bettors.
+ LICENSE view
@@ -0,0 +1,20 @@+Copyright (c) 2026 Jared Tobin++Permission is hereby granted, free of charge, to any person obtaining+a copy of this software and associated documentation files (the+"Software"), to deal in the Software without restriction, including+without limitation the rights to use, copy, modify, merge, publish,+distribute, sublicense, and/or sell copies of the Software, and to+permit persons to whom the Software is furnished to do so, subject to+the following conditions:++The above copyright notice and this permission notice shall be included+in all copies or substantial portions of the Software.++THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.+IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY+CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,+TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE+SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ bench/Main.hs view
@@ -0,0 +1,73 @@+{-# OPTIONS_GHC -fno-warn-orphans -fno-warn-type-defaults #-}+{-# LANGUAGE BangPatterns #-}++module Main where++import Control.DeepSeq+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 Bounded.State where rnf !_ = ()+instance NFData P.State where rnf !_ = ()+instance NFData Bounded.Verdict where rnf !_ = ()++main :: IO ()+main = defaultMain [+ update+ , decide+ , stream+ , twosample+ ]++update :: Benchmark+update =+ let !cfg_f = Bounded.config 0.5 0.0 1.0 1.0e-3 (Bounded.Fixed 0.5)+ !cfg_a = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Adaptive+ !cfg_o = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton+ !st_f = Bounded.initial cfg_f+ !st_a = Bounded.initial cfg_a+ !st_o = Bounded.initial cfg_o+ !x = 0.7+ in bgroup "Bounded.update (one step)" [+ bench "fixed" $ nf (Bounded.update cfg_f st_f) x+ , bench "adaptive" $ nf (Bounded.update cfg_a st_a) x+ , bench "newton" $ nf (Bounded.update cfg_o st_o) x+ ]++decide :: Benchmark+decide =+ let !cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton+ !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 = Bounded.config 0.5 0.0 1.0 1.0e-3 (Bounded.Fixed 0.5)+ !cfg_a = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Adaptive+ !cfg_o = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton+ 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 "adaptive" $ nf (run_m cfg_a) xs+ , bench "newton" $ nf (run_m cfg_o) xs+ ]++twosample :: Benchmark+twosample =+ let !ps = force (take 1000 (cycle [(0.3, 0.7), (0.7, 0.3)]))+ !cfg_f = P.config 0.0 1.0 1.0e-3 (Bounded.Fixed 0.5)+ !cfg_a = P.config 0.0 1.0 1.0e-3 Bounded.Adaptive+ !cfg_o = P.config 0.0 1.0 1.0e-3 Bounded.Newton+ run_t cfg = foldl' (P.update cfg) (P.initial cfg)+ in bgroup "Paired.update (1000-sample fold)" [+ bench "fixed" $ nf (run_t cfg_f) ps+ , bench "adaptive" $ nf (run_t cfg_a) ps+ , bench "newton" $ nf (run_t cfg_o) ps+ ]
+ bench/Weight.hs view
@@ -0,0 +1,65 @@+{-# OPTIONS_GHC -fno-warn-orphans -fno-warn-type-defaults #-}+{-# LANGUAGE BangPatterns #-}++module Main where++import Control.DeepSeq+import qualified Numeric.Eproc.Bounded as Bounded+import qualified Numeric.Eproc.Paired as P+import Weigh++instance NFData Bounded.State where rnf !_ = ()+instance NFData P.State where rnf !_ = ()+instance NFData Bounded.Verdict where rnf !_ = ()++-- note that 'weigh' doesn't work properly in a repl+main :: IO ()+main = mainWith $ do+ update+ decide+ stream+ twosample++update :: Weigh ()+update =+ let !cfg_f = Bounded.config 0.5 0.0 1.0 1.0e-3 (Bounded.Fixed 0.5)+ !cfg_a = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Adaptive+ !cfg_o = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton+ !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 "adaptive" (Bounded.update cfg_a st_a) 0.7+ func "newton" (Bounded.update cfg_o st_o) 0.7++decide :: Weigh ()+decide =+ let !cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton+ !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 = Bounded.config 0.5 0.0 1.0 1.0e-3 (Bounded.Fixed 0.5)+ !cfg_a = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Adaptive+ !cfg_o = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton+ 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 "adaptive" (run_m cfg_a) xs+ func "newton" (run_m cfg_o) xs++twosample :: Weigh ()+twosample =+ let !ps = force (take 1000 (cycle [(0.3, 0.7), (0.7, 0.3)]))+ !cfg_f = P.config 0.0 1.0 1.0e-3 (Bounded.Fixed 0.5)+ !cfg_a = P.config 0.0 1.0 1.0e-3 Bounded.Adaptive+ !cfg_o = P.config 0.0 1.0 1.0e-3 Bounded.Newton+ run_t cfg = foldl' (P.update cfg) (P.initial cfg)+ in wgroup "Paired.update (1000-sample fold)" $ do+ func "fixed" (run_t cfg_f) ps+ func "adaptive" (run_t cfg_a) ps+ func "newton" (run_t cfg_o) ps
+ lib/Numeric/Eproc/Bernoulli.hs view
@@ -0,0 +1,283 @@+{-# OPTIONS_HADDOCK prune #-}+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE RecordWildCards #-}++-- |+-- Module: Numeric.Eproc.Bernoulli+-- Copyright: (c) 2026 Jared Tobin+-- License: MIT+-- Maintainer: Jared Tobin <jared@ppad.tech>+--+-- One-sided Bernoulli rate anytime-valid test.+--+-- For samples @x_t@ in @{0, 1}@, tests @H_0: E[x] <= p_0@ against+-- @H_1: E[x] > p_0@.+--+-- A single wealth process is run:+--+-- @W_n = prod_{i=1..n} (1 + lambda_i * (x_i - p_0))@+--+-- where each per-step bet @lambda_i@ is chosen predictably (from+-- data observed strictly before step @i@) and clipped to+-- @[0, lambda_max]@ so that the wealth factor stays nonnegative for+-- every admissible observation. Under @H_0@ the wealth process is+-- a nonnegative supermartingale, so by Ville's inequality the+-- probability of @W_n@ ever crossing @1 \/ alpha@ is at most+-- @alpha@, regardless of when the user decides to stop streaming+-- samples.+--+-- Unlike "Numeric.Eproc.Bounded", the alternative here is one-sided,+-- so a single wealth process suffices and no Bonferroni adjustment+-- is needed -- the rejection threshold is @log(1 \/ alpha)@.+--+-- == Example+--+-- Test @H_0: E[x] <= 0.05@ at level @alpha = 1e-3@ against a stream+-- with empirical rate @~0.5@:+--+-- >>> let cfg = config 1.0e-3 0.05 Newton+-- >>> let xs = take 200 (cycle [True, False])+-- >>> decide cfg (foldl' (update cfg) (initial cfg) xs)+-- Reject++module Numeric.Eproc.Bernoulli (+ -- * Test configuration and state+ Config+ , State+ , Verdict(..)++ -- * Bettor strategies+ , Bettor(..)++ -- * Construction+ , config+ , initial++ -- * Streaming+ , update+ , decide++ -- * Inspection+ , log_wealth+ , samples+ ) where++import Numeric.Eproc.Common (Bettor(..), Verdict(..))++-- types ----------------------------------------------------------------------++-- here, the centred observation @z_t@ referenced in+-- "Numeric.Eproc.Common" is @x_t - p_0@; the safe-bet ceiling+-- @lambda_max@ is derived from @p_0@ (see 'config').++-- 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+ | SAdaptive+ {-# UNPACK #-} !Double -- sum of z (centred observation)+ {-# UNPACK #-} !Double -- sum of z^2 (for online variance)+ {-# UNPACK #-} !Int -- count+ | SNewton+ {-# UNPACK #-} !Double -- current bet lambda+ {-# UNPACK #-} !Double -- running sum of per-step squared gradients++-- | Bernoulli rate test configuration. Build with 'config'.+--+-- Carries the bettor strategy, the baseline rate, the significance+-- level, the precomputed log-wealth rejection threshold, and the+-- safe-bet ceiling derived from @p_0@.+data Config = Config {+ -- ^ bettor strategy+ cfg_bettor :: !Bettor+ -- ^ safe-bet ceiling+ , cfg_lam_max :: {-# UNPACK #-} !Double+ -- ^ baseline rate @p_0@+ , cfg_p0 :: {-# UNPACK #-} !Double+ -- ^ significance level @alpha@+ , cfg_alpha :: {-# UNPACK #-} !Double+ -- ^ rejection threshold @log(1 \/ alpha)@+ , cfg_log_thresh :: {-# UNPACK #-} !Double+ }++-- | Streaming test state. Construct with 'initial' and fold+-- observations through 'update'.+--+-- Carries the sample count, running log-wealth, and whatever+-- per-step state the chosen 'Bettor' needs.+data State = State {+ st_n :: {-# UNPACK #-} !Int -- ^ sample count+ , st_log_w :: {-# UNPACK #-} !Double -- ^ running log-wealth+ , st_bet :: !BetState -- ^ bettor state+ }++-- internal -------------------------------------------------------------------++-- per-bettor initial state.+init_bet :: Bettor -> BetState+init_bet b = case b of+ Fixed _ -> SFixed+ Adaptive -> SAdaptive 0 0 0+ Newton -> SNewton 0 1.0e-6 -- small acc seed avoids div-by-zero+{-# INLINE init_bet #-}++-- compute the next bet 'lambda' from the bettor and its current+-- state. for Adaptive we form a Kelly-style plug-in from the running+-- sample mean and variance; for Newton 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+ Adaptive -> case s of+ SAdaptive !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+ Newton -> case s of+ SNewton !lam _ -> lam+ _ -> 0+{-# INLINE bet_lambda #-}++-- update bettor state with newly observed centred value 'z'. for+-- Adaptive this is just accumulating sums; for Newton 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+ Adaptive -> case s of+ SAdaptive !sm !sm2 !n -> SAdaptive (sm + z) (sm2 + z * z) (n + 1)+ _ -> SAdaptive z (z * z) 1+ Newton -> case s of+ SNewton !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 SNewton clp acc'+ _ -> SNewton 0 1.0e-6+{-# INLINE step_bet #-}++-- construction ---------------------------------------------------------------++-- | Build a 'Config' for the Bernoulli rate test.+--+-- The safe-bet ceiling @lambda_max@ is set so that the wealth+-- factor @1 + lambda * (x - p_0)@ stays nonnegative for both+-- @x = 0@ and @x = 1@. The binding constraint is @x = 0@, which+-- requires @lambda <= 1 \/ p_0@; the ceiling stored is half this+-- to leave numerical margin -- the WSR safety recommendation.+--+-- @p_0@ must lie strictly in @(0, 1)@ and @alpha@ strictly in+-- @(0, 1)@. The degenerate case @p_0 = 0@ would make @lambda_max@+-- infinite (any divergence would reject immediately and the test+-- becomes uninteresting); the caller is expected to pass a small+-- positive baseline.+--+-- >>> let cfg = config 1.0e-3 0.05 Newton+config+ :: Double -- ^ significance level @alpha@, in @(0, 1)@+ -> Double -- ^ baseline rate @p_0@, in @(0, 1)@+ -> Bettor -- ^ bettor strategy+ -> Config+config !alpha !p0 !b = Config {+ cfg_bettor = b+ , cfg_lam_max = 0.5 / p0+ , cfg_p0 = p0+ , cfg_alpha = alpha+ , cfg_log_thresh = log (1 / alpha)+ }+{-# INLINE config #-}++-- | The initial 'State' for a fresh streaming test.+--+-- Log-wealth starts at @0@ (i.e., wealth @1@) and the bettor+-- starts in the per-strategy initial state appropriate for the+-- 'Bettor' chosen in the 'Config'.+--+-- >>> let s0 = initial cfg+initial :: Config -> State+initial Config{..} = State {+ st_n = 0+ , st_log_w = 0+ , st_bet = init_bet cfg_bettor+ }+{-# INLINE initial #-}++-- streaming ------------------------------------------------------------------++-- | Fold one observation into the running 'State'.+--+-- @True@ means @x_t = 1@ (the event of interest occurred -- e.g.,+-- two readings diverged); @False@ means @x_t = 0@ (they matched).+-- The caller decides what \"matched\" means at the application+-- level.+--+-- Computes the centred observation @z = x - p_0@, queries the+-- bettor for its predictable bet, accumulates log-wealth via+--+-- @log_w' = log_w + log (1 + lambda * z)@+--+-- and then steps the bettor state given the newly observed @z@.+--+-- >>> let s1 = update cfg s0 True+update :: Config -> State -> Bool -> State+update Config{..} State{..} !x =+ let !xd = if x then 1 else 0+ !z = xd - cfg_p0+ !lam = bet_lambda cfg_bettor cfg_lam_max st_bet+ !fac = 1 + lam * z+ !logw' = st_log_w + log fac+ !s' = step_bet cfg_bettor cfg_lam_max st_bet z+ in State (st_n + 1) logw' s'+{-# INLINE update #-}++-- | Compute the current 'Verdict' from the running 'State'.+--+-- 'Reject' iff log-wealth has crossed the threshold+-- @log(1 \/ alpha)@; equivalently, wealth has exceeded+-- @1 \/ 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 >= cfg_log_thresh = Reject+ | otherwise = Continue+{-# INLINE decide #-}++-- inspection -----------------------------------------------------------------++-- | The current log-wealth.+--+-- This is the natural \"test statistic\": it is monotone (in+-- expectation under @H_1@) in the evidence against @H_0@+-- accumulated so far, and the test rejects exactly when it crosses+-- @log(1 \/ alpha)@.+--+-- >>> log_wealth s0+-- 0.0+log_wealth :: State -> Double+log_wealth = st_log_w+{-# INLINE log_wealth #-}++-- | The number of samples consumed so far.+--+-- >>> samples s0+-- 0+samples :: State -> Int+samples = st_n+{-# INLINE samples #-}
+ lib/Numeric/Eproc/Bounded.hs view
@@ -0,0 +1,321 @@+{-# 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.+--+-- == Example+--+-- Test @H_0: E[x] = 0.5@ for @x@ in @[0, 1]@ at level @alpha = 1e-3@+-- against a stream with empirical mean @0.8@:+--+-- >>> let cfg = config 0.5 0.0 1.0 1.0e-3 Newton+-- >>> let xs = concat (replicate 30 [1, 1, 0, 1, 1, 0, 1, 1, 1, 1])+-- >>> decide cfg (foldl' (update cfg) (initial cfg) xs)+-- Reject++module Numeric.Eproc.Bounded (+ -- * Test configuration and state+ Config+ , State+ , Verdict(..)++ -- * Bettor strategies+ , Bettor(..)++ -- * Construction+ , config+ , initial++ -- * Streaming+ , update+ , decide++ -- * Inspection+ , log_wealth+ , samples+ ) where++import GHC.Exts (Double(D#))+import Numeric.Eproc.Common (Bettor(..), Verdict(..))++-- types ----------------------------------------------------------------------++-- here, the centred observation @z_t@ referenced in+-- "Numeric.Eproc.Common" is @x_t - m@; the per-direction safe-bet+-- ceilings @lambda_max@ are derived from the sample bounds (see+-- 'config').++-- 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+ | SAdaptive+ {-# UNPACK #-} !Double -- sum of z (centred observation)+ {-# UNPACK #-} !Double -- sum of z^2 (for online variance)+ {-# UNPACK #-} !Int -- count+ | SNewton+ {-# 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 {+ -- ^ bettor strategy+ cfg_bettor :: !Bettor+ -- ^ positive-direction safe-bet ceiling+ , cfg_lam_max_pos :: {-# UNPACK #-} !Double+ -- ^ negative-direction safe-bet ceiling+ , cfg_lam_max_neg :: {-# UNPACK #-} !Double+ -- ^ null mean @m@+ , cfg_null_mean :: {-# UNPACK #-} !Double+ -- ^ significance level @alpha@+ , cfg_alpha :: {-# UNPACK #-} !Double+ -- ^ rejection threshold @log(2 \/ alpha)@+ , cfg_log_thresh :: {-# UNPACK #-} !Double+ }++-- | 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-dir process+ , st_log_w_neg :: {-# UNPACK #-} !Double -- ^ log-wealth, neg-dir 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+ Adaptive -> SAdaptive 0 0 0+ Newton -> SNewton 0 1.0e-6 -- small acc seed avoids div-by-zero+{-# INLINE init_bet #-}++-- compute the next bet 'lambda' from the bettor and its current+-- state; 'lam_max' is the direction-specific safety bound. for+-- Adaptive we form a Kelly-style plug-in from the running sample+-- mean and variance; for Newton 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+ Adaptive -> case s of+ SAdaptive !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+ Newton -> case s of+ SNewton !lam _ -> lam+ _ -> 0+{-# INLINE bet_lambda #-}++-- update bettor state with newly observed centred value 'z'. for+-- Adaptive this is just accumulating sums; for Newton 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+ Adaptive -> case s of+ SAdaptive !sm !sm2 !n -> SAdaptive (sm + z) (sm2 + z * z) (n + 1)+ _ -> SAdaptive z (z * z) 1+ Newton -> case s of+ SNewton !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 SNewton clp acc'+ _ -> SNewton 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.+--+-- >>> let cfg = config 0.5 0.0 1.0 1.0e-3 Newton+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 #-}
+ lib/Numeric/Eproc/Common.hs view
@@ -0,0 +1,70 @@+{-# OPTIONS_HADDOCK prune #-}++-- |+-- Module: Numeric.Eproc.Common+-- Copyright: (c) 2026 Jared Tobin+-- License: MIT+-- Maintainer: Jared Tobin <jared@ppad.tech>+--+-- Shared vocabulary for the eproc tests: the predictable bettor+-- strategies and the test verdict type. Re-exported from each test+-- module ("Numeric.Eproc.Bounded", "Numeric.Eproc.Paired",+-- "Numeric.Eproc.Bernoulli"); import this module directly only if+-- you need the types without picking a particular test.++module Numeric.Eproc.Common (+ Bettor(..)+ , Verdict(..)+ ) where++-- | A predictable bettor.+--+-- A bettor describes how, given the history of centred+-- observations @z_t@ (each test module specifies its own centring;+-- see the per-module documentation), the next predictable bet+-- @lambda_t@ is chosen. Predictability -- that is, @lambda_t@+-- depends only on data observed strictly before step @t@ -- is+-- what makes the resulting wealth process a nonnegative+-- supermartingale under @H_0@.+--+-- For 'Adaptive' and 'Newton', a safe-bet ceiling @lambda_max@+-- derived from the test's admissible-observation range is enforced+-- by clipping @lambda@ to @[0, lambda_max]@, so the wealth factor+-- stays nonnegative.+--+-- * 'Fixed' always bets the supplied constant @lambda@. The wager+-- does not respond to observed data; this strategy is useful+-- only as a baseline.+--+-- * 'Adaptive' is the aGRAPA (approximate growth-rate adaptive+-- predictable plug-in) bettor of Waudby-Smith & Ramdas (2024).+-- It tracks the empirical mean @mu@ and variance @sigma^2@ of+-- centred observations and bets the Kelly-optimal plug-in+-- @lambda* = mu \/ (sigma^2 + mu^2)@ clipped to+-- @[0, lambda_max]@. Fast to compute and competitive in+-- practice.+--+-- * 'Newton' is the online Newton step (ONS) bettor. The per-step+-- log-wealth loss @-log(1 + lambda * z)@ is convex in @lambda@;+-- ONS performs one Newton step per observation, accumulating+-- squared gradients to scale the update. Achieves logarithmic+-- regret against the best constant bet in hindsight and is in+-- practice the strongest of the three bettors under most signal+-- regimes.+data Bettor =+ Fixed {-# UNPACK #-} !Double+ | Adaptive+ | Newton+ deriving (Eq, Show)++-- | Test outcome at the current sample count.+--+-- 'Reject' means the wealth process has crossed the rejection+-- 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)
+ lib/Numeric/Eproc/Paired.hs view
@@ -0,0 +1,144 @@+{-# OPTIONS_HADDOCK prune #-}+{-# LANGUAGE BangPatterns #-}++-- |+-- Module: Numeric.Eproc.Paired+-- Copyright: (c) 2026 Jared Tobin+-- License: MIT+-- Maintainer: Jared Tobin <jared@ppad.tech>+--+-- Paired two-sample anytime-valid mean-equality test.+--+-- For paired observations @(a_t, b_t)@ where both samples lie in+-- @[lo, hi]@, tests @H_0: E[a] = E[b]@ against+-- @H_1: E[a] /= E[b]@.+--+-- 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.Bounded") on @d_t@ with+-- null mean @0@ and sample bounds @[lo - hi, hi - lo]@.+--+-- Pairing is required: independent two-sample testing without+-- alignment would need to bet against a richer alternative (the+-- joint distribution rather than the marginal difference) and is+-- beyond the scope of this module.+--+-- == Example+--+-- Test @H_0: E[a] = E[b]@ for samples in @[0, 1]@ at level+-- @alpha = 1e-3@ against a stream of paired observations where @a@+-- runs systematically higher than @b@:+--+-- >>> let cfg = config 0.0 1.0 1.0e-3 Newton+-- >>> let ps = take 1000 (cycle [(1, 0), (1, 0), (0, 0), (1, 1)])+-- >>> decide cfg (foldl' (update cfg) (initial cfg) ps)+-- Reject++module Numeric.Eproc.Paired (+ -- * Test configuration and state+ Config+ , State+ , Verdict(..)++ -- * Bettor strategies+ , Bettor(..)++ -- * Construction+ , config+ , initial++ -- * Streaming+ , update+ , decide++ -- * Inspection+ , log_wealth+ , samples+ ) where++import qualified Numeric.Eproc.Bounded as Bounded+import Numeric.Eproc.Common (Bettor(..), Verdict(..))++-- types ----------------------------------------------------------------------++-- | Paired two-sample test configuration. Build with 'config'. Wraps+-- a 'Numeric.Eproc.Bounded.Config' for the underlying+-- difference test.+newtype Config = Config Bounded.Config++-- | Streaming paired two-sample test state. Construct with 'initial'+-- and fold paired observations through 'update'.+newtype State = State Bounded.State++-- construction ---------------------------------------------------------------++-- | Build a 'Config' for the paired two-sample test.+--+-- Bounds @lo@ and @hi@ are the (shared) bounds on the individual+-- @a@ and @b@ samples; the underlying mean test is then configured+-- on the differences, which lie in @[lo - hi, hi - lo]@ with null+-- mean @0@.+--+-- >>> let cfg = config 0.0 1.0 1.0e-3 Newton+config+ :: Double -- ^ sample lower bound @lo@+ -> Double -- ^ sample upper bound @hi@+ -> Double -- ^ significance level @alpha@+ -> Bettor -- ^ bettor strategy+ -> Config+config !lo !hi !alpha b =+ let !d = hi - lo+ 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 (Bounded.initial c)+{-# INLINE initial #-}++-- streaming ------------------------------------------------------------------++-- | Fold one paired observation @(a, b)@ into the running 'State'.+--+-- Equivalent to feeding the difference @a - b@ into the underlying+-- bounded-mean test.+--+-- >>> let s1 = update cfg s0 (0.3, 0.7)+update :: Config -> State -> (Double, Double) -> State+update (Config c) (State s) (!a, !b) =+ State (Bounded.update c s (a - b))+{-# INLINE update #-}++-- | Compute the current 'Verdict' from the running 'State'.+--+-- 'Reject' iff either directional log-wealth of the underlying+-- bounded-mean test on the differences has crossed+-- @log(2 \/ alpha)@.+--+-- >>> decide cfg s0+-- Continue+decide :: Config -> State -> Verdict+decide (Config c) (State s) = Bounded.decide c s+{-# INLINE decide #-}++-- inspection -----------------------------------------------------------------++-- | The current log-wealth of the underlying bounded-mean test on+-- the differences.+--+-- >>> log_wealth s0+-- 0.0+log_wealth :: State -> Double+log_wealth (State s) = Bounded.log_wealth s+{-# INLINE log_wealth #-}++-- | The number of paired observations consumed so far.+--+-- >>> samples s0+-- 0+samples :: State -> Int+samples (State s) = Bounded.samples s+{-# INLINE samples #-}
+ ppad-eproc.cabal view
@@ -0,0 +1,92 @@+cabal-version: 3.0+name: ppad-eproc+version: 0.1.0+synopsis: Anytime-valid sequential testing via e-processes.+license: MIT+license-file: LICENSE+author: Jared Tobin+maintainer: jared@ppad.tech+category: Statistics+build-type: Simple+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-sided Bernoulli rate tests with fixed, adaptive (aGRAPA), and+ online Newton bettors.++flag llvm+ description: Use GHC's LLVM backend.+ default: False+ manual: True++source-repository head+ type: git+ location: git.ppad.tech/eproc.git++library+ default-language: Haskell2010+ hs-source-dirs: lib+ ghc-options:+ -Wall+ if flag(llvm)+ ghc-options: -fllvm -O2+ exposed-modules:+ Numeric.Eproc.Bernoulli+ Numeric.Eproc.Bounded+ Numeric.Eproc.Common+ Numeric.Eproc.Paired+ build-depends:+ base >= 4.9 && < 5++test-suite eproc-tests+ type: exitcode-stdio-1.0+ default-language: Haskell2010+ hs-source-dirs: test+ main-is: Main.hs++ ghc-options:+ -rtsopts -Wall -O2++ build-depends:+ base+ , ppad-eproc+ , tasty+ , tasty-hunit+ , tasty-quickcheck++benchmark eproc-bench+ type: exitcode-stdio-1.0+ default-language: Haskell2010+ hs-source-dirs: bench+ main-is: Main.hs++ ghc-options:+ -rtsopts -O2 -Wall -fno-warn-orphans+ if flag(llvm)+ ghc-options: -fllvm++ build-depends:+ base+ , criterion+ , deepseq+ , ppad-eproc++benchmark eproc-weigh+ type: exitcode-stdio-1.0+ default-language: Haskell2010+ hs-source-dirs: bench+ main-is: Weight.hs++ ghc-options:+ -rtsopts -O2 -Wall -fno-warn-orphans+ if flag(llvm)+ ghc-options: -fllvm++ build-depends:+ base+ , deepseq+ , ppad-eproc+ , weigh
+ test/Main.hs view
@@ -0,0 +1,288 @@+{-# LANGUAGE BangPatterns #-}++module Main where++import Data.Bits+import Data.Word+import qualified Numeric.Eproc.Bernoulli as Bern+import qualified Numeric.Eproc.Bounded as Bounded+import qualified Numeric.Eproc.Paired as P+import Test.Tasty+import Test.Tasty.HUnit++main :: IO ()+main = defaultMain $ testGroup "ppad-eproc" [+ sanity_tests+ , calibration_tests+ , power_tests+ , two_sample_tests+ , bernoulli_tests+ , bettor_smoke_tests+ ]++-- prng -----------------------------------------------------------------------++-- inline PCG-style PRNG, no external deps.++newtype Gen = Gen Word64++mk_gen :: Word64 -> Gen+mk_gen = Gen++step_gen :: Gen -> (Word64, Gen)+step_gen (Gen s) =+ let !s' = s * 6364136223846793005 + 1442695040888963407+ in (s', Gen s')++next_double :: Gen -> (Double, Gen)+next_double g =+ let (w, g') = step_gen g+ !x = fromIntegral (w `shiftR` 11 .&. 0x1FFFFFFFFFFFFF) /+ 9007199254740992+ in (x, g')++bernoulli :: Double -> Gen -> (Double, Gen)+bernoulli !p g =+ let (u, g') = next_double g+ in (if u < p then 1.0 else 0.0, g')++-- per-trial independent seeds via a splitmix-style finalizer.+-- previously this just stepped the prng once per trial, which made+-- consecutive trials share all but one observation -- fine under a+-- symmetric H_0 (rare streaks cancel), catastrophic under a skewed+-- one (rare streaks dominate all overlapping trials).+gen_seq :: Gen -> [Gen]+gen_seq (Gen s0) =+ [Gen (mix64 (s0 + fromIntegral i)) | i <- [(0 :: Word64) ..]]+ where+ mix64 x =+ let !y = (x `xor` (x `shiftR` 30)) * 0xbf58476d1ce4e5b9+ !z = (y `xor` (y `shiftR` 27)) * 0x94d049bb133111eb+ in z `xor` (z `shiftR` 31)++-- harness --------------------------------------------------------------------++-- 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_bounded_bernoulli+ :: Bounded.Config+ -> Double -- ^ p+ -> Int -- ^ budget+ -> Gen+ -> (Bounded.Verdict, Int)+run_bounded_bernoulli cfg p budget g0 = go 0 g0 (Bounded.initial cfg)+ where+ go !n !g !st+ | 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' = Bounded.update cfg st x+ in go (n + 1) g' st'++-- fraction of trials that rejected.+rejection_rate+ :: Bounded.Config+ -> Double -- ^ true bernoulli p+ -> Int -- ^ budget per trial+ -> Int -- ^ number of trials+ -> Word64 -- ^ seed+ -> Double+rejection_rate cfg p budget trials seed =+ let gens = take trials (gen_seq (mk_gen seed))+ rejects = length+ [ () | g <- gens+ , let (v, _) = run_bounded_bernoulli cfg p budget g+ , v == Bounded.Reject ]+ in fromIntegral rejects / fromIntegral trials++run_paired+ :: P.Config+ -> Double+ -> Double -- ^ p for A and B+ -> Int+ -> Gen+ -> (P.Verdict, Int)+run_paired cfg pa pb budget g0 = go 0 g0 (P.initial cfg)+ where+ go !n !g !st+ | n >= budget = (P.decide cfg st, n)+ | otherwise = case P.decide cfg st of+ 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)+ in go (n + 1) g2 st'++paired_avg_rate+ :: P.Config+ -> Double+ -> Double+ -> Int+ -> Int+ -> Word64+ -> Double+paired_avg_rate cfg pa pb budget trials seed =+ let gens = take trials (gen_seq (mk_gen seed))+ rejects = length+ [ () | g <- gens+ , let (v, _) = run_paired cfg pa pb budget g+ , v == Bounded.Reject ]+ in fromIntegral rejects / fromIntegral trials++-- sanity ---------------------------------------------------------------------++-- with all-zero deviations from the null mean, no rejection.+sanity_tests :: TestTree+sanity_tests = testGroup "sanity" [+ testCase "degenerate input never rejects" $ do+ let cfg = Bounded.config 0.5 0.0 1.0 1.0e-6 Bounded.Newton+ xs = replicate 5000 0.5+ st = foldl' (Bounded.update cfg) (Bounded.initial cfg) xs+ Bounded.decide cfg st @?= Bounded.Continue+ , testCase "two-sided thresholds applied symmetrically" $ do+ let cfg = Bounded.config 0.5 0.0 1.0 1.0e-6 Bounded.Newton+ Bounded.decide cfg (Bounded.initial cfg) @?= Bounded.Continue+ ]++-- null calibration -----------------------------------------------------------++-- under H_0, with optional stopping, the empirical rejection rate should be+-- bounded by alpha. ville's inequality is typically conservative on bernoulli,+-- so the slack is small.+calibration_tests :: TestTree+calibration_tests = testGroup "null calibration" [+ testCase "Newton, Bernoulli(0.5), m=0.5, alpha=0.05" $ do+ let cfg = Bounded.config 0.5 0.0 1.0 0.05 Bounded.Newton+ 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 "Adaptive, Bernoulli(0.5), m=0.5, alpha=0.05" $ do+ let cfg = Bounded.config 0.5 0.0 1.0 0.05 Bounded.Adaptive+ rate = rejection_rate cfg 0.5 2000 200 67890+ assertBool ("FPR " ++ show rate ++ " exceeded slack") $+ rate <= 0.10+ ]++-- power ----------------------------------------------------------------------++-- under a clear shift, all (or nearly all) trials reject within budget.+power_tests :: TestTree+power_tests = testGroup "power" [+ testCase "Newton detects Bernoulli(0.7) vs m=0.5" $ do+ let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton+ rate = rejection_rate cfg 0.7 5000 100 11111+ assertBool ("power " ++ show rate ++ " too low") $+ rate >= 0.95+ , testCase "Adaptive detects Bernoulli(0.7) vs m=0.5" $ do+ let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Adaptive+ rate = rejection_rate cfg 0.7 5000 100 22222+ assertBool ("power " ++ show rate ++ " too low") $+ rate >= 0.95+ ]++-- two-sample paired test -----------------------------------------------------++two_sample_tests :: TestTree+two_sample_tests = testGroup "two-sample" [+ testCase "identical distributions don't reject" $ do+ let cfg = P.config 0.0 1.0 1.0e-3 Bounded.Newton+ rate = paired_avg_rate cfg 0.5 0.5 2000 100 33333+ assertBool ("FPR " ++ show rate) $ rate <= 0.05+ , testCase "different distributions reject" $ do+ let cfg = P.config 0.0 1.0 1.0e-3 Bounded.Newton+ rate = paired_avg_rate cfg 0.3 0.7 5000 100 44444+ assertBool ("power " ++ show rate) $ rate >= 0.95+ ]++-- bernoulli (one-sided rate) -------------------------------------------------++run_bernoulli+ :: Bern.Config+ -> Double -- ^ true rate p+ -> Int -- ^ budget+ -> Gen+ -> (Bern.Verdict, Int)+run_bernoulli cfg p budget g0 = go 0 g0 (Bern.initial cfg)+ where+ go !n !g !st+ | n >= budget = (Bern.decide cfg st, n)+ | otherwise = case Bern.decide cfg st of+ Bern.Reject -> (Bern.Reject, n)+ Bern.Continue ->+ let (u, g') = next_double g+ !x = u < p+ st' = Bern.update cfg st x+ in go (n + 1) g' st'++bernoulli_rate+ :: Bern.Config+ -> Double -- ^ true rate p+ -> Int -- ^ budget per trial+ -> Int -- ^ number of trials+ -> Word64 -- ^ seed+ -> Double+bernoulli_rate cfg p budget trials seed =+ let gens = take trials (gen_seq (mk_gen seed))+ rejects = length+ [ () | g <- gens+ , let (v, _) = run_bernoulli cfg p budget g+ , v == Bern.Reject ]+ in fromIntegral rejects / fromIntegral trials++bernoulli_tests :: TestTree+bernoulli_tests = testGroup "bernoulli" [+ testCase "all-zero stream never rejects" $ do+ let cfg = Bern.config 1.0e-6 0.05 Bern.Newton+ xs = replicate 5000 False+ st = foldl' (Bern.update cfg) (Bern.initial cfg) xs+ Bern.decide cfg st @?= Bern.Continue+ , testCase "Newton FPR under H_0 (p = p_0 = 0.05)" $ do+ let cfg = Bern.config 0.05 0.05 Bern.Newton+ rate = bernoulli_rate cfg 0.05 2000 200 55555+ assertBool ("FPR " ++ show rate ++ " exceeded slack") $+ rate <= 0.10+ , testCase "Adaptive FPR under H_0 (p = p_0 = 0.05)" $ do+ let cfg = Bern.config 0.05 0.05 Bern.Adaptive+ rate = bernoulli_rate cfg 0.05 2000 200 66666+ assertBool ("FPR " ++ show rate ++ " exceeded slack") $+ rate <= 0.10+ , testCase "Newton detects p = 0.3 vs p_0 = 0.05" $ do+ let cfg = Bern.config 1.0e-3 0.05 Bern.Newton+ rate = bernoulli_rate cfg 0.3 5000 100 77777+ assertBool ("power " ++ show rate ++ " too low") $+ rate >= 0.95+ , testCase "Adaptive detects p = 0.3 vs p_0 = 0.05" $ do+ let cfg = Bern.config 1.0e-3 0.05 Bern.Adaptive+ rate = bernoulli_rate cfg 0.3 5000 100 88888+ assertBool ("power " ++ show rate ++ " too low") $+ rate >= 0.95+ ]++-- bettor smoke tests ---------------------------------------------------------++-- each bettor produces a well-defined state and decision when run on a small+-- deterministic stream.+bettor_smoke_tests :: TestTree+bettor_smoke_tests = testGroup "bettor smoke" [+ testCase "fixed bettor runs without error" $ do+ let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 (Bounded.Fixed 0.5)+ xs = take 100 (cycle [0.0, 1.0])+ st = foldl' (Bounded.update cfg) (Bounded.initial cfg) xs+ assertBool "samples advanced" (Bounded.samples st == 100)+ , testCase "Newton bettor runs without error" $ do+ let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton+ xs = take 100 (cycle [0.0, 1.0])+ st = foldl' (Bounded.update cfg) (Bounded.initial cfg) xs+ assertBool "samples advanced" (Bounded.samples st == 100)+ , testCase "Adaptive bettor runs without error" $ do+ let cfg = Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Adaptive+ xs = take 100 (cycle [0.0, 1.0])+ st = foldl' (Bounded.update cfg) (Bounded.initial cfg) xs+ assertBool "samples advanced" (Bounded.samples st == 100)+ ]