Bernoulli.hs (10535B)
1 {-# OPTIONS_HADDOCK prune #-} 2 {-# LANGUAGE BangPatterns #-} 3 {-# LANGUAGE RecordWildCards #-} 4 5 -- | 6 -- Module: Numeric.Eproc.Bernoulli 7 -- Copyright: (c) 2026 Jared Tobin 8 -- License: MIT 9 -- Maintainer: Jared Tobin <jared@ppad.tech> 10 -- 11 -- One-sided Bernoulli rate anytime-valid test. See 12 -- "Numeric.Eproc.Bernoulli.TwoSided" for the two-sided companion 13 -- (used for the sign test at @p_0 = 1\/2@, among other things). 14 -- 15 -- For samples @x_t@ in @{0, 1}@, tests 16 -- 17 -- @H_0: E[x_t | F_{t-1}] <= p_0 for all t@ 18 -- 19 -- against @H_1: E[x_t | F_{t-1}] > p_0@ (at some @t@). Here 20 -- @F_{t-1}@ is the filtration generated by everything observed 21 -- strictly before time @t@; the conditional form is what anytime 22 -- validity actually requires. For i.i.d. samples this reduces to 23 -- the usual marginal statement @E[x] <= p_0@. 24 -- 25 -- A single wealth process is run: 26 -- 27 -- @W_n = prod_{i=1..n} (1 + lambda_i * (x_i - p_0))@ 28 -- 29 -- where each per-step bet @lambda_i@ is chosen predictably (from 30 -- data observed strictly before step @i@) and clipped to 31 -- @[0, lambda_max]@ so that the wealth factor stays nonnegative for 32 -- every admissible observation. Under @H_0@ the wealth process is 33 -- a nonnegative supermartingale, so by Ville's inequality the 34 -- probability of @W_n@ ever crossing @1 \/ alpha@ is at most 35 -- @alpha@, regardless of when the user decides to stop streaming 36 -- samples. Rejection is /latched/ in the running state: once the 37 -- wealth has crossed threshold, 'decide' continues to return 38 -- 'Reject' even if subsequent observations drive the current 39 -- wealth back below threshold. 40 -- 41 -- The alternative here is one-sided, so a single wealth process 42 -- suffices and no Bonferroni or hedge adjustment is needed -- the 43 -- rejection threshold is @log(1 \/ alpha)@. 44 -- 45 -- == Example 46 -- 47 -- Test @H_0: E[x] <= 0.05@ at level @alpha = 1e-3@ against a stream 48 -- with empirical rate @~0.5@: 49 -- 50 -- >>> let Right cfg = config 0.05 1.0e-3 Newton 51 -- >>> let xs = take 200 (cycle [True, False]) 52 -- >>> decide cfg (foldl' (update cfg) (initial cfg) xs) 53 -- Reject 54 55 module Numeric.Eproc.Bernoulli ( 56 -- * Test configuration and state 57 Config 58 , State 59 , Verdict(..) 60 , ConfigError(..) 61 62 -- * Bettor strategies 63 , Bettor(..) 64 65 -- * Construction 66 , config 67 , initial 68 69 -- * Streaming 70 , update 71 , decide 72 73 -- * Inspection 74 , log_wealth 75 , log_wealth_sup 76 , log_evalue 77 , log_evalue_sup 78 , p_value 79 , samples 80 ) where 81 82 import GHC.Float (log1p) 83 import Numeric.Eproc.Common ( 84 Bettor(..), Verdict(..), ConfigError(..) 85 , BetState, init_bet, bet_lambda, step_bet 86 , finite 87 ) 88 89 -- types ---------------------------------------------------------------------- 90 91 -- here, the centred observation @z_t@ referenced in 92 -- "Numeric.Eproc.Common" is @x_t - p_0@; the safe-bet ceiling 93 -- @lambda_max@ is derived from @p_0@ (see 'config'). 94 95 -- | Bernoulli rate test configuration. Build with 'config'. 96 -- 97 -- Carries the bettor strategy, the baseline rate, the significance 98 -- level, the precomputed log-wealth rejection threshold, and the 99 -- safe-bet ceiling derived from @p_0@. 100 data Config = Config { 101 -- ^ bettor strategy 102 cfg_bettor :: !Bettor 103 -- ^ safe-bet ceiling 104 , cfg_lam_max :: {-# UNPACK #-} !Double 105 -- ^ baseline rate @p_0@ 106 , cfg_p0 :: {-# UNPACK #-} !Double 107 -- ^ significance level @alpha@ 108 , cfg_alpha :: {-# UNPACK #-} !Double 109 -- ^ rejection threshold @log(1 \/ alpha)@ 110 , cfg_log_thresh :: {-# UNPACK #-} !Double 111 } 112 113 -- | Streaming test state. Construct with 'initial' and fold 114 -- observations through 'update'. 115 -- 116 -- Carries the sample count, current and supremum-so-far running 117 -- log-wealth, and whatever per-step state the chosen 'Bettor' 118 -- needs. The supremum field is what 'decide' tests against the 119 -- rejection threshold; this is the supremum-style event Ville's 120 -- inequality actually bounds. 121 data State = State { 122 st_n :: {-# UNPACK #-} !Int -- ^ sample count 123 , st_log_w :: {-# UNPACK #-} !Double -- ^ running log-wealth 124 , st_sup_log_w :: {-# UNPACK #-} !Double -- ^ sup log-wealth so far 125 , st_bet :: !BetState -- ^ bettor state 126 } 127 128 -- construction --------------------------------------------------------------- 129 130 -- | Build a 'Config' for the Bernoulli rate test. 131 -- 132 -- The safe-bet ceiling @lambda_max@ is set so that the wealth 133 -- factor @1 + lambda * (x - p_0)@ stays nonnegative for both 134 -- @x = 0@ and @x = 1@. The binding constraint is @x = 0@, which 135 -- requires @lambda <= 1 \/ p_0@; the ceiling stored is half this 136 -- to leave numerical margin -- the WSR safety recommendation. 137 -- 138 -- Returns 'Left' with a 'ConfigError' on inputs that would leave 139 -- the mathematical regime: either of @p_0@ or @alpha@ non-finite 140 -- (NaN or infinite); @p_0@ outside @(0, 1)@ (the degenerate case 141 -- @p_0 = 0@ would make @lambda_max@ infinite, and @p_0 = 1@ 142 -- leaves no room for an alternative); or @alpha@ outside 143 -- @(0, 1)@. 144 -- 145 -- >>> let Right cfg = config 0.05 1.0e-3 Newton 146 config 147 :: Double -- ^ baseline rate @p_0@, in @(0, 1)@ 148 -> Double -- ^ significance level @alpha@, in @(0, 1)@ 149 -> Bettor -- ^ bettor strategy 150 -> Either ConfigError Config 151 config !p0 !alpha !b 152 | not (finite p0 && p0 > 0 && p0 < 1) = 153 Left (InvalidBaselineRate p0) 154 | not (finite alpha && alpha > 0 && alpha < 1) = 155 Left (InvalidAlpha alpha) 156 | otherwise = Right Config { 157 cfg_bettor = b 158 , cfg_lam_max = 0.5 / p0 159 , cfg_p0 = p0 160 , cfg_alpha = alpha 161 , cfg_log_thresh = log (1 / alpha) 162 } 163 {-# INLINE config #-} 164 165 -- | The initial 'State' for a fresh streaming test. 166 -- 167 -- Both log-wealth fields start at @0@ (i.e., wealth @1@) and the 168 -- bettor starts in the per-strategy initial state appropriate 169 -- for the 'Bettor' chosen in the 'Config'. 170 -- 171 -- >>> let s0 = initial cfg 172 initial :: Config -> State 173 initial Config{..} = State { 174 st_n = 0 175 , st_log_w = 0 176 , st_sup_log_w = 0 177 , st_bet = init_bet cfg_bettor 178 } 179 {-# INLINE initial #-} 180 181 -- streaming ------------------------------------------------------------------ 182 183 -- | Fold one observation into the running 'State'. 184 -- 185 -- @True@ means @x_t = 1@ (the event of interest occurred -- e.g., 186 -- two readings diverged); @False@ means @x_t = 0@ (they matched). 187 -- The caller decides what \"matched\" means at the application 188 -- level. 189 -- 190 -- Computes the centred observation @z = x - p_0@, queries the 191 -- bettor for its predictable bet, accumulates log-wealth via 192 -- 193 -- @log_w' = log_w + log (1 + lambda * z)@ 194 -- 195 -- updates the running supremum log-wealth, then steps the bettor 196 -- state given the newly observed @z@. 197 -- 198 -- /Precondition/: @True@ and @False@ both /must/ be admissible 199 -- under the test (this holds vacuously for the @{0, 1}@ support). 200 -- The function is total. 201 -- 202 -- >>> let s1 = update cfg s0 True 203 update :: Config -> State -> Bool -> State 204 update Config{..} State{..} !x = 205 let !xd = if x then 1 else 0 206 !z = xd - cfg_p0 207 !lam = bet_lambda cfg_bettor cfg_lam_max st_bet 208 !logw' = st_log_w + log1p (lam * z) 209 !supw' = max st_sup_log_w logw' 210 !s' = step_bet cfg_bettor cfg_lam_max st_bet z 211 in State (st_n + 1) logw' supw' s' 212 {-# INLINE update #-} 213 214 -- | Compute the current 'Verdict' from the running 'State'. 215 -- 216 -- 'Reject' iff log-wealth has /ever/ crossed the threshold 217 -- @log(1 \/ alpha)@; equivalently, wealth has exceeded 218 -- @1 \/ alpha@ at some point in the stream so far. Under @H_0@, 219 -- by Ville's inequality, the probability of this ever happening 220 -- is at most @alpha@ -- and crucially this bound holds at /every/ 221 -- sample size simultaneously, so the user is free to peek at the 222 -- verdict as often as they like and stop on the first 'Reject'. 223 -- 224 -- >>> decide cfg s0 225 -- Continue 226 decide :: Config -> State -> Verdict 227 decide Config{..} State{..} 228 | st_sup_log_w >= cfg_log_thresh = Reject 229 | otherwise = Continue 230 {-# INLINE decide #-} 231 232 -- inspection ----------------------------------------------------------------- 233 234 -- | The current running log-wealth @log W_n@ at the present sample 235 -- count. 236 -- 237 -- Unlike 'log_wealth_sup' this is not monotone: adverse 238 -- observations decrease it. It is bounded above by 239 -- 'log_wealth_sup', which is what 'decide' tests against the 240 -- rejection threshold. 241 -- 242 -- >>> log_wealth s0 243 -- 0.0 244 log_wealth :: State -> Double 245 log_wealth = st_log_w 246 {-# INLINE log_wealth #-} 247 248 -- | The supremum-so-far log-wealth, across all sample counts up to 249 -- the current one. 250 -- 251 -- This is the natural \"test statistic\": it is monotone 252 -- nondecreasing in the sample count, and 'decide' rejects exactly 253 -- when it crosses @log(1 \/ alpha)@. 254 -- 255 -- >>> log_wealth_sup s0 256 -- 0.0 257 log_wealth_sup :: State -> Double 258 log_wealth_sup = st_sup_log_w 259 {-# INLINE log_wealth_sup #-} 260 261 -- | The current log e-value. For this one-sided test the single 262 -- wealth process is itself the e-process (a fresh state already 263 -- sits at wealth @1@), so this coincides with 'log_wealth'; the 264 -- accessor exists so that e-values read uniformly across test 265 -- modules regardless of their internal hedging, e.g. when 266 -- convex-combining several e-processes. Not monotone; bounded 267 -- above by 'log_evalue_sup'. 268 -- 269 -- >>> log_evalue s0 270 -- 0.0 271 log_evalue :: State -> Double 272 log_evalue = st_log_w 273 {-# INLINE log_evalue #-} 274 275 -- | The supremum-so-far of the log e-value; coincides with 276 -- 'log_wealth_sup' for this one-sided test. Monotone 277 -- nondecreasing, starting at @0@; 'decide' rejects exactly when 278 -- it crosses @log(1 \/ alpha)@. 279 -- 280 -- >>> log_evalue_sup s0 281 -- 0.0 282 log_evalue_sup :: State -> Double 283 log_evalue_sup = st_sup_log_w 284 {-# INLINE log_evalue_sup #-} 285 286 -- | The anytime-valid p-value: the reciprocal of the largest 287 -- e-value attained so far, @min 1 (exp (negate (log_evalue_sup 288 -- s)))@. 289 -- 290 -- Monotone nonincreasing in the sample count, and valid under 291 -- optional stopping: under @H_0@, 292 -- @P(exists t: p_t <= alpha) <= alpha@ for every @alpha@ 293 -- simultaneously. 'decide' returns 'Reject' exactly when this 294 -- value has reached the configured @alpha@ or below. 295 -- 296 -- >>> p_value s0 297 -- 1.0 298 p_value :: State -> Double 299 p_value s = min 1 (exp (negate (log_evalue_sup s))) 300 {-# INLINE p_value #-} 301 302 -- | The number of samples consumed so far. 303 -- 304 -- >>> samples s0 305 -- 0 306 samples :: State -> Int 307 samples = st_n 308 {-# INLINE samples #-}