Main.hs (45925B)
1 {-# LANGUAGE BangPatterns #-} 2 3 module Main where 4 5 import Data.Bits 6 import qualified Data.List as L 7 import Data.Word 8 import qualified Numeric.Eproc.Bernoulli as Bern 9 import qualified Numeric.Eproc.Bernoulli.TwoSided as BernTS 10 import qualified Numeric.Eproc.Bounded as Bounded 11 import qualified Numeric.Eproc.Common as C 12 import qualified Numeric.Eproc.ConfSeq as CS 13 import qualified Numeric.Eproc.Mixture as Mix 14 import qualified Numeric.Eproc.Paired as P 15 import Test.Tasty 16 import Test.Tasty.HUnit 17 import qualified Test.Tasty.QuickCheck as QC 18 19 main :: IO () 20 main = defaultMain $ testGroup "ppad-eproc" [ 21 sanity_tests 22 , calibration_tests 23 , power_tests 24 , interval_tests 25 , two_sample_tests 26 , bernoulli_tests 27 , bettor_smoke_tests 28 , latched_rejection_tests 29 , config_validation_tests 30 , safety_property_tests 31 , two_sided_bernoulli_tests 32 , evalue_accessor_tests 33 , mixture_tests 34 , confseq_tests 35 ] 36 37 -- partial helper: tests below hardcode valid configs. 38 ok :: Either e a -> a 39 ok (Right x) = x 40 ok (Left _) = error "test: invalid config" 41 42 -- prng ----------------------------------------------------------------------- 43 44 -- inline PCG-style PRNG, no external deps. 45 46 newtype Gen = Gen Word64 47 48 mk_gen :: Word64 -> Gen 49 mk_gen = Gen 50 51 step_gen :: Gen -> (Word64, Gen) 52 step_gen (Gen s) = 53 let !s' = s * 6364136223846793005 + 1442695040888963407 54 in (s', Gen s') 55 56 next_double :: Gen -> (Double, Gen) 57 next_double g = 58 let (w, g') = step_gen g 59 !x = fromIntegral (w `shiftR` 11 .&. 0x1FFFFFFFFFFFFF) / 60 9007199254740992 61 in (x, g') 62 63 bernoulli :: Double -> Gen -> (Double, Gen) 64 bernoulli !p g = 65 let (u, g') = next_double g 66 in (if u < p then 1.0 else 0.0, g') 67 68 -- per-trial independent seeds via a splitmix-style finalizer. 69 -- previously this just stepped the prng once per trial, which made 70 -- consecutive trials share all but one observation -- fine under a 71 -- symmetric H_0 (rare streaks cancel), catastrophic under a skewed 72 -- one (rare streaks dominate all overlapping trials). 73 gen_seq :: Gen -> [Gen] 74 gen_seq (Gen s0) = 75 [Gen (mix64 (s0 + fromIntegral i)) | i <- [(0 :: Word64) ..]] 76 where 77 mix64 x = 78 let !y = (x `xor` (x `shiftR` 30)) * 0xbf58476d1ce4e5b9 79 !z = (y `xor` (y `shiftR` 27)) * 0x94d049bb133111eb 80 in z `xor` (z `shiftR` 31) 81 82 -- harness -------------------------------------------------------------------- 83 84 -- run a sequential mean test on a stream of n bernoulli(p) samples, 85 -- with the early-stopping rule built in. returns (verdict, samples 86 -- consumed). 87 run_bounded_bernoulli 88 :: Bounded.Config 89 -> Double -- ^ p 90 -> Int -- ^ budget 91 -> Gen 92 -> (Bounded.Verdict, Int) 93 run_bounded_bernoulli cfg p budget g0 = go 0 g0 (Bounded.initial cfg) 94 where 95 go !n !g !st 96 | n >= budget = (Bounded.decide cfg st, n) 97 | otherwise = case Bounded.decide cfg st of 98 Bounded.Reject -> (Bounded.Reject, n) 99 Bounded.Continue -> 100 let (x, g') = bernoulli p g 101 st' = Bounded.update cfg st x 102 in go (n + 1) g' st' 103 104 -- fraction of trials that rejected. 105 rejection_rate 106 :: Bounded.Config 107 -> Double -- ^ true bernoulli p 108 -> Int -- ^ budget per trial 109 -> Int -- ^ number of trials 110 -> Word64 -- ^ seed 111 -> Double 112 rejection_rate cfg p budget trials seed = 113 let gens = take trials (gen_seq (mk_gen seed)) 114 rejects = length 115 [ () | g <- gens 116 , let (v, _) = run_bounded_bernoulli cfg p budget g 117 , v == Bounded.Reject ] 118 in fromIntegral rejects / fromIntegral trials 119 120 run_paired 121 :: P.Config 122 -> Double 123 -> Double -- ^ p for A and B 124 -> Int 125 -> Gen 126 -> (P.Verdict, Int) 127 run_paired cfg pa pb budget g0 = go 0 g0 (P.initial cfg) 128 where 129 go !n !g !st 130 | n >= budget = (P.decide cfg st, n) 131 | otherwise = case P.decide cfg st of 132 P.Reject -> (P.Reject, n) 133 P.Continue -> 134 let (a, g1) = bernoulli pa g 135 (b, g2) = bernoulli pb g1 136 st' = P.update cfg st (a, b) 137 in go (n + 1) g2 st' 138 139 paired_avg_rate 140 :: P.Config 141 -> Double 142 -> Double 143 -> Int 144 -> Int 145 -> Word64 146 -> Double 147 paired_avg_rate cfg pa pb budget trials seed = 148 let gens = take trials (gen_seq (mk_gen seed)) 149 rejects = length 150 [ () | g <- gens 151 , let (v, _) = run_paired cfg pa pb budget g 152 , v == P.Reject ] 153 in fromIntegral rejects / fromIntegral trials 154 155 -- sanity --------------------------------------------------------------------- 156 157 -- with all-zero deviations from the null mean, no rejection. 158 sanity_tests :: TestTree 159 sanity_tests = testGroup "sanity" [ 160 testCase "degenerate input never rejects" $ do 161 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-6 Bounded.Newton) 162 xs = replicate 5000 0.5 163 st = L.foldl' (Bounded.update cfg) (Bounded.initial cfg) xs 164 Bounded.decide cfg st @?= Bounded.Continue 165 , testCase "two-sided thresholds applied symmetrically" $ do 166 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-6 Bounded.Newton) 167 Bounded.decide cfg (Bounded.initial cfg) @?= Bounded.Continue 168 ] 169 170 -- null calibration ----------------------------------------------------------- 171 172 -- under H_0, with optional stopping, the empirical rejection rate should be 173 -- bounded by alpha. ville's inequality is typically conservative on bernoulli, 174 -- so the slack is small. 175 calibration_tests :: TestTree 176 calibration_tests = testGroup "null calibration" [ 177 testCase "Newton, Bernoulli(0.5), m=0.5, alpha=0.05" $ do 178 let cfg = ok (Bounded.config 0.5 0.0 1.0 0.05 Bounded.Newton) 179 rate = rejection_rate cfg 0.5 2000 200 12345 180 -- expected rate <= 0.05; allow up to ~0.08 slack for sampling 181 -- variability over 200 trials (sigma ~ 0.015). 182 assertBool ("FPR " ++ show rate ++ " exceeded slack") $ 183 rate <= 0.08 184 , testCase "Adaptive, Bernoulli(0.5), m=0.5, alpha=0.05" $ do 185 let cfg = ok (Bounded.config 0.5 0.0 1.0 0.05 Bounded.Adaptive) 186 rate = rejection_rate cfg 0.5 2000 200 67890 187 assertBool ("FPR " ++ show rate ++ " exceeded slack") $ 188 rate <= 0.08 189 ] 190 191 -- power ---------------------------------------------------------------------- 192 193 -- under a clear shift, all (or nearly all) trials reject within budget. 194 power_tests :: TestTree 195 power_tests = testGroup "power" [ 196 testCase "Newton detects Bernoulli(0.7) vs m=0.5" $ do 197 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton) 198 rate = rejection_rate cfg 0.7 5000 100 11111 199 assertBool ("power " ++ show rate ++ " too low") $ 200 rate >= 0.95 201 , testCase "Adaptive detects Bernoulli(0.7) vs m=0.5" $ do 202 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Adaptive) 203 rate = rejection_rate cfg 0.7 5000 100 22222 204 assertBool ("power " ++ show rate ++ " too low") $ 205 rate >= 0.95 206 ] 207 208 -- interval null -------------------------------------------------------------- 209 210 interval_tests :: TestTree 211 interval_tests = testGroup "interval null" [ 212 QC.testProperty "configInterval m m matches config m exactly" $ 213 QC.forAll arb_bettor $ \b -> 214 QC.forAll (QC.listOf unit_double) $ \xs -> 215 let cfgP = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 b) 216 cfgI = ok (Bounded.configInterval 0.5 0.5 0.0 1.0 1.0e-3 b) 217 stP = L.foldl' (Bounded.update cfgP) (Bounded.initial cfgP) xs 218 stI = L.foldl' (Bounded.update cfgI) (Bounded.initial cfgI) xs 219 in Bounded.log_wealth stP == Bounded.log_wealth stI && 220 Bounded.log_wealth_sup stP == Bounded.log_wealth_sup stI 221 , testCase "interior systematic tolerated where point null rejects" $ do 222 -- a constant stream at 0.52: inside [0.45, 0.55], outside the 223 -- point null at 0.5. The interval test accumulates nothing; 224 -- the point test rejects. 225 let cfgI = ok (Bounded.configInterval 0.45 0.55 0.0 1.0 226 1.0e-6 Bounded.Newton) 227 cfgP = ok (Bounded.config 0.5 0.0 1.0 1.0e-6 Bounded.Newton) 228 xs = replicate 5000 0.52 229 stI = L.foldl' (Bounded.update cfgI) (Bounded.initial cfgI) xs 230 stP = L.foldl' (Bounded.update cfgP) (Bounded.initial cfgP) xs 231 Bounded.decide cfgI stI @?= Bounded.Continue 232 Bounded.decide cfgP stP @?= Bounded.Reject 233 , testCase "calibration at the null boundary (p = m_hi)" $ do 234 -- the guarantee is tightest at the endpoints; the empirical 235 -- rate there should still be bounded by alpha (same slack 236 -- convention as the point-null calibration tests). 237 let cfg = ok (Bounded.configInterval 0.45 0.55 0.0 1.0 238 0.05 Bounded.Newton) 239 rate = rejection_rate cfg 0.55 2000 200 34567 240 assertBool ("FPR " ++ show rate ++ " exceeded slack") $ 241 rate <= 0.08 242 , testCase "calibration in the interior (p = 0.5)" $ do 243 let cfg = ok (Bounded.configInterval 0.45 0.55 0.0 1.0 244 0.05 Bounded.Newton) 245 rate = rejection_rate cfg 0.5 2000 200 45678 246 assertBool ("FPR " ++ show rate ++ " exceeded slack") $ 247 rate <= 0.08 248 , testCase "detects Bernoulli(0.7) above [0.45, 0.55]" $ do 249 let cfg = ok (Bounded.configInterval 0.45 0.55 0.0 1.0 250 1.0e-3 Bounded.Newton) 251 rate = rejection_rate cfg 0.7 5000 100 56789 252 assertBool ("power " ++ show rate ++ " too low") $ 253 rate >= 0.95 254 , testCase "detects Bernoulli(0.3) below [0.45, 0.55]" $ do 255 let cfg = ok (Bounded.configInterval 0.45 0.55 0.0 1.0 256 1.0e-3 Bounded.Newton) 257 rate = rejection_rate cfg 0.3 5000 100 67890 258 assertBool ("power " ++ show rate ++ " too low") $ 259 rate >= 0.95 260 , testCase "m_lo > m_hi rejected as InvalidNullInterval" $ 261 case Bounded.configInterval 0.6 0.4 0.0 1.0 0.01 Bounded.Newton of 262 Left (C.InvalidNullInterval _ _ _ _) -> pure () 263 Left e -> assertFailure ("wrong error: " ++ show e) 264 Right _ -> assertFailure "expected Left" 265 , testCase "m_lo == lo rejected" $ 266 case Bounded.configInterval 0.0 0.5 0.0 1.0 0.01 Bounded.Newton of 267 Left (C.InvalidNullInterval _ _ _ _) -> pure () 268 Left e -> assertFailure ("wrong error: " ++ show e) 269 Right _ -> assertFailure "expected Left" 270 , testCase "m_hi == hi rejected" $ 271 case Bounded.configInterval 0.5 1.0 0.0 1.0 0.01 Bounded.Newton of 272 Left (C.InvalidNullInterval _ _ _ _) -> pure () 273 Left e -> assertFailure ("wrong error: " ++ show e) 274 Right _ -> assertFailure "expected Left" 275 , testCase "NaN endpoint rejected" $ 276 case Bounded.configInterval (0 / 0) 0.5 0.0 1.0 0.01 Bounded.Newton of 277 Left (C.InvalidNullInterval _ _ _ _) -> pure () 278 Left e -> assertFailure ("wrong error: " ++ show e) 279 Right _ -> assertFailure "expected Left" 280 ] 281 282 -- two-sample paired test ----------------------------------------------------- 283 284 two_sample_tests :: TestTree 285 two_sample_tests = testGroup "two-sample" [ 286 testCase "identical distributions don't reject" $ do 287 let cfg = ok (P.config 0.0 1.0 1.0e-3 Bounded.Newton) 288 rate = paired_avg_rate cfg 0.5 0.5 2000 100 33333 289 assertBool ("FPR " ++ show rate) $ rate <= 0.05 290 , testCase "different distributions reject" $ do 291 let cfg = ok (P.config 0.0 1.0 1.0e-3 Bounded.Newton) 292 rate = paired_avg_rate cfg 0.3 0.7 5000 100 44444 293 assertBool ("power " ++ show rate) $ rate >= 0.95 294 ] 295 296 -- bernoulli (one-sided rate) ------------------------------------------------- 297 298 run_bernoulli 299 :: Bern.Config 300 -> Double -- ^ true rate p 301 -> Int -- ^ budget 302 -> Gen 303 -> (Bern.Verdict, Int) 304 run_bernoulli cfg p budget g0 = go 0 g0 (Bern.initial cfg) 305 where 306 go !n !g !st 307 | n >= budget = (Bern.decide cfg st, n) 308 | otherwise = case Bern.decide cfg st of 309 Bern.Reject -> (Bern.Reject, n) 310 Bern.Continue -> 311 let (u, g') = next_double g 312 !x = u < p 313 st' = Bern.update cfg st x 314 in go (n + 1) g' st' 315 316 bernoulli_rate 317 :: Bern.Config 318 -> Double -- ^ true rate p 319 -> Int -- ^ budget per trial 320 -> Int -- ^ number of trials 321 -> Word64 -- ^ seed 322 -> Double 323 bernoulli_rate cfg p budget trials seed = 324 let gens = take trials (gen_seq (mk_gen seed)) 325 rejects = length 326 [ () | g <- gens 327 , let (v, _) = run_bernoulli cfg p budget g 328 , v == Bern.Reject ] 329 in fromIntegral rejects / fromIntegral trials 330 331 bernoulli_tests :: TestTree 332 bernoulli_tests = testGroup "bernoulli" [ 333 testCase "all-zero stream never rejects" $ do 334 let cfg = ok (Bern.config 0.05 1.0e-6 Bern.Newton) 335 xs = replicate 5000 False 336 st = L.foldl' (Bern.update cfg) (Bern.initial cfg) xs 337 Bern.decide cfg st @?= Bern.Continue 338 , testCase "Newton FPR under H_0 (p = p_0 = 0.05)" $ do 339 let cfg = ok (Bern.config 0.05 0.05 Bern.Newton) 340 rate = bernoulli_rate cfg 0.05 2000 200 55555 341 assertBool ("FPR " ++ show rate ++ " exceeded slack") $ 342 rate <= 0.08 343 , testCase "Adaptive FPR under H_0 (p = p_0 = 0.05)" $ do 344 let cfg = ok (Bern.config 0.05 0.05 Bern.Adaptive) 345 rate = bernoulli_rate cfg 0.05 2000 200 66666 346 assertBool ("FPR " ++ show rate ++ " exceeded slack") $ 347 rate <= 0.08 348 , testCase "Newton detects p = 0.3 vs p_0 = 0.05" $ do 349 let cfg = ok (Bern.config 0.05 1.0e-3 Bern.Newton) 350 rate = bernoulli_rate cfg 0.3 5000 100 77777 351 assertBool ("power " ++ show rate ++ " too low") $ 352 rate >= 0.95 353 , testCase "Adaptive detects p = 0.3 vs p_0 = 0.05" $ do 354 let cfg = ok (Bern.config 0.05 1.0e-3 Bern.Adaptive) 355 rate = bernoulli_rate cfg 0.3 5000 100 88888 356 assertBool ("power " ++ show rate ++ " too low") $ 357 rate >= 0.95 358 ] 359 360 -- bettor smoke tests --------------------------------------------------------- 361 362 -- each bettor produces a well-defined state and decision when run on a small 363 -- deterministic stream. 364 bettor_smoke_tests :: TestTree 365 bettor_smoke_tests = testGroup "bettor smoke" [ 366 testCase "fixed bettor runs without error (bounded)" $ do 367 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 (Bounded.Fixed 0.5)) 368 xs = take 100 (cycle [0.0, 1.0]) 369 st = L.foldl' (Bounded.update cfg) (Bounded.initial cfg) xs 370 assertBool "samples advanced" (Bounded.samples st == 100) 371 , testCase "Newton bettor runs without error (bounded)" $ do 372 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton) 373 xs = take 100 (cycle [0.0, 1.0]) 374 st = L.foldl' (Bounded.update cfg) (Bounded.initial cfg) xs 375 assertBool "samples advanced" (Bounded.samples st == 100) 376 , testCase "Adaptive bettor runs without error (bounded)" $ do 377 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Adaptive) 378 xs = take 100 (cycle [0.0, 1.0]) 379 st = L.foldl' (Bounded.update cfg) (Bounded.initial cfg) xs 380 assertBool "samples advanced" (Bounded.samples st == 100) 381 , testCase "fixed bettor runs without error (bernoulli)" $ do 382 let cfg = ok (Bern.config 0.5 1.0e-3 (Bern.Fixed 0.5)) 383 xs = take 100 (cycle [True, False]) 384 st = L.foldl' (Bern.update cfg) (Bern.initial cfg) xs 385 assertBool "samples advanced" (Bern.samples st == 100) 386 , testCase "Newton bettor runs without error (bernoulli)" $ do 387 let cfg = ok (Bern.config 0.5 1.0e-3 Bern.Newton) 388 xs = take 100 (cycle [True, False]) 389 st = L.foldl' (Bern.update cfg) (Bern.initial cfg) xs 390 assertBool "samples advanced" (Bern.samples st == 100) 391 , testCase "Adaptive bettor runs without error (bernoulli)" $ do 392 let cfg = ok (Bern.config 0.5 1.0e-3 Bern.Adaptive) 393 xs = take 100 (cycle [True, False]) 394 st = L.foldl' (Bern.update cfg) (Bern.initial cfg) xs 395 assertBool "samples advanced" (Bern.samples st == 100) 396 ] 397 398 -- latched rejection ---------------------------------------------------------- 399 400 -- once the wealth crosses threshold, subsequent observations driving the 401 -- current wealth back below threshold must not unrejection the test. 402 latched_rejection_tests :: TestTree 403 latched_rejection_tests = testGroup "latched rejection" [ 404 testCase "bounded: cross then drown stays rejected" $ do 405 -- alpha = 0.5 => threshold log(2/0.5) = log 4 ~ 1.386. 406 -- Fixed 1.0 with x=1 grows log_w_pos by log 1.5 ~ 0.405/step; 407 -- five 1s push it past threshold. Then forty 0s drop it well 408 -- below. 409 let cfg = ok (Bounded.config 0.5 0.0 1.0 0.5 (Bounded.Fixed 1.0)) 410 xs1 = replicate 5 1.0 411 xs2 = replicate 40 0.0 412 st1 = L.foldl' (Bounded.update cfg) (Bounded.initial cfg) xs1 413 st2 = L.foldl' (Bounded.update cfg) st1 xs2 414 Bounded.decide cfg st1 @?= Bounded.Reject 415 Bounded.decide cfg st2 @?= Bounded.Reject 416 , testCase "bernoulli: cross then drown stays rejected" $ do 417 let cfg = ok (Bern.config 0.05 0.5 (Bern.Fixed 1.0)) 418 xs1 = replicate 5 True 419 xs2 = replicate 200 False 420 st1 = L.foldl' (Bern.update cfg) (Bern.initial cfg) xs1 421 st2 = L.foldl' (Bern.update cfg) st1 xs2 422 Bern.decide cfg st1 @?= Bern.Reject 423 Bern.decide cfg st2 @?= Bern.Reject 424 ] 425 426 -- config validation ---------------------------------------------------------- 427 428 config_validation_tests :: TestTree 429 config_validation_tests = testGroup "config validation" [ 430 testCase "Bounded: alpha <= 0 rejected" $ 431 assertLeft (Bounded.config 0.5 0.0 1.0 0.0 Bounded.Newton) 432 , testCase "Bounded: alpha >= 1 rejected" $ 433 assertLeft (Bounded.config 0.5 0.0 1.0 1.5 Bounded.Newton) 434 , testCase "Bounded: lo >= hi rejected" $ 435 assertLeft (Bounded.config 0.5 1.0 0.0 0.01 Bounded.Newton) 436 , testCase "Bounded: m == lo rejected" $ 437 assertLeft (Bounded.config 0.0 0.0 1.0 0.01 Bounded.Newton) 438 , testCase "Bounded: m == hi rejected" $ 439 assertLeft (Bounded.config 1.0 0.0 1.0 0.01 Bounded.Newton) 440 , testCase "Bounded: m outside [lo, hi] rejected" $ 441 assertLeft (Bounded.config 2.0 0.0 1.0 0.01 Bounded.Newton) 442 , testCase "Bernoulli: alpha <= 0 rejected" $ 443 assertLeft (Bern.config 0.5 0.0 Bern.Newton) 444 , testCase "Bernoulli: alpha >= 1 rejected" $ 445 assertLeft (Bern.config 0.5 1.0 Bern.Newton) 446 , testCase "Bernoulli: p0 == 0 rejected" $ 447 assertLeft (Bern.config 0.0 0.05 Bern.Newton) 448 , testCase "Bernoulli: p0 == 1 rejected" $ 449 assertLeft (Bern.config 1.0 0.05 Bern.Newton) 450 , testCase "Paired: alpha out of range rejected" $ 451 assertLeft (P.config 0.0 1.0 0.0 Bounded.Newton) 452 , testCase "Paired: lo >= hi rejected" $ 453 assertLeft (P.config 1.0 0.0 0.01 Bounded.Newton) 454 , testCase "Bounded: infinite bounds rejected" $ 455 assertLeft (Bounded.config 0.0 nInf pInf 0.01 Bounded.Newton) 456 , testCase "Bounded: NaN m rejected" $ 457 assertLeft (Bounded.config nan 0.0 1.0 0.01 Bounded.Newton) 458 , testCase "Bounded: NaN alpha rejected" $ 459 assertLeft (Bounded.config 0.5 0.0 1.0 nan Bounded.Newton) 460 , testCase "Bernoulli: NaN p0 rejected" $ 461 assertLeft (Bern.config nan 0.01 Bern.Newton) 462 , testCase "Bernoulli: infinite alpha rejected" $ 463 assertLeft (Bern.config 0.05 pInf Bern.Newton) 464 , testCase "Paired: infinite hi rejected" $ 465 assertLeft (P.config 0.0 pInf 0.01 Bounded.Newton) 466 ] 467 where 468 nan, pInf, nInf :: Double 469 nan = 0 / 0 470 pInf = 1 / 0 471 nInf = negate (1 / 0) 472 assertLeft :: Either C.ConfigError a -> Assertion 473 assertLeft e = case e of 474 Left _ -> pure () 475 Right _ -> assertFailure "expected Left" 476 477 -- two-sided bernoulli -------------------------------------------------------- 478 479 run_ts_bernoulli 480 :: BernTS.Config 481 -> Double -- ^ true rate p 482 -> Int -- ^ budget 483 -> Gen 484 -> (BernTS.Verdict, Int) 485 run_ts_bernoulli cfg p budget g0 = 486 go 0 g0 (BernTS.initial cfg) 487 where 488 go !n !g !st 489 | n >= budget = (BernTS.decide cfg st, n) 490 | otherwise = case BernTS.decide cfg st of 491 BernTS.Reject -> (BernTS.Reject, n) 492 BernTS.Continue -> 493 let (u, g') = next_double g 494 !x = u < p 495 st' = BernTS.update cfg st x 496 in go (n + 1) g' st' 497 498 ts_bernoulli_rate 499 :: BernTS.Config 500 -> Double 501 -> Int 502 -> Int 503 -> Word64 504 -> Double 505 ts_bernoulli_rate cfg p budget trials seed = 506 let gens = take trials (gen_seq (mk_gen seed)) 507 rejects = length 508 [ () | g <- gens 509 , let (v, _) = run_ts_bernoulli cfg p budget g 510 , v == BernTS.Reject ] 511 in fromIntegral rejects / fromIntegral trials 512 513 two_sided_bernoulli_tests :: TestTree 514 two_sided_bernoulli_tests = testGroup "two-sided bernoulli" [ 515 testCase "constant at p_0 doesn't reject" $ do 516 -- Bernoulli(0.5) with p_0 = 0.5 is under the null. 517 let cfg = ok (BernTS.config 0.5 1.0e-6 BernTS.Newton) 518 -- alternating True/False keeps the empirical rate at 0.5. 519 xs = take 5000 (cycle [True, False]) 520 st = L.foldl' (BernTS.update cfg) (BernTS.initial cfg) xs 521 BernTS.decide cfg st @?= BernTS.Continue 522 , testCase "detects upward shift (p = 0.7 vs p_0 = 0.5)" $ do 523 let cfg = ok (BernTS.config 0.5 1.0e-3 BernTS.Newton) 524 rate = ts_bernoulli_rate cfg 0.7 5000 100 111222 525 assertBool ("power " ++ show rate ++ " too low") $ 526 rate >= 0.95 527 , testCase "detects downward shift (p = 0.3 vs p_0 = 0.5)" $ do 528 let cfg = ok (BernTS.config 0.5 1.0e-3 BernTS.Newton) 529 rate = ts_bernoulli_rate cfg 0.3 5000 100 333444 530 assertBool ("power " ++ show rate ++ " too low") $ 531 rate >= 0.95 532 , testCase "Adaptive detects shift (p = 0.7 vs p_0 = 0.5)" $ do 533 let cfg = ok (BernTS.config 0.5 1.0e-3 BernTS.Adaptive) 534 rate = ts_bernoulli_rate cfg 0.7 5000 100 777888 535 assertBool ("power " ++ show rate ++ " too low") $ 536 rate >= 0.95 537 , testCase "FPR at p = p_0 = 0.5 within slack" $ do 538 let cfg = ok (BernTS.config 0.5 0.05 BernTS.Newton) 539 rate = ts_bernoulli_rate cfg 0.5 2000 200 555666 540 assertBool ("FPR " ++ show rate ++ " exceeded slack") $ 541 rate <= 0.08 542 , testCase "latched: cross then drown stays rejected" $ do 543 let cfg = ok (BernTS.config 0.5 0.5 (BernTS.Fixed 1.0)) 544 -- ten 1s push the positive side well past threshold. 545 xs1 = replicate 10 True 546 -- then two hundred 0s drop the current wealth, but the 547 -- latch must hold. 548 xs2 = replicate 200 False 549 st1 = L.foldl' (BernTS.update cfg) (BernTS.initial cfg) xs1 550 st2 = L.foldl' (BernTS.update cfg) st1 xs2 551 BernTS.decide cfg st1 @?= BernTS.Reject 552 BernTS.decide cfg st2 @?= BernTS.Reject 553 , testCase "config: NaN p0 rejected" $ do 554 let nan = 0/0 :: Double 555 case BernTS.config nan 0.05 BernTS.Newton of 556 Left _ -> pure () 557 Right _ -> assertFailure "expected Left" 558 , testCase "config: alpha out of range rejected" $ 559 case BernTS.config 0.5 1.5 BernTS.Newton of 560 Left _ -> pure () 561 Right _ -> assertFailure "expected Left" 562 ] 563 564 -- safety properties ---------------------------------------------------------- 565 566 unit_double :: QC.Gen Double 567 unit_double = QC.choose (0, 1) 568 569 arb_bettor :: QC.Gen C.Bettor 570 arb_bettor = QC.oneof [ 571 pure C.Adaptive 572 , pure C.Newton 573 , C.Fixed <$> QC.choose (-10, 10) -- intentionally include unsafe values 574 ] 575 576 finite :: Double -> Bool 577 finite x = not (isNaN x) && not (isInfinite x) 578 579 monotone_reject_bounded :: [Bounded.Verdict] -> Bool 580 monotone_reject_bounded [] = True 581 monotone_reject_bounded (Bounded.Continue : rest) = monotone_reject_bounded rest 582 monotone_reject_bounded (Bounded.Reject : rest) = all (== Bounded.Reject) rest 583 584 monotone_reject_bern :: [Bern.Verdict] -> Bool 585 monotone_reject_bern [] = True 586 monotone_reject_bern (Bern.Continue : rest) = monotone_reject_bern rest 587 monotone_reject_bern (Bern.Reject : rest) = all (== Bern.Reject) rest 588 589 monotone_reject_bern_ts :: [BernTS.Verdict] -> Bool 590 monotone_reject_bern_ts [] = True 591 monotone_reject_bern_ts (BernTS.Continue : rest) = monotone_reject_bern_ts rest 592 monotone_reject_bern_ts (BernTS.Reject : rest) = all (== BernTS.Reject) rest 593 594 safety_property_tests :: TestTree 595 safety_property_tests = testGroup "safety properties" [ 596 QC.testProperty "Bounded: log_wealth finite after any admissible stream" $ 597 QC.forAll arb_bettor $ \b -> 598 QC.forAll (QC.listOf unit_double) $ \xs -> 599 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 b) 600 st = L.foldl' (Bounded.update cfg) (Bounded.initial cfg) xs 601 in finite (Bounded.log_wealth st) && 602 finite (Bounded.log_wealth_sup st) 603 604 , QC.testProperty "Bernoulli: log_wealth finite after any admissible stream" $ 605 QC.forAll arb_bettor $ \b -> 606 QC.forAll QC.arbitrary $ \xs -> 607 let cfg = ok (Bern.config 0.05 1.0e-3 b) 608 st = L.foldl' (Bern.update cfg) (Bern.initial cfg) (xs :: [Bool]) 609 in finite (Bern.log_wealth st) && finite (Bern.log_wealth_sup st) 610 611 , QC.testProperty "Bounded: Fixed with arbitrary lambda is safe" $ 612 QC.forAll (QC.choose (-1000, 1000)) $ \lam -> 613 QC.forAll (QC.listOf unit_double) $ \xs -> 614 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 (C.Fixed lam)) 615 st = L.foldl' (Bounded.update cfg) (Bounded.initial cfg) xs 616 in finite (Bounded.log_wealth st) 617 618 , QC.testProperty "Bernoulli: Fixed with arbitrary lambda is safe" $ 619 QC.forAll (QC.choose (-1000, 1000)) $ \lam -> 620 QC.forAll QC.arbitrary $ \xs -> 621 let cfg = ok (Bern.config 0.05 1.0e-3 (C.Fixed lam)) 622 st = L.foldl' (Bern.update cfg) (Bern.initial cfg) (xs :: [Bool]) 623 in finite (Bern.log_wealth st) 624 625 , QC.testProperty "Bounded: log_wealth_sup is monotone nondecreasing" $ 626 QC.forAll arb_bettor $ \b -> 627 QC.forAll (QC.listOf unit_double) $ \xs -> 628 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 b) 629 sts = scanl (Bounded.update cfg) (Bounded.initial cfg) xs 630 lws = map Bounded.log_wealth_sup sts 631 in and (zipWith (<=) lws (drop 1 lws)) 632 633 , QC.testProperty "Bernoulli: log_wealth_sup is monotone nondecreasing" $ 634 QC.forAll arb_bettor $ \b -> 635 QC.forAll QC.arbitrary $ \xs -> 636 let cfg = ok (Bern.config 0.05 1.0e-3 b) 637 sts = scanl (Bern.update cfg) (Bern.initial cfg) (xs :: [Bool]) 638 lws = map Bern.log_wealth_sup sts 639 in and (zipWith (<=) lws (drop 1 lws)) 640 641 , QC.testProperty "Bounded: log_wealth bounded above by log_wealth_sup" $ 642 QC.forAll arb_bettor $ \b -> 643 QC.forAll (QC.listOf unit_double) $ \xs -> 644 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 b) 645 sts = scanl (Bounded.update cfg) (Bounded.initial cfg) xs 646 in all (\s -> Bounded.log_wealth s <= Bounded.log_wealth_sup s) sts 647 648 , QC.testProperty "Bernoulli: log_wealth bounded above by log_wealth_sup" $ 649 QC.forAll arb_bettor $ \b -> 650 QC.forAll QC.arbitrary $ \xs -> 651 let cfg = ok (Bern.config 0.05 1.0e-3 b) 652 sts = scanl (Bern.update cfg) (Bern.initial cfg) (xs :: [Bool]) 653 in all (\s -> Bern.log_wealth s <= Bern.log_wealth_sup s) sts 654 655 , QC.testProperty "Bounded: rejection is latched" $ 656 QC.forAll arb_bettor $ \b -> 657 QC.forAll (QC.listOf unit_double) $ \xs -> 658 let cfg = ok (Bounded.config 0.5 0.0 1.0 0.5 b) 659 sts = scanl (Bounded.update cfg) (Bounded.initial cfg) xs 660 vs = map (Bounded.decide cfg) sts 661 in monotone_reject_bounded vs 662 663 , QC.testProperty "Bernoulli: rejection is latched" $ 664 QC.forAll arb_bettor $ \b -> 665 QC.forAll QC.arbitrary $ \xs -> 666 let cfg = ok (Bern.config 0.5 0.5 b) 667 sts = scanl (Bern.update cfg) (Bern.initial cfg) (xs :: [Bool]) 668 vs = map (Bern.decide cfg) sts 669 in monotone_reject_bern vs 670 671 , QC.testProperty "BernTS: log_wealth finite after any admissible stream" $ 672 QC.forAll arb_bettor $ \b -> 673 QC.forAll QC.arbitrary $ \xs -> 674 let cfg = ok (BernTS.config 0.5 1.0e-3 b) 675 st = L.foldl' (BernTS.update cfg) (BernTS.initial cfg) 676 (xs :: [Bool]) 677 in finite (BernTS.log_wealth st) && finite (BernTS.log_wealth_sup st) 678 679 , QC.testProperty "BernTS: Fixed with arbitrary lambda is safe" $ 680 QC.forAll (QC.choose (-1000, 1000)) $ \lam -> 681 QC.forAll QC.arbitrary $ \xs -> 682 let cfg = ok (BernTS.config 0.5 1.0e-3 (C.Fixed lam)) 683 st = L.foldl' (BernTS.update cfg) (BernTS.initial cfg) 684 (xs :: [Bool]) 685 in finite (BernTS.log_wealth st) 686 687 , QC.testProperty "BernTS: log_wealth_sup is monotone nondecreasing" $ 688 QC.forAll arb_bettor $ \b -> 689 QC.forAll QC.arbitrary $ \xs -> 690 let cfg = ok (BernTS.config 0.5 1.0e-3 b) 691 sts = scanl (BernTS.update cfg) (BernTS.initial cfg) (xs :: [Bool]) 692 lws = map BernTS.log_wealth_sup sts 693 in and (zipWith (<=) lws (drop 1 lws)) 694 695 , QC.testProperty "BernTS: log_wealth bounded above by log_wealth_sup" $ 696 QC.forAll arb_bettor $ \b -> 697 QC.forAll QC.arbitrary $ \xs -> 698 let cfg = ok (BernTS.config 0.5 1.0e-3 b) 699 sts = scanl (BernTS.update cfg) (BernTS.initial cfg) (xs :: [Bool]) 700 in all (\s -> BernTS.log_wealth s <= BernTS.log_wealth_sup s) sts 701 702 , QC.testProperty "BernTS: rejection is latched" $ 703 QC.forAll arb_bettor $ \b -> 704 QC.forAll QC.arbitrary $ \xs -> 705 let cfg = ok (BernTS.config 0.5 0.5 b) 706 sts = scanl (BernTS.update cfg) (BernTS.initial cfg) (xs :: [Bool]) 707 vs = map (BernTS.decide cfg) sts 708 in monotone_reject_bern_ts vs 709 ] 710 711 712 unit_pair :: QC.Gen (Double, Double) 713 unit_pair = (,) <$> unit_double <*> unit_double 714 715 evalue_accessor_tests :: TestTree 716 evalue_accessor_tests = testGroup "e-value accessors" [ 717 testCase "fresh states normalize to e-value 1, p-value 1" $ do 718 let bcfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton) 719 ncfg = ok (Bern.config 0.05 1.0e-3 Bern.Newton) 720 tcfg = ok (BernTS.config 0.5 1.0e-3 BernTS.Newton) 721 pcfg = ok (P.config 0.0 1.0 1.0e-3 Bounded.Newton) 722 Bounded.log_evalue (Bounded.initial bcfg) @?= 0.0 723 Bounded.log_evalue_sup (Bounded.initial bcfg) @?= 0.0 724 Bounded.p_value (Bounded.initial bcfg) @?= 1.0 725 Bern.log_evalue (Bern.initial ncfg) @?= 0.0 726 Bern.p_value (Bern.initial ncfg) @?= 1.0 727 BernTS.log_evalue (BernTS.initial tcfg) @?= 0.0 728 BernTS.p_value (BernTS.initial tcfg) @?= 1.0 729 P.log_evalue (P.initial pcfg) @?= 0.0 730 P.p_value (P.initial pcfg) @?= 1.0 731 732 , QC.testProperty "Bounded: log_evalue is log_wealth less log 2" $ 733 QC.forAll arb_bettor $ \b -> 734 QC.forAll (QC.listOf unit_double) $ \xs -> 735 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 b) 736 st = L.foldl' (Bounded.update cfg) (Bounded.initial cfg) xs 737 in Bounded.log_evalue st == Bounded.log_wealth st - C.log2_dbl 738 739 , QC.testProperty "Bernoulli: log_evalue coincides with log_wealth" $ 740 QC.forAll arb_bettor $ \b -> 741 QC.forAll QC.arbitrary $ \xs -> 742 let cfg = ok (Bern.config 0.05 1.0e-3 b) 743 st = L.foldl' (Bern.update cfg) (Bern.initial cfg) (xs :: [Bool]) 744 in Bern.log_evalue st == Bern.log_wealth st 745 746 , QC.testProperty "Bounded: decide agrees with p_value at alpha" $ 747 QC.forAll arb_bettor $ \b -> 748 QC.forAll (QC.listOf unit_double) $ \xs -> 749 let alpha = 0.5 750 cfg = ok (Bounded.config 0.5 0.0 1.0 alpha b) 751 sts = scanl (Bounded.update cfg) (Bounded.initial cfg) xs 752 in all (\s -> (Bounded.decide cfg s == Bounded.Reject) 753 == (Bounded.p_value s <= alpha)) sts 754 755 , QC.testProperty "Bernoulli: decide agrees with p_value at alpha" $ 756 QC.forAll arb_bettor $ \b -> 757 QC.forAll QC.arbitrary $ \xs -> 758 let alpha = 0.5 759 cfg = ok (Bern.config 0.5 alpha b) 760 sts = scanl (Bern.update cfg) (Bern.initial cfg) 761 (xs :: [Bool]) 762 in all (\s -> (Bern.decide cfg s == Bern.Reject) 763 == (Bern.p_value s <= alpha)) sts 764 765 , QC.testProperty "BernTS: decide agrees with p_value at alpha" $ 766 QC.forAll arb_bettor $ \b -> 767 QC.forAll QC.arbitrary $ \xs -> 768 let alpha = 0.5 769 cfg = ok (BernTS.config 0.5 alpha b) 770 sts = scanl (BernTS.update cfg) (BernTS.initial cfg) 771 (xs :: [Bool]) 772 in all (\s -> (BernTS.decide cfg s == BernTS.Reject) 773 == (BernTS.p_value s <= alpha)) sts 774 775 , QC.testProperty "Bounded: p_value monotone nonincreasing" $ 776 QC.forAll arb_bettor $ \b -> 777 QC.forAll (QC.listOf unit_double) $ \xs -> 778 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 b) 779 sts = scanl (Bounded.update cfg) (Bounded.initial cfg) xs 780 ps = map Bounded.p_value sts 781 in and (zipWith (>=) ps (drop 1 ps)) 782 783 , QC.testProperty "Paired: p_value monotone nonincreasing" $ 784 QC.forAll arb_bettor $ \b -> 785 QC.forAll (QC.listOf unit_pair) $ \ps -> 786 let cfg = ok (P.config 0.0 1.0 1.0e-3 b) 787 sts = scanl (P.update cfg) (P.initial cfg) ps 788 pv = map P.p_value sts 789 in and (zipWith (>=) pv (drop 1 pv)) 790 791 , QC.testProperty "Bounded: p_value in [0, 1], evalue below sup" $ 792 QC.forAll arb_bettor $ \b -> 793 QC.forAll (QC.listOf unit_double) $ \xs -> 794 let cfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 b) 795 sts = scanl (Bounded.update cfg) (Bounded.initial cfg) xs 796 in all (\s -> let p = Bounded.p_value s 797 in p >= 0 && p <= 1 && 798 Bounded.log_evalue s 799 <= Bounded.log_evalue_sup s) sts 800 801 , QC.testProperty "Bernoulli: p_value in [0, 1], evalue below sup" $ 802 QC.forAll arb_bettor $ \b -> 803 QC.forAll QC.arbitrary $ \xs -> 804 let cfg = ok (Bern.config 0.05 1.0e-3 b) 805 sts = scanl (Bern.update cfg) (Bern.initial cfg) 806 (xs :: [Bool]) 807 in all (\s -> let p = Bern.p_value s 808 in p >= 0 && p <= 1 && 809 Bern.log_evalue s 810 <= Bern.log_evalue_sup s) sts 811 ] 812 813 -- mixture -------------------------------------------------------------------- 814 815 approx_eq :: Double -> Double -> Bool 816 approx_eq a b = abs (a - b) <= 1.0e-9 * max 1 (max (abs a) (abs b)) 817 818 -- step a censor-style two-component hedge (sign + magnitude) over a 819 -- shared bernoulli stream, feeding the mixture the components' 820 -- current log e-values, with the early-stopping rule built in. 821 run_mixture 822 :: Mix.Config 823 -> BernTS.Config 824 -> Bounded.Config 825 -> Double -- ^ true bernoulli p 826 -> Int -- ^ budget 827 -> Gen 828 -> (Mix.Verdict, Int) 829 run_mixture xc sc mc p budget g0 = 830 go 0 g0 (BernTS.initial sc) (Bounded.initial mc) (Mix.initial xc) 831 where 832 go !n !g !s !m !x 833 | n >= budget = (Mix.decide xc x, n) 834 | otherwise = case Mix.decide xc x of 835 Mix.Reject -> (Mix.Reject, n) 836 Mix.Continue -> 837 let (v, g') = bernoulli p g 838 s' = BernTS.update sc s (v == 1.0) 839 m' = Bounded.update mc m v 840 x' = Mix.update xc x 841 [BernTS.log_evalue s', Bounded.log_evalue m'] 842 in go (n + 1) g' s' m' x' 843 844 mixture_rate :: Double -> Double -> Int -> Int -> Word64 -> Double 845 mixture_rate alpha p budget trials seed = 846 let xc = ok (Mix.config 2 alpha) 847 sc = ok (BernTS.config 0.5 alpha BernTS.Newton) 848 mc = ok (Bounded.config 0.5 0.0 1.0 alpha Bounded.Newton) 849 gens = take trials (gen_seq (mk_gen seed)) 850 rejects = length 851 [ () | g <- gens 852 , let (v, _) = run_mixture xc sc mc p budget g 853 , v == Mix.Reject ] 854 in fromIntegral rejects / fromIntegral trials 855 856 mixture_tests :: TestTree 857 mixture_tests = testGroup "mixture" [ 858 testCase "fresh mixture sits at log K, p-value 1" $ do 859 let cfg = ok (Mix.config 4 1.0e-3) 860 s0 = Mix.initial cfg 861 assertBool "log_wealth is log K" $ 862 approx_eq (Mix.log_wealth s0) (log 4) 863 Mix.log_evalue cfg s0 @?= 0.0 864 Mix.log_evalue_sup cfg s0 @?= 0.0 865 Mix.p_value cfg s0 @?= 1.0 866 Mix.decide cfg s0 @?= Mix.Continue 867 868 , testCase "latch is on the mixture sup, not per-component sups" $ do 869 -- two components peak at different times, each attaining log 870 -- e-value 1.0. A bogus combination of per-component suprema, 871 -- log_sum_exp 1 1 ~ 1.69, crosses the K = 2, alpha = 0.5 872 -- threshold log 4 ~ 1.39; the mixture itself never exceeds 873 -- ~1.003 and must not reject. 874 let cfg = ok (Mix.config 2 0.5) 875 s1 = Mix.update cfg (Mix.initial cfg) [1.0, -5.0] 876 s2 = Mix.update cfg s1 [-5.0, 1.0] 877 Mix.decide cfg s2 @?= Mix.Continue 878 assertBool "mixture sup below threshold" $ 879 Mix.log_wealth_sup s2 < log 4 880 assertBool "per-component-sup combination would cross" $ 881 C.log_sum_exp 1.0 1.0 >= log 4 882 883 , testCase "empty update vector is a no-op" $ do 884 let cfg = ok (Mix.config 2 1.0e-3) 885 s0 = Mix.initial cfg 886 s1 = Mix.update cfg s0 [] 887 Mix.samples s1 @?= 0 888 Mix.log_wealth s1 @?= Mix.log_wealth s0 889 890 , testCase "config validation" $ do 891 let assert_left :: Either C.ConfigError Mix.Config -> Assertion 892 assert_left e = case e of 893 Left _ -> pure () 894 Right _ -> assertFailure "expected Left" 895 assert_left (Mix.config 0 0.05) 896 assert_left (Mix.config (-3) 0.05) 897 assert_left (Mix.config 4 0.0) 898 assert_left (Mix.config 4 1.5) 899 assert_left (Mix.config 4 (0 / 0)) 900 901 , QC.testProperty "K identical components track the component" $ 902 QC.forAll (QC.choose (1, 6)) $ \k -> 903 QC.forAll (QC.listOf unit_double) $ \xs -> 904 let bcfg = ok (Bounded.config 0.5 0.0 1.0 1.0e-3 Bounded.Newton) 905 xcfg = ok (Mix.config k 1.0e-3) 906 sts = drop 1 (scanl (Bounded.update bcfg) 907 (Bounded.initial bcfg) xs) 908 les = map Bounded.log_evalue sts 909 mix = L.foldl' 910 (\acc l -> Mix.update xcfg acc (replicate k l)) 911 (Mix.initial xcfg) les 912 cfin = L.foldl' (Bounded.update bcfg) (Bounded.initial bcfg) xs 913 in approx_eq (Mix.log_evalue xcfg mix) 914 (Bounded.log_evalue cfin) 915 && approx_eq (Mix.log_evalue_sup xcfg mix) 916 (Bounded.log_evalue_sup cfin) 917 918 , QC.testProperty "decide agrees with p_value at alpha" $ 919 QC.forAll (QC.choose (1, 6)) $ \k -> 920 QC.forAll (QC.listOf (QC.vectorOf k (QC.choose (-5, 5)))) $ \vs -> 921 let alpha = 0.5 922 cfg = ok (Mix.config k alpha) 923 sts = scanl (Mix.update cfg) (Mix.initial cfg) vs 924 in all (\s -> (Mix.decide cfg s == Mix.Reject) 925 == (Mix.p_value cfg s <= alpha)) sts 926 927 , QC.testProperty "sup monotone nondecreasing, verdict latched" $ 928 QC.forAll (QC.choose (1, 6)) $ \k -> 929 QC.forAll (QC.listOf (QC.vectorOf k (QC.choose (-5, 5)))) $ \vs -> 930 let cfg = ok (Mix.config k 0.5) 931 sts = scanl (Mix.update cfg) (Mix.initial cfg) vs 932 sups = map Mix.log_wealth_sup sts 933 in and (zipWith (<=) sups (drop 1 sups)) 934 && monotone_reject_bounded (map (Mix.decide cfg) sts) 935 936 , testCase "FPR under H_0 within slack (sign + magnitude hedge)" $ do 937 let rate = mixture_rate 0.05 0.5 2000 200 424242 938 assertBool ("FPR " ++ show rate ++ " exceeded slack") $ 939 rate <= 0.08 940 941 , testCase "power against p = 0.7 (sign + magnitude hedge)" $ do 942 let rate = mixture_rate 1.0e-3 0.7 5000 100 434343 943 assertBool ("power " ++ show rate ++ " too low") $ 944 rate >= 0.95 945 ] 946 -- confidence sequences ------------------------------------------------------- 947 -- a finite stream of bernoulli(p) samples. 948 cs_stream :: Double -> Int -> Gen -> [Double] 949 cs_stream !p n g0 = go n g0 950 where 951 go 0 _ = [] 952 go !k !g = 953 let (x, g') = bernoulli p g 954 in x : go (k - 1) g' 955 956 -- do the intervals nest: each contained in its predecessor, with 957 -- Nothing (empty) absorbing? 958 cs_nested :: [Maybe (Double, Double)] -> Bool 959 cs_nested ivs = and (zipWith shrink ivs (drop 1 ivs)) 960 where 961 shrink (Just (l1, u1)) (Just (l2, u2)) = l2 >= l1 && u2 <= u1 962 shrink (Just _) Nothing = True 963 shrink Nothing Nothing = True 964 shrink Nothing (Just _) = False 965 966 -- fraction of trials in which the true mean ever escapes the running 967 -- interval (or the interval goes empty), checked after every 968 -- observation. 969 cs_miscoverage_rate 970 :: CS.Config 971 -> Double -- ^ true mean 972 -> Int -- ^ budget per trial 973 -> Int -- ^ number of trials 974 -> Word64 -- ^ seed 975 -> Double 976 cs_miscoverage_rate cfg p budget trials seed = 977 let gens = take trials (gen_seq (mk_gen seed)) 978 misses = length [ () | g <- gens, cs_trial_missed g ] 979 in fromIntegral misses / fromIntegral trials 980 where 981 cs_trial_missed g0 = go budget g0 (CS.initial cfg) 982 where 983 go !k !g !st 984 | k == 0 = False 985 | otherwise = 986 let (x, g') = bernoulli p g 987 st' = CS.update cfg st x 988 in case CS.interval cfg st' of 989 Nothing -> True 990 Just (l, u) 991 | p < l || p > u -> True 992 | otherwise -> go (k - 1) g' st' 993 994 confseq_tests :: TestTree 995 confseq_tests = testGroup "confidence sequences" [ 996 testCase "initial interval is the full range" $ do 997 let cfg = ok (CS.config 0.0 1.0 0.05 100) 998 CS.interval cfg (CS.initial cfg) @?= Just (0.0, 1.0) 999 , testCase "intervals nest along a deterministic stream" $ do 1000 let cfg = ok (CS.config 0.0 1.0 0.05 50) 1001 xs = take 500 (cycle [1.0, 1.0, 0.0, 1.0]) 1002 sts = scanl (CS.update cfg) (CS.initial cfg) xs 1003 ivs = map (CS.interval cfg) sts 1004 assertBool "nesting violated" (cs_nested ivs) 1005 -- the stream has empirical mean 0.75; the final interval must 1006 -- be a strict refinement of the initial one. 1007 case (ivs, reverse ivs) of 1008 (iv0 : _, ivn : _) -> assertBool "no shrinkage" (iv0 /= ivn) 1009 _ -> assertFailure "no intervals" 1010 , QC.testProperty "intervals nest along any admissible stream" $ 1011 QC.forAll (QC.listOf unit_double) $ \xs -> 1012 let cfg = ok (CS.config 0.0 1.0 0.05 25) 1013 sts = scanl (CS.update cfg) (CS.initial cfg) xs 1014 in cs_nested (map (CS.interval cfg) sts) 1015 , testCase "coverage: off-grid Bernoulli(0.437) at alpha = 0.05" $ do 1016 let cfg = ok (CS.config 0.0 1.0 0.05 100) 1017 rate = cs_miscoverage_rate cfg 0.437 1500 200 991199 1018 -- expected miscoverage <= 0.05; allow up to 0.08 slack for 1019 -- sampling variability over 200 trials. 1020 assertBool ("miscoverage " ++ show rate ++ " exceeded slack") $ 1021 rate <= 0.08 1022 , testCase "consistency: Bernoulli(0.3) interval shrinks onto mean" $ do 1023 let cfg = ok (CS.config 0.0 1.0 1.0e-3 200) 1024 xs = cs_stream 0.3 5000 (mk_gen 424242) 1025 st = L.foldl' (CS.update cfg) (CS.initial cfg) xs 1026 case CS.interval cfg st of 1027 Nothing -> assertFailure "interval empty" 1028 Just (l, u) -> do 1029 assertBool ("interval " ++ show (l, u) ++ " misses mean") $ 1030 l <= 0.3 && 0.3 <= u 1031 assertBool ("width " ++ show (u - l) ++ " too wide") $ 1032 u - l < 0.2 1033 , testCase "affine: mean recovered on [-5, 5]" $ do 1034 -- x = 4 w.p. 0.7, x = -4 w.p. 0.3: true mean 1.6, interior 1035 -- to the sample bounds and asymmetric about zero. 1036 let cfg = ok (CS.config (-5.0) 5.0 0.05 100) 1037 xs = [ if b == 1.0 then 4.0 else (-4.0) 1038 | b <- cs_stream 0.7 3000 (mk_gen 232323) ] 1039 st = L.foldl' (CS.update cfg) (CS.initial cfg) xs 1040 case CS.interval cfg st of 1041 Nothing -> assertFailure "interval empty" 1042 Just (l, u) -> do 1043 assertBool ("interval " ++ show (l, u) ++ " misses mean") $ 1044 l <= 1.6 && 1.6 <= u 1045 assertBool ("interval " ++ show (l, u) ++ " not refined") $ 1046 l > -5.0 && u < 5.0 1047 , testCase "config: grid size 0 rejected" $ 1048 assertLeftCS (CS.config 0.0 1.0 0.05 0) 1049 , testCase "config: negative grid size rejected" $ 1050 assertLeftCS (CS.config 0.0 1.0 0.05 (-3)) 1051 , testCase "config: alpha out of range rejected" $ do 1052 assertLeftCS (CS.config 0.0 1.0 0.0 100) 1053 assertLeftCS (CS.config 0.0 1.0 1.5 100) 1054 , testCase "config: lo >= hi rejected" $ 1055 assertLeftCS (CS.config 1.0 0.0 0.05 100) 1056 , testCase "config: non-finite inputs rejected" $ do 1057 let nan = 0 / 0 :: Double 1058 pInf = 1 / 0 :: Double 1059 assertLeftCS (CS.config nan 1.0 0.05 100) 1060 assertLeftCS (CS.config 0.0 pInf 0.05 100) 1061 assertLeftCS (CS.config 0.0 1.0 nan 100) 1062 , QC.testProperty "interval endpoints well-formed on any stream" $ 1063 QC.forAll (QC.listOf unit_double) $ \xs -> 1064 let cfg = ok (CS.config 0.0 1.0 0.05 25) 1065 st = L.foldl' (CS.update cfg) (CS.initial cfg) xs 1066 in case CS.interval cfg st of 1067 Nothing -> True 1068 Just (l, u) -> 1069 finite l && finite u && 0 <= l && l <= u && u <= 1 1070 ] 1071 where 1072 assertLeftCS :: Either C.ConfigError a -> Assertion 1073 assertLeftCS e = case e of 1074 Left _ -> pure () 1075 Right _ -> assertFailure "expected Left"