eproc

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

commit f84c84ac03843da025ad73533629af10e6867691
parent e168b262bac7bd5ae011314dab46c37ffee12055
Author: Jared Tobin <jared@jtobin.io>
Date:   Wed,  3 Jun 2026 10:20:37 -0230

rename Statistics.EProcess.* to Numeric.Eproc.*

Module path rename only. Each library module moves to its new path
(git tracks the renames), module declarations and internal cross-refs
updated, Haddock cross-links ('Statistics.EProcess.Mean.config' etc.)
follow. ppad-eproc.cabal exposed-modules, test/Main.hs imports,
bench/Main.hs and bench/Weight.hs imports, and the README's GHCi
example all updated in lockstep.

No semantic / API changes; all 11 tests still pass.

Diffstat:
MREADME.md | 6+++---
Mbench/Main.hs | 6+++---
Mbench/Weight.hs | 6+++---
Alib/Numeric/Eproc/Bettor.hs | 62++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/Numeric/Eproc/Mean.hs | 310+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/Numeric/Eproc/TwoSample.hs | 132+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Dlib/Statistics/EProcess/Bettor.hs | 69---------------------------------------------------------------------
Dlib/Statistics/EProcess/Mean.hs | 310-------------------------------------------------------------------------------
Dlib/Statistics/EProcess/TwoSample.hs | 132-------------------------------------------------------------------------------
Mppad-eproc.cabal | 6+++---
Mtest/Main.hs | 6+++---
11 files changed, 519 insertions(+), 526 deletions(-)

