Mixture.hs (10505B)
1 {-# OPTIONS_HADDOCK prune #-} 2 {-# LANGUAGE BangPatterns #-} 3 {-# LANGUAGE RecordWildCards #-} 4 5 -- | 6 -- Module: Numeric.Eproc.Mixture 7 -- Copyright: (c) 2026 Jared Tobin 8 -- License: MIT 9 -- Maintainer: Jared Tobin <jared@ppad.tech> 10 -- 11 -- Uniform convex mixture of e-processes. 12 -- 13 -- Given @K@ component e-processes @E^1_t, ..., E^K_t@ adapted to a 14 -- common filtration -- each testing (its facet of) a shared null 15 -- @H_0@ -- their arithmetic mean 16 -- 17 -- @M_t = (E^1_t + ... + E^K_t) \/ K@ 18 -- 19 -- is itself an e-process with @M_0 = 1@: convex combinations 20 -- preserve the nonnegative-supermartingale property. By Ville's 21 -- inequality @P(sup_t M_t >= 1 \/ alpha) <= alpha@ under @H_0@, so a 22 -- level-@alpha@ test of the /combined/ null rejects when 23 -- @sup_t log(E^1_t + ... + E^K_t)@ crosses @log(K \/ alpha)@ -- no 24 -- Bonferroni correction, and strictly more powerful than one, since 25 -- the sum dominates the max. Use a mixture when the alternative has 26 -- several qualitatively different faces (a location shift, a shape 27 -- change, a rare-outlier channel, ...) and you want a single test 28 -- with power against their union. 29 -- 30 -- This module does not own or update the components: they may be 31 -- heterogeneous (different test modules, different observation 32 -- transformations), so the caller steps each component itself and 33 -- feeds 'update' the vector of their current log e-values, as 34 -- reported by each module's @log_evalue@ accessor, one entry per 35 -- component in a fixed order. 36 -- 37 -- Two preconditions are the caller's responsibility, and the 38 -- type-I guarantee depends on both: 39 -- 40 -- 1. Each entry must be the current log e-value of a genuine 41 -- e-process for @H_0@, and all components must be adapted to 42 -- the same filtration and stepped in lockstep -- 'update' is 43 -- called exactly once per underlying observation, after all 44 -- components have absorbed it. 45 -- 46 -- 2. The vector must have exactly the @K@ entries declared in 47 -- 'config', always in the same order. 48 -- 49 -- The rejection latch is kept on the supremum of the /mixture's/ 50 -- log-wealth. Latching (or summing) per-component suprema instead 51 -- would combine peaks attained at different times -- a quantity 52 -- that can exceed anything the mixture ever reached, silently 53 -- inflating the effective alpha. Ville's inequality bounds the 54 -- mixture's own supremum; that is the only sound latch, and it is 55 -- the one this module maintains. 56 -- 57 -- == Example 58 -- 59 -- Combine a sign test and a magnitude test running against the same 60 -- stream of differences @d_t@ (the shape used for two-channel 61 -- symmetry testing): 62 -- 63 -- >>> import qualified Numeric.Eproc.Bernoulli.TwoSided as Sign 64 -- >>> import qualified Numeric.Eproc.Bounded as Magn 65 -- >>> import qualified Numeric.Eproc.Mixture as Mix 66 -- >>> let Right sc = Sign.config 0.5 1.0e-3 Sign.Newton 67 -- >>> let Right mc = Magn.config 0.0 (-1.0) 1.0 1.0e-3 Magn.Newton 68 -- >>> let Right xc = Mix.config 2 1.0e-3 69 -- >>> :{ 70 -- let step (s, m, x) d = 71 -- let s' = Sign.update sc s (d > 0) 72 -- m' = Magn.update mc m d 73 -- in (s', m', Mix.update xc x 74 -- [Sign.log_evalue s', Magn.log_evalue m']) 75 -- :} 76 -- >>> let ds = take 400 (cycle [0.6, 0.7, -0.2, 0.8]) 77 -- >>> let z0 = (Sign.initial sc, Magn.initial mc, Mix.initial xc) 78 -- >>> let (_, _, xf) = foldl' step z0 ds 79 -- >>> Mix.decide xc xf 80 -- Reject 81 -- >>> Mix.p_value xc xf 82 -- 9.482234479673792e-34 83 84 module Numeric.Eproc.Mixture ( 85 -- * Mixture configuration and state 86 Config 87 , State 88 , Verdict(..) 89 , ConfigError(..) 90 91 -- * Construction 92 , config 93 , initial 94 95 -- * Streaming 96 , update 97 , decide 98 99 -- * Inspection 100 , log_wealth 101 , log_wealth_sup 102 , log_evalue 103 , log_evalue_sup 104 , p_value 105 , samples 106 ) where 107 108 import Numeric.Eproc.Common (Verdict(..), ConfigError(..), finite) 109 110 -- types ---------------------------------------------------------------------- 111 112 -- | Mixture configuration. Build with 'config'. 113 -- 114 -- Carries the component count @K@, the significance level, the 115 -- precomputed rejection threshold @log(K \/ alpha)@, and @log K@ 116 -- (the mixture log-wealth of a fresh state). 117 data Config = Config { 118 -- ^ component count @K@ 119 cfg_k :: {-# UNPACK #-} !Int 120 -- ^ significance level @alpha@ 121 , cfg_alpha :: {-# UNPACK #-} !Double 122 -- ^ rejection threshold @log(K \/ alpha)@ 123 , cfg_log_thresh :: {-# UNPACK #-} !Double 124 -- ^ @log K@ 125 , cfg_log_k :: {-# UNPACK #-} !Double 126 } 127 128 -- | Streaming mixture state. Construct with 'initial' and fold 129 -- per-step component log e-value vectors through 'update'. 130 -- 131 -- Tracks the current mixture log-wealth @log(sum_i E^i_t)@ and 132 -- its latched supremum, which is what 'decide' tests against the 133 -- rejection threshold. 134 data State = State { 135 st_n :: {-# UNPACK #-} !Int -- ^ update count 136 , st_log_sum :: {-# UNPACK #-} !Double -- ^ log(sum_i E^i) 137 , st_sup_log_sum :: {-# UNPACK #-} !Double -- ^ sup of the above 138 } 139 140 -- construction --------------------------------------------------------------- 141 142 -- | Build a 'Config' for a @K@-component uniform mixture at level 143 -- @alpha@. 144 -- 145 -- The rejection threshold is precomputed as @log(K \/ alpha)@: 146 -- the mixture @M_t = (sum_i E^i_t) \/ K@ crosses @1 \/ alpha@ 147 -- exactly when the sum crosses @K \/ alpha@. 148 -- 149 -- Returns 'Left' with a 'ConfigError' on inputs outside the 150 -- mathematical regime: @K < 1@, or @alpha@ non-finite or outside 151 -- @(0, 1)@. 152 -- 153 -- >>> let Right cfg = config 4 1.0e-3 154 config 155 :: Int -- ^ component count @K@ 156 -> Double -- ^ significance level @alpha@ 157 -> Either ConfigError Config 158 config !k !alpha 159 | k < 1 = 160 Left (InvalidComponentCount k) 161 | not (finite alpha && alpha > 0 && alpha < 1) = 162 Left (InvalidAlpha alpha) 163 | otherwise = 164 let !kd = fromIntegral k 165 in Right Config { 166 cfg_k = k 167 , cfg_alpha = alpha 168 , cfg_log_thresh = log (kd / alpha) 169 , cfg_log_k = log kd 170 } 171 {-# INLINE config #-} 172 173 -- | The initial 'State' for a fresh mixture. 174 -- 175 -- Every component starts at e-value @1@, so the mixture log-sum 176 -- (and its supremum) starts at @log K@. 177 -- 178 -- >>> let s0 = initial cfg 179 initial :: Config -> State 180 initial Config{..} = State { 181 st_n = 0 182 , st_log_sum = cfg_log_k 183 , st_sup_log_sum = cfg_log_k 184 } 185 {-# INLINE initial #-} 186 187 -- streaming ------------------------------------------------------------------ 188 189 -- | Fold one step's component log e-values into the running 190 -- 'State': computes the current mixture log-sum via a numerically 191 -- stable log-sum-exp and latches its supremum. 192 -- 193 -- /Preconditions/ (documented in the module header, unchecked 194 -- here): the vector holds exactly the @K@ log e-values of 195 -- components adapted to a common filtration, in a fixed order, 196 -- with 'update' called once per underlying observation. The 197 -- degenerate empty vector leaves the state unchanged. 198 -- 199 -- >>> let s1 = update cfg s0 [0.1, -0.2, 0.0, 0.4] 200 update :: Config -> State -> [Double] -> State 201 update _ st@State{..} les = case les of 202 [] -> st 203 (l : ls) -> 204 let !m = foldl' max l ls 205 !s = foldl' (\ !acc v -> acc + exp (v - m)) 0 les 206 -- all components at e-value zero: the mixture log-sum is 207 -- -Infinity, and (m +) would poison it into NaN. 208 !cur | isInfinite m && m < 0 = m 209 | otherwise = m + log s 210 in State { 211 st_n = st_n + 1 212 , st_log_sum = cur 213 , st_sup_log_sum = max st_sup_log_sum cur 214 } 215 {-# INLINE update #-} 216 217 -- | Compute the current 'Verdict' from the running 'State'. 218 -- 219 -- 'Reject' iff the supremum-so-far of @log(sum_i E^i_t)@ has ever 220 -- crossed @log(K \/ alpha)@ -- equivalently, the mixture 221 -- e-process @M_t@ has exceeded @1 \/ alpha@ at some point in the 222 -- stream so far. Under the combined @H_0@, by Ville's inequality, 223 -- the probability of this ever happening is at most @alpha@, 224 -- simultaneously over all sample sizes: peek and stop freely. 225 -- 226 -- >>> decide cfg s0 227 -- Continue 228 decide :: Config -> State -> Verdict 229 decide Config{..} State{..} 230 | st_sup_log_sum >= cfg_log_thresh = Reject 231 | otherwise = Continue 232 {-# INLINE decide #-} 233 234 -- inspection ----------------------------------------------------------------- 235 236 -- | The current mixture log-wealth @log(sum_i E^i_t)@, before 237 -- normalization by @K@. Not monotone; bounded above by 238 -- 'log_wealth_sup'. Starts at @log K@. 239 -- 240 -- >>> log_wealth s0 241 -- 1.3862943611198906 242 log_wealth :: State -> Double 243 log_wealth = st_log_sum 244 {-# INLINE log_wealth #-} 245 246 -- | The supremum-so-far of @log(sum_i E^i_t)@. Monotone 247 -- nondecreasing; 'decide' rejects exactly when it crosses 248 -- @log(K \/ alpha)@. Starts at @log K@. 249 -- 250 -- >>> log_wealth_sup s0 251 -- 1.3862943611198906 252 log_wealth_sup :: State -> Double 253 log_wealth_sup = st_sup_log_sum 254 {-# INLINE log_wealth_sup #-} 255 256 -- | The current log e-value of the mixture: the log of 257 -- @M_t = (sum_i E^i_t) \/ K@, i.e. 'log_wealth' minus @log K@, 258 -- normalized so a fresh state sits at @0@. This is itself a 259 -- component-shaped quantity: mixtures nest, so it can in turn be 260 -- fed to an outer mixture. Not monotone; bounded above by 261 -- 'log_evalue_sup'. 262 -- 263 -- >>> log_evalue s0 264 -- 0.0 265 log_evalue :: Config -> State -> Double 266 log_evalue Config{..} State{..} = st_log_sum - cfg_log_k 267 {-# INLINE log_evalue #-} 268 269 -- | The supremum-so-far of the log e-value: 'log_wealth_sup' minus 270 -- @log K@. Monotone nondecreasing, starting at @0@; 'decide' 271 -- rejects exactly when it crosses @log(1 \/ alpha)@. 272 -- 273 -- >>> log_evalue_sup s0 274 -- 0.0 275 log_evalue_sup :: Config -> State -> Double 276 log_evalue_sup Config{..} State{..} = st_sup_log_sum - cfg_log_k 277 {-# INLINE log_evalue_sup #-} 278 279 -- | The anytime-valid p-value: the reciprocal of the largest 280 -- mixture e-value attained so far. Monotone nonincreasing; under 281 -- the combined @H_0@, @P(exists t: p_t <= alpha) <= alpha@ for 282 -- every @alpha@ simultaneously. 'decide' returns 'Reject' exactly 283 -- when this value has reached the configured @alpha@ or below. 284 -- 285 -- >>> p_value cfg s0 286 -- 1.0 287 p_value :: Config -> State -> Double 288 p_value cfg s = min 1 (exp (negate (log_evalue_sup cfg s))) 289 {-# INLINE p_value #-} 290 291 -- | The number of 'update' steps consumed so far. 292 -- 293 -- >>> samples s0 294 -- 0 295 samples :: State -> Int 296 samples = st_n 297 {-# INLINE samples #-}