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