diff --git a/CHANGELOG b/CHANGELOG new file mode 100644 --- /dev/null +++ b/CHANGELOG @@ -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. diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/bench/Main.hs b/bench/Main.hs new file mode 100644 --- /dev/null +++ b/bench/Main.hs @@ -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 + ] diff --git a/bench/Weight.hs b/bench/Weight.hs new file mode 100644 --- /dev/null +++ b/bench/Weight.hs @@ -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 diff --git a/lib/Numeric/Eproc/Bernoulli.hs b/lib/Numeric/Eproc/Bernoulli.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Eproc/Bernoulli.hs @@ -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 +-- +-- 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 #-} diff --git a/lib/Numeric/Eproc/Bounded.hs b/lib/Numeric/Eproc/Bounded.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Eproc/Bounded.hs @@ -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 +-- +-- 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 #-} diff --git a/lib/Numeric/Eproc/Common.hs b/lib/Numeric/Eproc/Common.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Eproc/Common.hs @@ -0,0 +1,70 @@ +{-# OPTIONS_HADDOCK prune #-} + +-- | +-- Module: Numeric.Eproc.Common +-- Copyright: (c) 2026 Jared Tobin +-- License: MIT +-- Maintainer: Jared Tobin +-- +-- 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) diff --git a/lib/Numeric/Eproc/Paired.hs b/lib/Numeric/Eproc/Paired.hs new file mode 100644 --- /dev/null +++ b/lib/Numeric/Eproc/Paired.hs @@ -0,0 +1,144 @@ +{-# OPTIONS_HADDOCK prune #-} +{-# LANGUAGE BangPatterns #-} + +-- | +-- Module: Numeric.Eproc.Paired +-- Copyright: (c) 2026 Jared Tobin +-- License: MIT +-- Maintainer: Jared Tobin +-- +-- 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 #-} diff --git a/ppad-eproc.cabal b/ppad-eproc.cabal new file mode 100644 --- /dev/null +++ b/ppad-eproc.cabal @@ -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 diff --git a/test/Main.hs b/test/Main.hs new file mode 100644 --- /dev/null +++ b/test/Main.hs @@ -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) + ]