diff --git a/README.md b/README.md @@ -17,8 +17,8 @@ A sample GHCi session: ``` > -- import qualified - > import qualified Statistics.EProcess.Bettor as B - > import qualified Statistics.EProcess.Mean as M + > import qualified Numeric.Eproc.Bettor as B + > import qualified Numeric.Eproc.Mean as M > > -- test H_0: E[X] = 0.5 for samples in [0, 1] at alpha = 1e-3, > -- with the ONS bettor @@ -46,7 +46,7 @@ A sample GHCi session: ``` For the paired two-sample mean-equality test, see -`Statistics.EProcess.TwoSample`. +`Numeric.Eproc.TwoSample`. ## Documentation diff --git a/bench/Main.hs b/bench/Main.hs @@ -4,9 +4,9 @@ module Main where import Control.DeepSeq -import qualified Statistics.EProcess.Bettor as B -import qualified Statistics.EProcess.Mean as M -import qualified Statistics.EProcess.TwoSample as TS +import qualified Numeric.Eproc.Bettor as B +import qualified Numeric.Eproc.Mean as M +import qualified Numeric.Eproc.TwoSample as TS import Criterion.Main -- all relevant fields are strict (and UNPACK'd for the doubles), so diff --git a/bench/Weight.hs b/bench/Weight.hs @@ -4,9 +4,9 @@ module Main where import Control.DeepSeq -import qualified Statistics.EProcess.Bettor as B -import qualified Statistics.EProcess.Mean as M -import qualified Statistics.EProcess.TwoSample as TS +import qualified Numeric.Eproc.Bettor as B +import qualified Numeric.Eproc.Mean as M +import qualified Numeric.Eproc.TwoSample as TS import Weigh instance NFData M.State where rnf !_ = () diff --git a/lib/Numeric/Eproc/Bettor.hs b/lib/Numeric/Eproc/Bettor.hs @@ -0,0 +1,62 @@ +{-# OPTIONS_HADDOCK prune #-} +{-# LANGUAGE BangPatterns #-} + +-- | +-- Module: Numeric.Eproc.Bettor +-- Copyright: (c) 2026 Jared Tobin +-- License: MIT +-- Maintainer: Jared Tobin <jared@ppad.tech> +-- +-- Bettor strategies for the e-process framework. +-- +-- A bettor describes how, given the history of centred observations +-- @z_t = x_t - m@ (where @x_t@ is the new observation and @m@ is the +-- null mean), the next predictable bet @lambda_t@ is chosen. The +-- wealth process is the running product of per-step factors +-- +-- @W_t = prod_{s <= t} (1 + lambda_s * z_s)@ +-- +-- and the test rejects when @W_t@ crosses @1\/alpha@. Predictability +-- -- that is, @lambda_t@ depends only on data observed strictly +-- before step @t@ -- is what makes @W@ a nonnegative supermartingale +-- under @H_0@, so that Ville's inequality applies and the resulting +-- test is anytime-valid. + +module Numeric.Eproc.Bettor ( + -- * Bettor strategies + Bettor(..) + ) where + +-- bettor strategies ---------------------------------------------------------- + +-- | A predictable 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 +-- @[0, lambda_max]@ so that the wealth factor @1 + lambda * z@ +-- stays nonnegative for every admissible observation. +-- +-- * 'Fixed' always bets the supplied constant @lambda@. The wager +-- does not respond to observed data; this strategy is useful only +-- as a baseline. +-- +-- * 'Agrapa' 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. +-- +-- * 'Ons' is the online Newton step 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 + | Agrapa + | Ons + deriving (Eq, Show) diff --git a/lib/Numeric/Eproc/Mean.hs b/lib/Numeric/Eproc/Mean.hs @@ -0,0 +1,310 @@ +{-# 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/TwoSample.hs b/lib/Numeric/Eproc/TwoSample.hs @@ -0,0 +1,132 @@ +{-# OPTIONS_HADDOCK prune #-} +{-# LANGUAGE BangPatterns #-} + +-- | +-- Module: Numeric.Eproc.TwoSample +-- 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.Mean") 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. + +module Numeric.Eproc.TwoSample ( + -- * Test configuration and state + Config + , State + , Verdict(..) + + -- * Construction + , config + , initial + + -- * Streaming + , update + , decide + + -- * Inspection + , log_wealth + , samples + ) where + +import qualified Numeric.Eproc.Mean as M +import Numeric.Eproc.Mean (Verdict(..)) +import Numeric.Eproc.Bettor (Bettor) + +-- types ---------------------------------------------------------------------- + +-- | Paired two-sample test configuration. Build with 'config'. Wraps +-- a 'Numeric.Eproc.Mean.Config' for the underlying +-- difference test. +newtype Config = Config M.Config + +-- | Streaming paired two-sample test state. Construct with 'initial' +-- and fold paired observations through 'update'. +newtype State = State M.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@. +-- +-- >>> import qualified Numeric.Eproc.Bettor as B +-- >>> let cfg = config 0.0 1.0 1.0e-3 B.Ons +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 (M.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) +{-# 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 (M.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) = M.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) = M.log_wealth s +{-# INLINE log_wealth #-} + +-- | The number of paired observations consumed so far. +-- +-- >>> samples s0 +-- 0 +samples :: State -> Int +samples (State s) = M.samples s +{-# INLINE samples #-} diff --git a/lib/Statistics/EProcess/Bettor.hs b/lib/Statistics/EProcess/Bettor.hs @@ -1,69 +0,0 @@ -{-# OPTIONS_HADDOCK prune #-} -{-# LANGUAGE BangPatterns #-} - --- | --- Module: Statistics.EProcess.Bettor --- Copyright: (c) 2026 Jared Tobin --- License: MIT --- Maintainer: Jared Tobin <jared@ppad.tech> --- --- Bettor strategies for the e-process framework. --- --- A bettor describes how, given the history of centred observations --- @z_t = x_t - m@ (where @x_t@ is the new observation and @m@ is the --- null mean), the next predictable bet @lambda_t@ is chosen. The --- wealth process is the running product of per-step factors --- --- @W_t = prod_{s <= t} (1 + lambda_s * z_s)@ --- --- and the test rejects when @W_t@ crosses @1\/alpha@. Predictability --- -- that is, @lambda_t@ depends only on data observed strictly --- before step @t@ -- is what makes @W@ a nonnegative supermartingale --- under @H_0@, so that Ville's inequality applies and the resulting --- test is anytime-valid. - -module Statistics.EProcess.Bettor ( - -- * Bettor strategies - Bettor(..) - ) where - --- bettor strategies ---------------------------------------------------------- - --- | A predictable 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. --- 'Statistics.EProcess.Mean.config') -- bets get clipped to --- @[0, lambda_max]@ so that the wealth factor @1 + lambda * z@ --- stays nonnegative for every admissible observation. --- --- * 'Fixed' always bets the supplied constant @lambda@. The wager --- does not respond to observed data; this strategy is useful only --- as a baseline or for smoke testing the framework. --- --- * 'Agrapa' 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. --- --- * 'Ons' is the online Newton step 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. --- --- >>> Fixed 0.5 --- Fixed 0.5 --- >>> Agrapa --- Agrapa --- >>> Ons --- Ons -data Bettor = - Fixed {-# UNPACK #-} !Double - | Agrapa - | Ons - deriving (Eq, Show) diff --git a/lib/Statistics/EProcess/Mean.hs b/lib/Statistics/EProcess/Mean.hs @@ -1,310 +0,0 @@ -{-# OPTIONS_HADDOCK prune #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE MagicHash #-} -{-# LANGUAGE RecordWildCards #-} - --- | --- Module: Statistics.EProcess.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 Statistics.EProcess.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 Statistics.EProcess.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 Statistics.EProcess.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/Statistics/EProcess/TwoSample.hs b/lib/Statistics/EProcess/TwoSample.hs @@ -1,132 +0,0 @@ -{-# OPTIONS_HADDOCK prune #-} -{-# LANGUAGE BangPatterns #-} - --- | --- Module: Statistics.EProcess.TwoSample --- 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 ("Statistics.EProcess.Mean") 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. - -module Statistics.EProcess.TwoSample ( - -- * Test configuration and state - Config - , State - , Verdict(..) - - -- * Construction - , config - , initial - - -- * Streaming - , update - , decide - - -- * Inspection - , log_wealth - , samples - ) where - -import qualified Statistics.EProcess.Mean as M -import Statistics.EProcess.Mean (Verdict(..)) -import Statistics.EProcess.Bettor (Bettor) - --- types ---------------------------------------------------------------------- - --- | Paired two-sample test configuration. Build with 'config'. Wraps --- a 'Statistics.EProcess.Mean.Config' for the underlying --- difference test. -newtype Config = Config M.Config - --- | Streaming paired two-sample test state. Construct with 'initial' --- and fold paired observations through 'update'. -newtype State = State M.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@. --- --- >>> import qualified Statistics.EProcess.Bettor as B --- >>> let cfg = config 0.0 1.0 1.0e-3 B.Ons -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 (M.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) -{-# 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 (M.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) = M.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) = M.log_wealth s -{-# INLINE log_wealth #-} - --- | The number of paired observations consumed so far. --- --- >>> samples s0 --- 0 -samples :: State -> Int -samples (State s) = M.samples s -{-# INLINE samples #-} diff --git a/ppad-eproc.cabal b/ppad-eproc.cabal @@ -34,9 +34,9 @@ library if flag(llvm) ghc-options: -fllvm -O2 exposed-modules: - Statistics.EProcess.Bettor - Statistics.EProcess.Mean - Statistics.EProcess.TwoSample + Numeric.Eproc.Bettor + Numeric.Eproc.Mean + Numeric.Eproc.TwoSample build-depends: base >= 4.9 && < 5 diff --git a/test/Main.hs b/test/Main.hs @@ -4,9 +4,9 @@ module Main where import Data.Bits import Data.Word -import qualified Statistics.EProcess.Bettor as B -import qualified Statistics.EProcess.Mean as M -import qualified Statistics.EProcess.TwoSample as TS +import qualified Numeric.Eproc.Bettor as B +import qualified Numeric.Eproc.Mean as M +import qualified Numeric.Eproc.TwoSample as TS import Test.Tasty import Test.Tasty.HUnit