eproc

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

commit f47c6310e1bc44f8a0213beeac799ec75cc4d2b6
Author: Jared Tobin <jared@jtobin.io>
Date:   Sun, 31 May 2026 21:15:11 -0230

Initial release.

Anytime-valid sequential testing via e-processes: bounded-mean and
paired two-sample tests built on the WSR betting framework, with
fixed-lambda, aGRAPA, and ONS bettor strategies.

Diffstat:
A.gitignore | 6++++++
ACHANGELOG | 6++++++
ALICENSE | 20++++++++++++++++++++
AREADME.md | 40++++++++++++++++++++++++++++++++++++++++
Aflake.lock | 88+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Aflake.nix | 59+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/Statistics/EProcess.hs | 123+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/Statistics/EProcess/Bettor.hs | 123+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/Statistics/EProcess/Mean.hs | 131+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Alib/Statistics/EProcess/TwoSample.hs | 82+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Appad-eproc.cabal | 58++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Atest/Main.hs | 201+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12 files changed, 937 insertions(+), 0 deletions(-)

diff --git a/.gitignore b/.gitignore @@ -0,0 +1,6 @@ +dist +dist-newstyle +result +result-doc +sandbox +.ghc.environment.* diff --git a/CHANGELOG b/CHANGELOG @@ -0,0 +1,6 @@ +# Changelog + +- 0.1.0 (2026-05-31) + * Initial release. Anytime-valid sequential testing via + e-processes: bounded-mean and paired two-sample tests, with + fixed-lambda, aGRAPA, and ONS bettors. diff --git a/LICENSE 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/README.md b/README.md @@ -0,0 +1,40 @@ +# ppad-eproc + +Anytime-valid sequential testing for Haskell, via e-processes and the +betting framework. + +Implements the bounded-mean and paired two-sample tests of Waudby-Smith +& Ramdas (2024) using predictable-plug-in bettors (fixed-lambda, +aGRAPA, ONS). Tests are valid under optional stopping: reject as soon +as the wealth process exceeds `1/alpha`, with type-I error controlled +at `alpha` regardless of when you stop. + +## Use + +```haskell +import qualified Statistics.EProcess as E + +-- Test H0: E[X] = 0.5 against H1: E[X] != 0.5, +-- samples bounded in [0, 1], alpha = 1e-6. +let cfg = E.meanConfig 0.5 0.0 1.0 1.0e-6 E.ons + s0 = E.initMeanState cfg + +-- Stream samples through the test: +let s1 = E.updateMean cfg s0 x1 + s2 = E.updateMean cfg s1 x2 + ... + +case E.decideMean cfg sN of + E.Reject -> ... -- H0 falsified + E.Continue -> ... -- more data needed +``` + +For paired two-sample testing, see `Statistics.EProcess.TwoSample`. + +## Background + +- Waudby-Smith & Ramdas (2024), "Estimating means of bounded random + variables by betting." JRSS-B. +- Ramdas, Grunwald, Vovk, Shafer (2023), "Game-theoretic statistics + and safe anytime-valid inference." Statistical Science. +- Shafer (2021), "Testing by betting." JRSS-A. diff --git a/flake.lock b/flake.lock @@ -0,0 +1,88 @@ +{ + "nodes": { + "flake-utils": { + "inputs": { + "systems": "systems" + }, + "locked": { + "lastModified": 1731533236, + "narHash": "sha256-l0KFg5HjrsfsO/JpG+r7fRrqm12kzFHyUHqHCVpMMbI=", + "owner": "numtide", + "repo": "flake-utils", + "rev": "11707dc2f618dd54ca8739b309ec4fc024de578b", + "type": "github" + }, + "original": { + "owner": "numtide", + "repo": "flake-utils", + "type": "github" + } + }, + "nixpkgs": { + "locked": { + "lastModified": 1766840161, + "narHash": "sha256-Ss/LHpJJsng8vz1Pe33RSGIWUOcqM1fjrehjUkdrWio=", + "owner": "NixOS", + "repo": "nixpkgs", + "rev": "3edc4a30ed3903fdf6f90c837f961fa6b49582d1", + "type": "github" + }, + "original": { + "owner": "NixOS", + "ref": "nixpkgs-unstable", + "repo": "nixpkgs", + "type": "github" + } + }, + "ppad-nixpkgs": { + "inputs": { + "flake-utils": "flake-utils", + "nixpkgs": "nixpkgs" + }, + "locked": { + "lastModified": 1766932084, + "narHash": "sha256-GvVsbTfW+B7IQ9K/QP2xcXJAm1lhBin1jYZWNjOzT+o=", + "ref": "master", + "rev": "353e61763b959b960a55321a85423501e3e9ed7a", + "revCount": 2, + "type": "git", + "url": "git://git.ppad.tech/nixpkgs.git" + }, + "original": { + "ref": "master", + "type": "git", + "url": "git://git.ppad.tech/nixpkgs.git" + } + }, + "root": { + "inputs": { + "flake-utils": [ + "ppad-nixpkgs", + "flake-utils" + ], + "nixpkgs": [ + "ppad-nixpkgs", + "nixpkgs" + ], + "ppad-nixpkgs": "ppad-nixpkgs" + } + }, + "systems": { + "locked": { + "lastModified": 1681028828, + "narHash": "sha256-Vy1rq5AaRuLzOxct8nz4T6wlgyUR7zLU309k9mBC768=", + "owner": "nix-systems", + "repo": "default", + "rev": "da67096a3b9bf56a91d16901293e51ba5b49a27e", + "type": "github" + }, + "original": { + "owner": "nix-systems", + "repo": "default", + "type": "github" + } + } + }, + "root": "root", + "version": 7 +} diff --git a/flake.nix b/flake.nix @@ -0,0 +1,59 @@ +{ + description = "Anytime-valid sequential testing via e-processes."; + + inputs = { + ppad-nixpkgs = { + type = "git"; + url = "git://git.ppad.tech/nixpkgs.git"; + ref = "master"; + }; + flake-utils.follows = "ppad-nixpkgs/flake-utils"; + nixpkgs.follows = "ppad-nixpkgs/nixpkgs"; + }; + + outputs = { self, nixpkgs, flake-utils, ppad-nixpkgs }: + flake-utils.lib.eachDefaultSystem (system: + let + lib = "ppad-eproc"; + + pkgs = import nixpkgs { inherit system; }; + hlib = pkgs.haskell.lib; + llvm = pkgs.llvmPackages_19.llvm; + clang = pkgs.llvmPackages_19.clang; + + hpkgs = pkgs.haskell.packages.ghc910.extend (new: old: { + ${lib} = new.callCabal2nix lib ./. {}; + }); + + cc = pkgs.stdenv.cc; + ghc = hpkgs.ghc; + cabal = hpkgs.cabal-install; + in + { + packages.default = hpkgs.${lib}; + + packages.haddock = hpkgs.${lib}.doc; + + devShells.default = hpkgs.shellFor { + packages = p: [ + p.${lib} + ]; + + buildInputs = [ + cabal + cc + llvm + ]; + + shellHook = '' + PS1="[${lib}] \w$ " + echo "entering ${system} shell, using" + echo "cc: $(${cc}/bin/cc --version)" + echo "ghc: $(${ghc}/bin/ghc --version)" + echo "cabal: $(${cabal}/bin/cabal --version)" + echo "llc: $(${llvm}/bin/llc --version | head -2 | tail -1)" + ''; + }; + } + ); +} diff --git a/lib/Statistics/EProcess.hs b/lib/Statistics/EProcess.hs @@ -0,0 +1,123 @@ +{-# OPTIONS_HADDOCK prune #-} + +-- | +-- Module: Statistics.EProcess +-- Copyright: (c) 2026 Jared Tobin +-- License: MIT +-- Maintainer: Jared Tobin <jared@ppad.tech> +-- +-- Anytime-valid sequential hypothesis testing for bounded random +-- variables, via the e-process / betting framework of Waudby-Smith +-- and Ramdas (2024). +-- +-- A bettor places predictable wagers against the null; the wealth +-- process is a nonnegative supermartingale under @H_0@, and Ville's +-- inequality gives type-I error control at @alpha@ for the stopping +-- rule \"reject the first time wealth exceeds @1\/alpha@\" — +-- regardless of when the user stops streaming samples. +-- +-- This module re-exports the primary API. For finer control, see: +-- +-- * "Statistics.EProcess.Bettor" for bettor strategies. +-- +-- * "Statistics.EProcess.Mean" for the one-sample bounded-mean +-- test. +-- +-- * "Statistics.EProcess.TwoSample" for the paired two-sample +-- mean-equality test. + +module Statistics.EProcess ( + -- * Bettors + Bettor + , AGRAPA + , ONS + , fixed + , agrapa + , ons + + -- * Bounded-mean test + -- + -- $mean + , Mean.Verdict(..) + , meanConfig + , initMeanState + , updateMean + , decideMean + + -- * Paired two-sample test + -- + -- $twosample + , twoSampleConfig + , initTwoSampleState + , updateTwoSample + , decideTwoSample + + -- * Inspection + , logWealth + , samples + ) where + +import Statistics.EProcess.Bettor +import qualified Statistics.EProcess.Mean as Mean +import qualified Statistics.EProcess.TwoSample as TS + +-- $mean +-- +-- For samples in @[lo, hi]@, test @H_0: E[x] = m@ two-sidedly. + +-- | See 'Mean.config'. +meanConfig + :: Double -- ^ null mean @m@ + -> Double -- ^ sample lower bound + -> Double -- ^ sample upper bound + -> Double -- ^ significance level @alpha@ + -> (Double -> Bettor s) + -> Mean.Config s +meanConfig = Mean.config + +-- | See 'Mean.initial'. +initMeanState :: Mean.Config s -> Mean.State s +initMeanState = Mean.initial + +-- | See 'Mean.update'. +updateMean :: Mean.Config s -> Mean.State s -> Double -> Mean.State s +updateMean = Mean.update + +-- | See 'Mean.decide'. +decideMean :: Mean.Config s -> Mean.State s -> Mean.Verdict +decideMean = Mean.decide + +-- $twosample +-- +-- For paired observations @(a, b)@ both in @[lo, hi]@, test @H_0: +-- E[a] = E[b]@ two-sidedly. + +-- | See 'TS.config'. +twoSampleConfig + :: Double + -> Double + -> Double + -> (Double -> Bettor s) + -> TS.Config s +twoSampleConfig = TS.config + +-- | See 'TS.initial'. +initTwoSampleState :: TS.Config s -> TS.State s +initTwoSampleState = TS.initial + +-- | See 'TS.update'. +updateTwoSample + :: TS.Config s -> TS.State s -> (Double, Double) -> TS.State s +updateTwoSample = TS.update + +-- | See 'TS.decide'. +decideTwoSample :: TS.Config s -> TS.State s -> TS.Verdict +decideTwoSample = TS.decide + +-- | Current log-wealth of a 'Mean.State'. +logWealth :: Mean.State s -> Double +logWealth = Mean.logWealth + +-- | Sample count consumed so far. +samples :: Mean.State s -> Int +samples = Mean.samples diff --git a/lib/Statistics/EProcess/Bettor.hs b/lib/Statistics/EProcess/Bettor.hs @@ -0,0 +1,123 @@ +{-# OPTIONS_HADDOCK prune #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE RecordWildCards #-} + +-- | +-- 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 maintains +-- internal state, consumes centred observations @z = x - m@, and +-- produces a predictable bet @lambda@ for the next observation. +-- +-- The bet placed at step @t@ depends only on data observed through +-- step @t-1@; this predictability is what makes the resulting wealth +-- process a nonnegative supermartingale under the null hypothesis, +-- and hence anytime-valid via Ville's inequality. + +module Statistics.EProcess.Bettor ( + -- * Bettor type + Bettor(..) + + -- * Strategies + , fixed + , agrapa + , ons + + -- * Strategy state types (opaque) + , AGRAPA + , ONS + ) where + +-- | A predictable bettor. +-- +-- Parameterised over its internal state type @s@. +-- +-- * @bettorInit@: initial state. +-- * @bettorStep@: update state with a newly observed centred +-- value @z = x - m@. +-- * @bettorBet@: bet @lambda@ to use for the /next/ observation, +-- given the current state. +data Bettor s = Bettor + { bettorInit :: !s + , bettorStep :: !(s -> Double -> s) + , bettorBet :: !(s -> Double) + } + +-- | Fixed-lambda bettor. +-- +-- Always bets the same value. Useful for smoke-testing the +-- framework and as a numerical baseline. +fixed :: Double -> Bettor () +fixed !lam = Bettor + { bettorInit = () + , bettorStep = \_ _ -> () + , bettorBet = \_ -> lam + } + +-- | aGRAPA bettor state (opaque). +data AGRAPA = AGRAPA + { aSum :: {-# UNPACK #-} !Double + , aSum2 :: {-# UNPACK #-} !Double + , aN :: {-# UNPACK #-} !Int + , aMax :: {-# UNPACK #-} !Double + } + +-- | aGRAPA (approximate growth-rate adaptive predictable plug-in). +-- +-- Tracks empirical mean and variance of centred observations @z@, +-- and bets the Kelly-optimal @lambda* = mu_z / (sigma_z^2 + mu_z^2)@ +-- given the current point estimate, clipped to @[0, lambda_max]@. +-- +-- The argument is @lambda_max@, the largest safe bet. For +-- observations @z = x - m@ where @x@ lies in @[lo, hi]@ and we are +-- testing @E[x] <= m@, a safe choice is @lambda_max = 1 \/ (m - lo)@ +-- (so that the wealth factor @1 + lambda * z@ stays nonnegative). +agrapa :: Double -> Bettor AGRAPA +agrapa !lamMax = Bettor + { bettorInit = AGRAPA 0 0 0 lamMax + , bettorStep = \AGRAPA{..} !z -> + AGRAPA (aSum + z) (aSum2 + z * z) (aN + 1) aMax + , bettorBet = \AGRAPA{..} -> + if aN == 0 + then 0 + else + let !n = fromIntegral aN + !mu = aSum / n + !mu2 = mu * mu + !var = max 0 (aSum2 / n - mu2) + !den = var + mu2 + !raw = if den == 0 then 0 else mu / den + in max 0 (min aMax raw) + } + +-- | ONS bettor state (opaque). +data ONS = ONS + { onsLambda :: {-# UNPACK #-} !Double + , onsAcc :: {-# UNPACK #-} !Double + , onsMax :: {-# UNPACK #-} !Double + } + +-- | ONS (online Newton step) bettor. +-- +-- Maintains a running sum of squared gradients of the per-step +-- log-wealth loss and updates @lambda@ by a Newton step at each +-- observation. Achieves logarithmic regret against the best +-- constant bet in hindsight; in practice the strongest of the +-- three bettors here under most signal regimes. +-- +-- The argument is @lambda_max@; see 'agrapa' for the sizing rule. +ons :: Double -> Bettor ONS +ons !lamMax = Bettor + { bettorInit = ONS 0 1.0e-6 lamMax + , bettorStep = \ONS{..} !z -> + let !denom = 1 + onsLambda * z + !g = if denom == 0 then 0 else negate z / denom + !acc' = onsAcc + g * g + !lam' = onsLambda - g / acc' + !clp = max 0 (min onsMax lam') + in ONS clp acc' onsMax + , bettorBet = onsLambda + } diff --git a/lib/Statistics/EProcess/Mean.hs b/lib/Statistics/EProcess/Mean.hs @@ -0,0 +1,131 @@ +{-# OPTIONS_HADDOCK prune #-} +{-# LANGUAGE BangPatterns #-} +{-# 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@. Runs two e-processes simultaneously (one per +-- direction) and combines them by Bonferroni: reject if either +-- side's wealth crosses @2 \/ alpha@. +-- +-- The test is anytime-valid: type-I error is controlled at @alpha@ +-- regardless of when the user stops streaming samples. + +module Statistics.EProcess.Mean ( + -- * Types + Config + , State + , Verdict(..) + + -- * Construction + , config + + -- * Streaming interface + , initial + , update + , decide + + -- * Inspection + , logWealth + , samples + ) where + +import Statistics.EProcess.Bettor + +-- | Test outcome at the current sample count. +data Verdict = Reject | Continue + deriving (Eq, Show) + +-- | Test configuration. Constructed by 'config'. +data Config s = Config + { cfgBetPos :: !(Bettor s) + , cfgBetNeg :: !(Bettor s) + , cfgNullMean :: {-# UNPACK #-} !Double + , cfgAlpha :: {-# UNPACK #-} !Double + , cfgLogThresh :: {-# UNPACK #-} !Double + } + +-- | Test state. Two log-wealth processes (one per direction) and +-- per-direction bettor state. +data State s = State + { stN :: {-# UNPACK #-} !Int + , stLogWPos :: {-# UNPACK #-} !Double + , stLogWNeg :: {-# UNPACK #-} !Double + , stBetPos :: !s + , stBetNeg :: !s + } + +-- | Build a test configuration. +-- +-- The bettor argument is a function from @lambda_max@ to a bettor; +-- the same builder is used for both directions, with appropriate +-- bounds computed from @lo@, @hi@, and @m@. +-- +-- >>> import qualified Statistics.EProcess.Bettor as B +-- >>> let cfg = config 0.5 0.0 1.0 1.0e-6 B.ons +config + :: Double -- ^ null mean @m@ + -> Double -- ^ sample lower bound @lo@ + -> Double -- ^ sample upper bound @hi@ + -> Double -- ^ significance level @alpha@ + -> (Double -> Bettor s) -- ^ bettor builder, taking @lambda_max@ + -> Config s +config !m !lo !hi !alpha mk = Config + { cfgBetPos = mk (0.5 / (m - lo)) + , cfgBetNeg = mk (0.5 / (hi - m)) + , cfgNullMean = m + , cfgAlpha = alpha + , cfgLogThresh = log (2 / alpha) + } +-- NB. argument to @mk@ is half the geometric @lambda_max@; the 1/2 +-- margin keeps the wealth factor bounded away from zero at the +-- boundary, which is the WSR safety recommendation. + +-- | Initial state for streaming. +initial :: Config s -> State s +initial Config{..} = State + { stN = 0 + , stLogWPos = 0 + , stLogWNeg = 0 + , stBetPos = bettorInit cfgBetPos + , stBetNeg = bettorInit cfgBetNeg + } + +-- | Fold one observation into the state. +update :: Config s -> State s -> Double -> State s +update Config{..} State{..} !x = + let !z = x - cfgNullMean + !lamP = bettorBet cfgBetPos stBetPos + !lamN = bettorBet cfgBetNeg stBetNeg + !facP = 1 + lamP * z + !facN = 1 - lamN * z + !logWP' = stLogWPos + log (max 1.0e-300 facP) + !logWN' = stLogWNeg + log (max 1.0e-300 facN) + !sP' = bettorStep cfgBetPos stBetPos z + !sN' = bettorStep cfgBetNeg stBetNeg (negate z) + in State (stN + 1) logWP' logWN' sP' sN' + +-- | Decide based on current wealth. +-- +-- 'Reject' iff either directional log-wealth has crossed the +-- Bonferroni-adjusted threshold @log(2 \/ alpha)@. +decide :: Config s -> State s -> Verdict +decide Config{..} State{..} + | stLogWPos >= cfgLogThresh = Reject + | stLogWNeg >= cfgLogThresh = Reject + | otherwise = Continue + +-- | Current log-wealth (the larger of the two directional processes). +logWealth :: State s -> Double +logWealth State{..} = max stLogWPos stLogWNeg + +-- | Sample count consumed so far. +samples :: State s -> Int +samples = stN diff --git a/lib/Statistics/EProcess/TwoSample.hs b/lib/Statistics/EProcess/TwoSample.hs @@ -0,0 +1,82 @@ +{-# 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]@ +-- by running the bounded-mean test on the differences @d_t = a_t - +-- b_t@ with null mean 0. + +module Statistics.EProcess.TwoSample ( + -- * Types + Config + , State + , Verdict(..) + + -- * Construction + , config + + -- * Streaming interface + , initial + , update + , decide + + -- * Inspection + , logWealth + , samples + ) where + +import qualified Statistics.EProcess.Mean as M +import Statistics.EProcess.Mean (Verdict(..)) +import Statistics.EProcess.Bettor (Bettor) + +-- | Test configuration. +newtype Config s = Config (M.Config s) + +-- | Test state. +newtype State s = State (M.State s) + +-- | Build a paired two-sample test configuration. +-- +-- Bounds @lo@ and @hi@ are the (shared) bounds on the individual +-- samples; differences then lie in @[lo - hi, hi - lo]@. +-- +-- >>> import qualified Statistics.EProcess.Bettor as B +-- >>> let cfg = config 0.0 1.0 1.0e-6 B.ons +config + :: Double -- ^ sample lower bound @lo@ + -> Double -- ^ sample upper bound @hi@ + -> Double -- ^ significance level @alpha@ + -> (Double -> Bettor s) -- ^ bettor builder + -> Config s +config !lo !hi !alpha mk = + let !b = hi - lo + in Config (M.config 0 (negate b) b alpha mk) + +-- | Initial state for streaming. +initial :: Config s -> State s +initial (Config c) = State (M.initial c) + +-- | Fold one paired observation @(a, b)@ into the state. +update :: Config s -> State s -> (Double, Double) -> State s +update (Config c) (State s) (!a, !b) = + State (M.update c s (a - b)) + +-- | Decide based on current wealth. +decide :: Config s -> State s -> Verdict +decide (Config c) (State s) = M.decide c s + +-- | Current log-wealth. +logWealth :: State s -> Double +logWealth (State s) = M.logWealth s + +-- | Sample count consumed so far. +samples :: State s -> Int +samples (State s) = M.samples s diff --git a/ppad-eproc.cabal b/ppad-eproc.cabal @@ -0,0 +1,58 @@ +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 and paired two-sample tests + with fixed-lambda, aGRAPA, and ONS bettors. Tests are valid under + optional stopping. + +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: + Statistics.EProcess + Statistics.EProcess.Bettor + Statistics.EProcess.Mean + Statistics.EProcess.TwoSample + 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 diff --git a/test/Main.hs b/test/Main.hs @@ -0,0 +1,201 @@ +{-# LANGUAGE BangPatterns #-} + +module Main where + +import Data.Bits +import Data.List (foldl') +import Data.Word +import qualified Statistics.EProcess as E +import qualified Statistics.EProcess.Mean as M +import qualified Statistics.EProcess.TwoSample as TS +import Test.Tasty +import Test.Tasty.HUnit + +main :: IO () +main = defaultMain $ testGroup "ppad-eproc" [ + sanityTests + , calibrationTests + , powerTests + , twoSampleTests + , bettorSmokeTests + ] + +-- inline PCG-style PRNG, no external deps. + +newtype Gen = Gen Word64 + +mkGen :: Word64 -> Gen +mkGen = Gen + +stepGen :: Gen -> (Word64, Gen) +stepGen (Gen s) = + let !s' = s * 6364136223846793005 + 1442695040888963407 + in (s', Gen s') + +nextDouble :: Gen -> (Double, Gen) +nextDouble g = + let (w, g') = stepGen g + !x = fromIntegral (w `shiftR` 11 .&. 0x1FFFFFFFFFFFFF) / + 9007199254740992 + in (x, g') + +bernoulli :: Double -> Gen -> (Double, Gen) +bernoulli !p g = + let (u, g') = nextDouble g + in (if u < p then 1.0 else 0.0, g') + +-- run a sequential mean test on a stream of n bernoulli(p) samples, +-- with the early-stopping rule built in. returns (verdict, samples +-- consumed). +runMeanBernoulli + :: M.Config s + -> Double -- ^ p + -> Int -- ^ budget + -> Gen + -> (M.Verdict, Int) +runMeanBernoulli cfg p budget g0 = go 0 g0 (M.initial cfg) + where + go !n !g !st + | n >= budget = (M.decide cfg st, n) + | otherwise = case M.decide cfg st of + M.Reject -> (M.Reject, n) + M.Continue -> + let (x, g') = bernoulli p g + st' = M.update cfg st x + in go (n + 1) g' st' + +-- fraction of trials that rejected. +rejectionRate + :: M.Config s + -> Double -- ^ true bernoulli p + -> Int -- ^ budget per trial + -> Int -- ^ number of trials + -> Word64 -- ^ seed + -> Double +rejectionRate cfg p budget trials seed = + let gens = take trials (genSeq (mkGen seed)) + rejects = length + [ () | g <- gens + , let (v, _) = runMeanBernoulli cfg p budget g + , v == M.Reject ] + in fromIntegral rejects / fromIntegral trials + +genSeq :: Gen -> [Gen] +genSeq g = let (_, g') = stepGen g in g : genSeq g' + +-- sanity: with all-zero deviations from the null mean, no rejection. + +sanityTests :: TestTree +sanityTests = testGroup "sanity" [ + testCase "degenerate input never rejects" $ do + let cfg = M.config 0.5 0.0 1.0 1.0e-6 E.ons :: M.Config E.ONS + xs = replicate 5000 0.5 + st = foldl' (M.update cfg) (M.initial cfg) xs + M.decide cfg st @?= M.Continue + , testCase "two-sided thresholds applied symmetrically" $ do + let cfg = M.config 0.5 0.0 1.0 1.0e-6 E.ons :: M.Config E.ONS + M.decide cfg (M.initial cfg) @?= M.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. + +calibrationTests :: TestTree +calibrationTests = testGroup "null calibration" [ + testCase "ONS, Bernoulli(0.5), m=0.5, alpha=0.05" $ do + let cfg = M.config 0.5 0.0 1.0 0.05 E.ons :: M.Config E.ONS + rate = rejectionRate cfg 0.5 2000 200 12345 + -- expected rate ≤ 0.05; allow up to 0.10 slack for sampling + -- variability over 200 trials. + assertBool ("FPR " ++ show rate ++ " exceeded slack") $ + rate <= 0.10 + , testCase "aGRAPA, Bernoulli(0.5), m=0.5, alpha=0.05" $ do + let cfg = M.config 0.5 0.0 1.0 0.05 E.agrapa :: M.Config E.AGRAPA + rate = rejectionRate 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. + +powerTests :: TestTree +powerTests = testGroup "power" [ + testCase "ONS detects Bernoulli(0.7) vs m=0.5" $ do + let cfg = M.config 0.5 0.0 1.0 1.0e-3 E.ons :: M.Config E.ONS + rate = rejectionRate cfg 0.7 5000 100 11111 + assertBool ("power " ++ show rate ++ " too low") $ + rate >= 0.95 + , testCase "aGRAPA detects Bernoulli(0.7) vs m=0.5" $ do + let cfg = M.config 0.5 0.0 1.0 1.0e-3 E.agrapa + :: M.Config E.AGRAPA + rate = rejectionRate cfg 0.7 5000 100 22222 + assertBool ("power " ++ show rate ++ " too low") $ + rate >= 0.95 + ] + +-- two-sample paired test. + +runTSPaired + :: TS.Config s + -> Double + -> Double -- ^ p for A and B + -> Int + -> Gen + -> (TS.Verdict, Int) +runTSPaired cfg pA pB budget g0 = go 0 g0 (TS.initial cfg) + where + go !n !g !st + | n >= budget = (TS.decide cfg st, n) + | otherwise = case TS.decide cfg st of + M.Reject -> (M.Reject, n) + M.Continue -> + let (a, g1) = bernoulli pA g + (b, g2) = bernoulli pB g1 + st' = TS.update cfg st (a, b) + in go (n + 1) g2 st' + +twoSampleTests :: TestTree +twoSampleTests = testGroup "two-sample" [ + testCase "identical distributions don't reject" $ do + let cfg = TS.config 0.0 1.0 1.0e-3 E.ons :: TS.Config E.ONS + rate = avgRate cfg 0.5 0.5 2000 100 33333 + assertBool ("FPR " ++ show rate) $ rate <= 0.05 + , testCase "different distributions reject" $ do + let cfg = TS.config 0.0 1.0 1.0e-3 E.ons :: TS.Config E.ONS + rate = avgRate cfg 0.3 0.7 5000 100 44444 + assertBool ("power " ++ show rate) $ rate >= 0.95 + ] + where + avgRate cfg pA pB budget trials seed = + let gens = take trials (genSeq (mkGen seed)) + rejects = length + [ () | g <- gens + , let (v, _) = runTSPaired cfg pA pB budget g + , v == M.Reject ] + in fromIntegral rejects / fromIntegral trials + +-- bettor smoke tests: each bettor produces a well-defined state and +-- decision when run on a small deterministic stream. + +bettorSmokeTests :: TestTree +bettorSmokeTests = testGroup "bettor smoke" [ + testCase "fixed bettor runs without error" $ do + let cfg = M.config 0.5 0.0 1.0 1.0e-3 + (const (E.fixed 0.5)) :: M.Config () + xs = take 100 (cycle [0.0, 1.0]) + st = foldl' (M.update cfg) (M.initial cfg) xs + assertBool "samples advanced" (M.samples st == 100) + , testCase "ONS bettor runs without error" $ do + let cfg = M.config 0.5 0.0 1.0 1.0e-3 E.ons :: M.Config E.ONS + xs = take 100 (cycle [0.0, 1.0]) + st = foldl' (M.update cfg) (M.initial cfg) xs + assertBool "samples advanced" (M.samples st == 100) + , testCase "aGRAPA bettor runs without error" $ do + let cfg = M.config 0.5 0.0 1.0 1.0e-3 E.agrapa + :: M.Config E.AGRAPA + xs = take 100 (cycle [0.0, 1.0]) + st = foldl' (M.update cfg) (M.initial cfg) xs + assertBool "samples advanced" (M.samples st == 100) + ]