Bounded.hs (17317B)
1 {-# OPTIONS_HADDOCK prune #-} 2 {-# LANGUAGE BangPatterns #-} 3 {-# LANGUAGE RecordWildCards #-} 4 5 -- | 6 -- Module: Numeric.Eproc.Bounded 7 -- Copyright: (c) 2026 Jared Tobin 8 -- License: MIT 9 -- Maintainer: Jared Tobin <jared@ppad.tech> 10 -- 11 -- Two-sided bounded-mean anytime-valid test. 12 -- 13 -- For samples @x_t@ in @[lo, hi]@, tests 14 -- 15 -- @H_0: E[x_t | F_{t-1}] = m for all t@ 16 -- 17 -- against the negation. Here @F_{t-1}@ is the filtration generated 18 -- by everything observed strictly before time @t@; the conditional 19 -- form is what anytime validity actually requires. For i.i.d. 20 -- samples this reduces to the usual marginal statement 21 -- @E[x] = m@; for adaptively-collected or otherwise non-i.i.d. 22 -- streams the conditional statement is the right thing to think 23 -- about. 24 -- 25 -- 'configInterval' generalises the null from a point to an 26 -- interval: 27 -- 28 -- @H_0: m_lo <= E[x_t | F_{t-1}] <= m_hi for all t@ 29 -- 30 -- The positive-direction process is centred at @m_hi@ and the 31 -- negative-direction process at @m_lo@, so each side's wealth is a 32 -- nonnegative supermartingale under /every/ conditional mean in 33 -- the interval, and the convex-hedge guarantee carries over 34 -- unchanged. The point null is the degenerate case 35 -- @m_lo = m_hi@. An interval null buys tolerance: a systematic 36 -- effect smaller than the interval half-width accumulates no 37 -- evidence, which is the honest null when the measurement channel 38 -- itself carries a small, bounded systematic that is not the 39 -- effect under test. 40 -- 41 -- Internally two one-sided e-processes are run in parallel: a 42 -- /positive-direction/ process @K^+_t@ betting against the 43 -- alternative @E[x_t | F_{t-1}] > m@ (using centred observations 44 -- @z = x - m@), and a /negative-direction/ process @K^-_t@ betting 45 -- against @E[x_t | F_{t-1}] < m@ (using @-z@). Each maintains its 46 -- own log-wealth and bettor state. 47 -- 48 -- The two sides are combined via the /hedged capital process/ of 49 -- Waudby-Smith & Ramdas (2024) §4: their average 50 -- @K_t = (K^+_t + K^-_t) \/ 2@ is itself an e-process (convex 51 -- combinations preserve the supermartingale property), with 52 -- @E[K_0] = 1@. By Ville's inequality 53 -- @P(sup_t K_t >= 1 \/ alpha) <= alpha@, so the test rejects when 54 -- the supremum of @K^+_t + K^-_t@ has ever crossed @2 \/ alpha@. 55 -- 56 -- This is strictly more powerful than the naive Bonferroni union 57 -- (reject when @max(K^+_t, K^-_t) >= 2 \/ alpha@): the convex-hedge 58 -- rejection region contains Bonferroni's (since 59 -- @K^+ + K^- >= max(K^+, K^-)@), with the same alpha guarantee. 60 -- For one-sided alternatives the gap is small (the losing-direction 61 -- bettor stays near @1@); for genuinely two-sided alternatives it 62 -- can be substantial. 63 -- 64 -- The test is /anytime-valid/: under @H_0@ the wealth process is a 65 -- nonnegative supermartingale, so by Ville's inequality the 66 -- probability of /ever/ crossing the threshold is at most @alpha@, 67 -- regardless of when the user decides to stop streaming samples. 68 -- Rejection is /latched/ in the running state -- once a side has 69 -- crossed threshold, 'decide' continues to return 'Reject' even if 70 -- the current log-wealth has since dropped back below threshold. 71 -- 72 -- == Example 73 -- 74 -- Test @H_0: E[x] = 0.5@ for @x@ in @[0, 1]@ at level @alpha = 1e-3@ 75 -- against a stream with empirical mean @0.8@: 76 -- 77 -- >>> let Right cfg = config 0.5 0.0 1.0 1.0e-3 Newton 78 -- >>> let xs = concat (replicate 30 [1, 1, 0, 1, 1, 0, 1, 1, 1, 1]) 79 -- >>> decide cfg (foldl' (update cfg) (initial cfg) xs) 80 -- Reject 81 82 module Numeric.Eproc.Bounded ( 83 -- * Test configuration and state 84 Config 85 , State 86 , Verdict(..) 87 , ConfigError(..) 88 89 -- * Bettor strategies 90 , Bettor(..) 91 92 -- * Construction 93 , config 94 , configInterval 95 , initial 96 97 -- * Streaming 98 , update 99 , decide 100 101 -- * Inspection 102 , log_wealth 103 , log_wealth_sup 104 , log_evalue 105 , log_evalue_sup 106 , p_value 107 , samples 108 ) where 109 110 import GHC.Float (log1p) 111 import Numeric.Eproc.Common ( 112 Bettor(..), Verdict(..), ConfigError(..) 113 , BetState, init_bet, bet_lambda, step_bet 114 , finite, log_sum_exp, log2_dbl 115 ) 116 117 -- types ---------------------------------------------------------------------- 118 119 -- here, the centred observation @z_t@ referenced in 120 -- "Numeric.Eproc.Common" is per-direction: @x_t - m_hi@ for the 121 -- positive side and @x_t - m_lo@ for the negative side, with 122 -- @m_lo = m_hi@ under a point null ('config'). The per-direction 123 -- safe-bet ceilings @lambda_max@ are derived from the sample 124 -- bounds (see 'config' \/ 'configInterval'). 125 126 -- | Bounded-mean test configuration. Build with 'config' (point 127 -- null) or 'configInterval' (interval null). 128 -- 129 -- Carries the bettor strategy, the null-mean endpoints (equal 130 -- under a point null), the significance level, the precomputed 131 -- convex-hedge log-wealth threshold, and the per-direction 132 -- safe-bet ceilings (see 'config' for how the latter are derived 133 -- from the sample bounds). 134 data Config = Config { 135 -- ^ bettor strategy 136 cfg_bettor :: !Bettor 137 -- ^ positive-direction safe-bet ceiling 138 , cfg_lam_max_pos :: {-# UNPACK #-} !Double 139 -- ^ negative-direction safe-bet ceiling 140 , cfg_lam_max_neg :: {-# UNPACK #-} !Double 141 -- ^ lower null mean @m_lo@ (negative-direction centre) 142 , cfg_null_lo :: {-# UNPACK #-} !Double 143 -- ^ upper null mean @m_hi@ (positive-direction centre) 144 , cfg_null_hi :: {-# UNPACK #-} !Double 145 -- ^ significance level @alpha@ 146 , cfg_alpha :: {-# UNPACK #-} !Double 147 -- ^ rejection threshold @log(2 \/ alpha)@ 148 , cfg_log_thresh :: {-# UNPACK #-} !Double 149 } 150 151 -- | Streaming test state. Construct with 'initial' and fold 152 -- observations through 'update'. 153 -- 154 -- The two log-wealth fields track the running log-wealth of the 155 -- positive- and negative-direction e-processes separately; the 156 -- /sup log-sum/ field latches the supremum so far of 157 -- @log(K^+_t + K^-_t)@, which is the test statistic the 158 -- convex-hedge construction actually monitors. The per-direction 159 -- bettor states carry whatever the chosen 'Bettor' needs (running 160 -- sums, current bet, etc.). 161 data State = State { 162 st_n :: {-# UNPACK #-} !Int -- ^ sample count 163 , st_log_w_pos :: {-# UNPACK #-} !Double -- ^ log-wealth, pos 164 , st_log_w_neg :: {-# UNPACK #-} !Double -- ^ log-wealth, neg 165 , st_sup_log_sum :: {-# UNPACK #-} !Double -- ^ sup log(K^+ + K^-) 166 , st_bet_pos :: !BetState -- ^ bettor state, pos 167 , st_bet_neg :: !BetState -- ^ bettor state, neg 168 } 169 170 -- construction --------------------------------------------------------------- 171 172 -- | Build a 'Config' for the bounded-mean test. 173 -- 174 -- Each per-direction safe-bet ceiling @lambda_max@ is set so that 175 -- the wealth factor stays nonnegative for every admissible 176 -- observation: 177 -- 178 -- * The positive-direction factor is @1 + lambda_p * (x - m)@. 179 -- Since @x@ can dip to @lo@, @x - m@ can reach @lo - m@ (the 180 -- most negative value), so we need 181 -- @lambda_p <= 1 \/ (m - lo)@. The ceiling stored is half this 182 -- to leave numerical margin -- the WSR safety recommendation. 183 -- 184 -- * The negative-direction factor is @1 - lambda_n * (x - m)@. 185 -- Since @x@ can rise to @hi@, @x - m@ can reach @hi - m@, so we 186 -- need @lambda_n <= 1 \/ (hi - m)@; again the ceiling is set to 187 -- half this. 188 -- 189 -- The log-wealth rejection threshold is precomputed as 190 -- @log(2 \/ alpha)@; the 2 reflects that the convex-hedge test 191 -- monitors the sum @K^+_t + K^-_t@, whose initial value is @2@ 192 -- (each side starts at @K = 1@). 193 -- 194 -- Returns 'Left' with a 'ConfigError' on inputs that would leave 195 -- the mathematical regime: any of @m@, @lo@, @hi@, @alpha@ 196 -- non-finite (NaN or infinite); @alpha@ outside @(0, 1)@; 197 -- @lo >= hi@; or @m@ outside the open interval @(lo, hi)@ 198 -- (strict, to avoid the safe-bet ceilings dividing by zero). 199 -- 200 -- >>> let Right cfg = config 0.5 0.0 1.0 1.0e-3 Newton 201 config 202 :: Double -- ^ null mean @m@ 203 -> Double -- ^ sample lower bound @lo@ 204 -> Double -- ^ sample upper bound @hi@ 205 -> Double -- ^ significance level @alpha@ 206 -> Bettor -- ^ bettor strategy 207 -> Either ConfigError Config 208 config !m !lo !hi !alpha !b 209 | not (finite alpha && alpha > 0 && alpha < 1) = 210 Left (InvalidAlpha alpha) 211 | not (finite lo && finite hi && lo < hi) = 212 Left (InvalidBounds lo hi) 213 | not (finite m && lo < m && m < hi) = 214 Left (InvalidNullMean m lo hi) 215 | otherwise = Right (mk_config m m lo hi alpha b) 216 {-# INLINE config #-} 217 218 -- | Build a 'Config' for the interval-null bounded-mean test 219 -- 220 -- @H_0: m_lo <= E[x_t | F_{t-1}] <= m_hi for all t@ 221 -- 222 -- The positive-direction process centres its observations at 223 -- @m_hi@ (betting against @E[x] > m_hi@) and the 224 -- negative-direction process at @m_lo@ (betting against 225 -- @E[x] < m_lo@). Each side's wealth factor has conditional 226 -- expectation at most @1@ under every mean in the interval, so 227 -- the convex-hedge combination and its @log(2 \/ alpha)@ 228 -- threshold apply exactly as under a point null. The safe-bet 229 -- ceilings generalise accordingly: @lambda_p <= 1 \/ (m_hi - lo)@ 230 -- and @lambda_n <= 1 \/ (hi - m_lo)@, each halved for numerical 231 -- margin as in 'config'. 232 -- 233 -- A conditional mean strictly inside the interval makes both 234 -- sides strict supermartingales, so the test is conservative in 235 -- the interior and tightest at the endpoints — where the type-I 236 -- guarantee is still @alpha@. 237 -- 238 -- @configInterval m m@ is equivalent to @config m@; passing 239 -- @m_lo > m_hi@, or endpoints outside the open sample interval, 240 -- returns 'Left' 'InvalidNullInterval'. Other inputs are 241 -- validated as in 'config'. 242 -- 243 -- >>> let Right cfg = configInterval 0.45 0.55 0.0 1.0 1.0e-3 Newton 244 configInterval 245 :: Double -- ^ lower null mean @m_lo@ 246 -> Double -- ^ upper null mean @m_hi@ 247 -> Double -- ^ sample lower bound @lo@ 248 -> Double -- ^ sample upper bound @hi@ 249 -> Double -- ^ significance level @alpha@ 250 -> Bettor -- ^ bettor strategy 251 -> Either ConfigError Config 252 configInterval !mlo !mhi !lo !hi !alpha !b 253 | not (finite alpha && alpha > 0 && alpha < 1) = 254 Left (InvalidAlpha alpha) 255 | not (finite lo && finite hi && lo < hi) = 256 Left (InvalidBounds lo hi) 257 | not (finite mlo && finite mhi && lo < mlo && mlo <= mhi && mhi < hi) = 258 Left (InvalidNullInterval mlo mhi lo hi) 259 | otherwise = Right (mk_config mlo mhi lo hi alpha b) 260 {-# INLINE configInterval #-} 261 262 -- shared constructor behind 'config' and 'configInterval'; assumes 263 -- its arguments already validated. 264 mk_config :: Double -> Double -> Double -> Double -> Double -> Bettor -> Config 265 mk_config !mlo !mhi !lo !hi !alpha !b = Config { 266 cfg_bettor = b 267 , cfg_lam_max_pos = 0.5 / (mhi - lo) 268 , cfg_lam_max_neg = 0.5 / (hi - mlo) 269 , cfg_null_lo = mlo 270 , cfg_null_hi = mhi 271 , cfg_alpha = alpha 272 , cfg_log_thresh = log (2 / alpha) 273 } 274 {-# INLINE mk_config #-} 275 276 -- | The initial 'State' for a fresh streaming test. 277 -- 278 -- Both per-direction log-wealths start at @0@ (i.e., @K = 1@); 279 -- the sup-log-sum starts at @log 2@ (since @K^+_0 + K^-_0 = 2@); 280 -- both bettors start in the per-strategy initial state 281 -- appropriate for the 'Bettor' chosen in the 'Config'. 282 -- 283 -- >>> let s0 = initial cfg 284 initial :: Config -> State 285 initial Config{..} = 286 let !s0 = init_bet cfg_bettor 287 in State { 288 st_n = 0 289 , st_log_w_pos = 0 290 , st_log_w_neg = 0 291 , st_sup_log_sum = log2_dbl 292 , st_bet_pos = s0 293 , st_bet_neg = s0 294 } 295 {-# INLINE initial #-} 296 297 -- streaming ------------------------------------------------------------------ 298 299 -- | Fold one observation into the running 'State'. 300 -- 301 -- Computes the per-direction centred observations 302 -- @z_p = x - m_hi@ and @z_n = x - m_lo@ (identical under a point 303 -- null, where @m_lo = m_hi = m@), queries the two directional 304 -- bettors for their predictable bets, accumulates per-direction 305 -- log-wealth via 306 -- 307 -- @log_w' = log_w + log (1 + lambda * z)@ 308 -- 309 -- (with the symmetric @-lambda@ for the negative direction), then 310 -- updates the running supremum of @log(K^+ + K^-)@ via 311 -- log-sum-exp and steps the bettor states given the newly 312 -- observed centred values. 313 -- 314 -- /Precondition/: @x@ must lie in the @[lo, hi]@ interval given 315 -- to 'config' \/ 'configInterval'. The type-I error guarantee of 316 -- the test depends on this. Out-of-range observations can drive 317 -- the wealth factor negative, taking the construction out of the 318 -- supermartingale regime entirely; the function does not check 319 -- for this. 320 -- 321 -- >>> let s1 = update cfg s0 0.7 322 update :: Config -> State -> Double -> State 323 update Config{..} State{..} !x = 324 let !zp = x - cfg_null_hi 325 !zn = x - cfg_null_lo 326 !lam_p = bet_lambda cfg_bettor cfg_lam_max_pos st_bet_pos 327 !lam_n = bet_lambda cfg_bettor cfg_lam_max_neg st_bet_neg 328 !logw_p = st_log_w_pos + log1p (lam_p * zp) 329 !logw_n = st_log_w_neg + log1p (negate lam_n * zn) 330 -- Skip 'log_sum_exp' when the cheap upper bound 331 -- log_sum_exp a b <= max a b + log 2 332 -- already sits at or below the running max: no update can 333 -- move it. Under H_0 (calibration) this is the common case. 334 !cheap_ub = max logw_p logw_n + log2_dbl 335 !sup_sum 336 | cheap_ub <= st_sup_log_sum = st_sup_log_sum 337 | otherwise = 338 max st_sup_log_sum (log_sum_exp logw_p logw_n) 339 !sp = step_bet cfg_bettor cfg_lam_max_pos st_bet_pos zp 340 !sn = step_bet cfg_bettor cfg_lam_max_neg st_bet_neg (negate zn) 341 in State (st_n + 1) logw_p logw_n sup_sum sp sn 342 {-# INLINE update #-} 343 344 -- | Compute the current 'Verdict' from the running 'State'. 345 -- 346 -- 'Reject' iff the supremum-so-far of @log(K^+_t + K^-_t)@ has 347 -- ever crossed the threshold @log(2 \/ alpha)@ — equivalently, 348 -- the convex-hedge e-process @(K^+ + K^-) \/ 2@ has exceeded 349 -- @1 \/ alpha@ at some point in the stream so far. Under @H_0@, 350 -- by Ville's inequality, the probability of this ever happening 351 -- is at most @alpha@ -- and crucially this bound holds at /every/ 352 -- sample size simultaneously, so the user is free to peek at the 353 -- verdict as often as they like and stop on the first 'Reject'. 354 -- 355 -- >>> decide cfg s0 356 -- Continue 357 decide :: Config -> State -> Verdict 358 decide Config{..} State{..} 359 | st_sup_log_sum >= cfg_log_thresh = Reject 360 | otherwise = Continue 361 {-# INLINE decide #-} 362 363 -- inspection ----------------------------------------------------------------- 364 365 -- | The current @log(K^+_t + K^-_t)@ -- the running log-wealth of 366 -- the convex-hedge combination at the present sample count. 367 -- 368 -- Unlike 'log_wealth_sup' this is not monotone: adverse 369 -- observations decrease it. It is bounded above by 370 -- 'log_wealth_sup', which is what 'decide' tests against the 371 -- rejection threshold. 372 -- 373 -- Starts at @log 2@ (since @K^+_0 + K^-_0 = 2@). 374 -- 375 -- >>> log_wealth s0 376 -- 0.6931471805599453 377 log_wealth :: State -> Double 378 log_wealth State{..} = log_sum_exp st_log_w_pos st_log_w_neg 379 {-# INLINE log_wealth #-} 380 381 -- | The supremum-so-far of @log(K^+_t + K^-_t)@, taken across all 382 -- sample counts up to the current one. This is the test statistic 383 -- the convex-hedge construction actually monitors: it is monotone 384 -- nondecreasing in the sample count, and 'decide' rejects exactly 385 -- when it crosses @log(2 \/ alpha)@. 386 -- 387 -- Starts at @log 2@ (since @K^+_0 + K^-_0 = 2@). 388 -- 389 -- >>> log_wealth_sup s0 390 -- 0.6931471805599453 391 log_wealth_sup :: State -> Double 392 log_wealth_sup State{..} = st_sup_log_sum 393 {-# INLINE log_wealth_sup #-} 394 395 -- | The current log e-value of the convex-hedge e-process: the log 396 -- of @(K^+_t + K^-_t) \/ 2@, i.e. 'log_wealth' minus @log 2@. 397 -- 398 -- Unlike 'log_wealth', this is normalized so that a fresh state 399 -- sits at @0@ (e-value @1@): it is directly comparable across 400 -- test modules regardless of their internal hedging, and is the 401 -- form to use when convex-combining several e-processes. Not 402 -- monotone; bounded above by 'log_evalue_sup'. 403 -- 404 -- >>> log_evalue s0 405 -- 0.0 406 log_evalue :: State -> Double 407 log_evalue s = log_wealth s - log2_dbl 408 {-# INLINE log_evalue #-} 409 410 -- | The supremum-so-far of the log e-value: 'log_wealth_sup' minus 411 -- @log 2@. Monotone nondecreasing, starting at @0@; 'decide' 412 -- rejects exactly when it crosses @log(1 \/ alpha)@. 413 -- 414 -- >>> log_evalue_sup s0 415 -- 0.0 416 log_evalue_sup :: State -> Double 417 log_evalue_sup s = log_wealth_sup s - log2_dbl 418 {-# INLINE log_evalue_sup #-} 419 420 -- | The anytime-valid p-value: the reciprocal of the largest 421 -- e-value attained so far, @min 1 (exp (negate (log_evalue_sup 422 -- s)))@. 423 -- 424 -- Monotone nonincreasing in the sample count, and valid under 425 -- optional stopping: under @H_0@, 426 -- @P(exists t: p_t <= alpha) <= alpha@ for every @alpha@ 427 -- simultaneously. 'decide' returns 'Reject' exactly when this 428 -- value has reached the configured @alpha@ or below. 429 -- 430 -- >>> p_value s0 431 -- 1.0 432 p_value :: State -> Double 433 p_value s = min 1 (exp (negate (log_evalue_sup s))) 434 {-# INLINE p_value #-} 435 436 -- | The number of samples consumed so far. 437 -- 438 -- >>> samples s0 439 -- 0 440 samples :: State -> Int 441 samples = st_n 442 {-# INLINE samples #-}