Mixture.hs (10541B)
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 qualified Data.List as L 109 import Numeric.Eproc.Common (Verdict(..), ConfigError(..), finite) 110 111 -- types ---------------------------------------------------------------------- 112 113 -- | Mixture configuration. Build with 'config'. 114 -- 115 -- Carries the component count @K@, the significance level, the 116 -- precomputed rejection threshold @log(K \/ alpha)@, and @log K@ 117 -- (the mixture log-wealth of a fresh state). 118 data Config = Config { 119 -- ^ component count @K@ 120 cfg_k :: {-# UNPACK #-} !Int 121 -- ^ significance level @alpha@ 122 , cfg_alpha :: {-# UNPACK #-} !Double 123 -- ^ rejection threshold @log(K \/ alpha)@ 124 , cfg_log_thresh :: {-# UNPACK #-} !Double 125 -- ^ @log K@ 126 , cfg_log_k :: {-# UNPACK #-} !Double 127 } 128 129 -- | Streaming mixture state. Construct with 'initial' and fold 130 -- per-step component log e-value vectors through 'update'. 131 -- 132 -- Tracks the current mixture log-wealth @log(sum_i E^i_t)@ and 133 -- its latched supremum, which is what 'decide' tests against the 134 -- rejection threshold. 135 data State = State { 136 st_n :: {-# UNPACK #-} !Int -- ^ update count 137 , st_log_sum :: {-# UNPACK #-} !Double -- ^ log(sum_i E^i) 138 , st_sup_log_sum :: {-# UNPACK #-} !Double -- ^ sup of the above 139 } 140 141 -- construction --------------------------------------------------------------- 142 143 -- | Build a 'Config' for a @K@-component uniform mixture at level 144 -- @alpha@. 145 -- 146 -- The rejection threshold is precomputed as @log(K \/ alpha)@: 147 -- the mixture @M_t = (sum_i E^i_t) \/ K@ crosses @1 \/ alpha@ 148 -- exactly when the sum crosses @K \/ alpha@. 149 -- 150 -- Returns 'Left' with a 'ConfigError' on inputs outside the 151 -- mathematical regime: @K < 1@, or @alpha@ non-finite or outside 152 -- @(0, 1)@. 153 -- 154 -- >>> let Right cfg = config 4 1.0e-3 155 config 156 :: Int -- ^ component count @K@ 157 -> Double -- ^ significance level @alpha@ 158 -> Either ConfigError Config 159 config !k !alpha 160 | k < 1 = 161 Left (InvalidComponentCount k) 162 | not (finite alpha && alpha > 0 && alpha < 1) = 163 Left (InvalidAlpha alpha) 164 | otherwise = 165 let !kd = fromIntegral k 166 in Right Config { 167 cfg_k = k 168 , cfg_alpha = alpha 169 , cfg_log_thresh = log (kd / alpha) 170 , cfg_log_k = log kd 171 } 172 {-# INLINE config #-} 173 174 -- | The initial 'State' for a fresh mixture. 175 -- 176 -- Every component starts at e-value @1@, so the mixture log-sum 177 -- (and its supremum) starts at @log K@. 178 -- 179 -- >>> let s0 = initial cfg 180 initial :: Config -> State 181 initial Config{..} = State { 182 st_n = 0 183 , st_log_sum = cfg_log_k 184 , st_sup_log_sum = cfg_log_k 185 } 186 {-# INLINE initial #-} 187 188 -- streaming ------------------------------------------------------------------ 189 190 -- | Fold one step's component log e-values into the running 191 -- 'State': computes the current mixture log-sum via a numerically 192 -- stable log-sum-exp and latches its supremum. 193 -- 194 -- /Preconditions/ (documented in the module header, unchecked 195 -- here): the vector holds exactly the @K@ log e-values of 196 -- components adapted to a common filtration, in a fixed order, 197 -- with 'update' called once per underlying observation. The 198 -- degenerate empty vector leaves the state unchanged. 199 -- 200 -- >>> let s1 = update cfg s0 [0.1, -0.2, 0.0, 0.4] 201 update :: Config -> State -> [Double] -> State 202 update _ st@State{..} les = case les of 203 [] -> st 204 (l : ls) -> 205 let !m = L.foldl' max l ls 206 !s = L.foldl' (\ !acc v -> acc + exp (v - m)) 0 les 207 -- all components at e-value zero: the mixture log-sum is 208 -- -Infinity, and (m +) would poison it into NaN. 209 !cur | isInfinite m && m < 0 = m 210 | otherwise = m + log s 211 in State { 212 st_n = st_n + 1 213 , st_log_sum = cur 214 , st_sup_log_sum = max st_sup_log_sum cur 215 } 216 {-# INLINE update #-} 217 218 -- | Compute the current 'Verdict' from the running 'State'. 219 -- 220 -- 'Reject' iff the supremum-so-far of @log(sum_i E^i_t)@ has ever 221 -- crossed @log(K \/ alpha)@ -- equivalently, the mixture 222 -- e-process @M_t@ has exceeded @1 \/ alpha@ at some point in the 223 -- stream so far. Under the combined @H_0@, by Ville's inequality, 224 -- the probability of this ever happening is at most @alpha@, 225 -- simultaneously over all sample sizes: peek and stop freely. 226 -- 227 -- >>> decide cfg s0 228 -- Continue 229 decide :: Config -> State -> Verdict 230 decide Config{..} State{..} 231 | st_sup_log_sum >= cfg_log_thresh = Reject 232 | otherwise = Continue 233 {-# INLINE decide #-} 234 235 -- inspection ----------------------------------------------------------------- 236 237 -- | The current mixture log-wealth @log(sum_i E^i_t)@, before 238 -- normalization by @K@. Not monotone; bounded above by 239 -- 'log_wealth_sup'. Starts at @log K@. 240 -- 241 -- >>> log_wealth s0 242 -- 1.3862943611198906 243 log_wealth :: State -> Double 244 log_wealth = st_log_sum 245 {-# INLINE log_wealth #-} 246 247 -- | The supremum-so-far of @log(sum_i E^i_t)@. Monotone 248 -- nondecreasing; 'decide' rejects exactly when it crosses 249 -- @log(K \/ alpha)@. Starts at @log K@. 250 -- 251 -- >>> log_wealth_sup s0 252 -- 1.3862943611198906 253 log_wealth_sup :: State -> Double 254 log_wealth_sup = st_sup_log_sum 255 {-# INLINE log_wealth_sup #-} 256 257 -- | The current log e-value of the mixture: the log of 258 -- @M_t = (sum_i E^i_t) \/ K@, i.e. 'log_wealth' minus @log K@, 259 -- normalized so a fresh state sits at @0@. This is itself a 260 -- component-shaped quantity: mixtures nest, so it can in turn be 261 -- fed to an outer mixture. Not monotone; bounded above by 262 -- 'log_evalue_sup'. 263 -- 264 -- >>> log_evalue s0 265 -- 0.0 266 log_evalue :: Config -> State -> Double 267 log_evalue Config{..} State{..} = st_log_sum - cfg_log_k 268 {-# INLINE log_evalue #-} 269 270 -- | The supremum-so-far of the log e-value: 'log_wealth_sup' minus 271 -- @log K@. Monotone nondecreasing, starting at @0@; 'decide' 272 -- rejects exactly when it crosses @log(1 \/ alpha)@. 273 -- 274 -- >>> log_evalue_sup s0 275 -- 0.0 276 log_evalue_sup :: Config -> State -> Double 277 log_evalue_sup Config{..} State{..} = st_sup_log_sum - cfg_log_k 278 {-# INLINE log_evalue_sup #-} 279 280 -- | The anytime-valid p-value: the reciprocal of the largest 281 -- mixture e-value attained so far. Monotone nonincreasing; under 282 -- the combined @H_0@, @P(exists t: p_t <= alpha) <= alpha@ for 283 -- every @alpha@ simultaneously. 'decide' returns 'Reject' exactly 284 -- when this value has reached the configured @alpha@ or below. 285 -- 286 -- >>> p_value cfg s0 287 -- 1.0 288 p_value :: Config -> State -> Double 289 p_value cfg s = min 1 (exp (negate (log_evalue_sup cfg s))) 290 {-# INLINE p_value #-} 291 292 -- | The number of 'update' steps consumed so far. 293 -- 294 -- >>> samples s0 295 -- 0 296 samples :: State -> Int 297 samples = st_n 298 {-# INLINE samples #-}