fix(core): Caveat-A follow-up — HAC estimator + multi-block gate + subset machine-check (Soraya's holes)

src/Core/AntiSybil.fs

@@ -232,6 +232,23 @@ module AntiSybil =
                 let num = [ for t in 0 .. n - 2 -> (x.[t] - mean) * (x.[t + 1] - mean) ] |> List.sum
                 num / denom
 
+    /// **Lag-`k` autocorrelation** `ρ_k` (Pearson between `xₜ` and `xₜ₊ₖ`). `0.0` for `k < 1`, `k ≥ n`,
+    /// or a constant series. Generalises `lag1Autocorr` (= `lagKAutocorr series 1`) — the HAC estimator
+    /// sums over lags to catch dependence that lag-1 alone misses (Soraya's lag-2 hole, 2026-08-04).
+    let lagKAutocorr (series: float list) (k: int) : float =
+        let x = List.toArray series
+        let n = x.Length
+        if k < 1 || k >= n then
+            0.0
+        else
+            let mean = Array.average x
+            let denom = x |> Array.sumBy (fun v -> (v - mean) * (v - mean))
+            if denom <= 1e-12 then
+                0.0
+            else
+                let num = [ for t in 0 .. n - 1 - k -> (x.[t] - mean) * (x.[t + k] - mean) ] |> List.sum
+                num / denom
+
     /// The round-ordered **±1 outcome-product series** that feeds the CHSH buckets:
     /// `prodₜ = sign(aₜ.Outcome) · sign(bₜ.Outcome)` (truncated to the shorter stream). Its
     /// autocorrelation is what drives `n_eff`.
@@ -254,15 +271,44 @@ module AntiSybil =
         let r = min RhoMax (max 0.0 rho1)
         float n * (1.0 - r) / (1.0 + r)
 
-    /// The **autocorrelation-corrected CHSH margin**: substitute `n_eff` (from the pair's own
-    /// outcome-product autocorrelation) for `n` in the Hoeffding ε. On an i.i.d. stream (`ρ₁ ≤ 0`) it
-    /// **equals** `chshMargin delta n`; on a positively-autocorrelated stream it is strictly **larger**
-    /// (`n_eff < n`) — the sound correction for Caveat (a). Takes the actual streams (not just `n`)
-    /// because `ρ₁` depends on the outcomes. `n_eff < 1` ⇒ `infinity` (no effective power ⇒ never convict).
+    /// **Newey–West bandwidth** `L = ⌊n^(1/3)⌋` (≥ 1) — the lag horizon summed in the long-run variance.
+    let neweyWestBandwidth (n: int) : int =
+        if n < 2 then 1 else max 1 (int (floor (float n ** (1.0 / 3.0))))
+
+    /// **Effective sample size under a Bartlett-windowed long-run variance** (Newey–West 1987) — the HAC
+    /// generalisation of the AR(1) `effectiveSampleSize`. It sums dependence across lags, so a stream with
+    /// weak `ρ₁` but strong higher-lag structure still shrinks `n_eff` (**fixes Soraya's lag-2 hole**):
+    ///   `n_eff = n / (1 + 2·Σ_{k=1}^{L} (1 − k/(L+1))·ρ_k⁺)`,  `ρ_k⁺ = max(0, ρ_k)`, `L = bandwidth`.
+    /// The clamped `ρ_k⁺ ≥ 0` and Bartlett weights keep the factor `≥ 1` ⇒ **`n_eff ≤ n` always, equality
+    /// iff every `ρ_k⁺ = 0`** — the monotonicity obligation, generalised past lag 1 (so Soraya's (a)/(b)/(c)
+    /// proofs still hold: they need only `n_eff ≤ n`). AR(1) is the special case `ρ_k = ρ₁^k`.
+    /// Anchors: Newey–West 1987 (HAC long-run variance); Bartlett 1946 (the psd kernel that guarantees
+    /// the factor is well-behaved); Kontorovich–Ramanan 2008 (concentration under mixing).
+    let effectiveSampleSizeHAC (series: float list) (bandwidth: int) : float =
+        let n = List.length series
+        if n < 2 then
+            float n
+        else
+            let l = max 1 bandwidth
+            let factor =
+                1.0
+                + 2.0
+                  * ([ for k in 1 .. min l (n - 1) ->
+                           let w = 1.0 - float k / float (l + 1)
+                           w * max 0.0 (lagKAutocorr series k) ]
+                     |> List.sum)
+            float n / factor
+
+    /// The **autocorrelation-corrected CHSH margin**: substitute `n_eff` (a Bartlett-windowed HAC estimate
+    /// over the pair's own outcome-product series, `effectiveSampleSizeHAC`) for `n` in the Hoeffding ε.
+    /// On a stream with no positive autocorrelation at any lag it **equals** `chshMargin delta n`; on any
+    /// positively-autocorrelated stream (at any lag ≤ bandwidth) it is strictly **larger** (`n_eff < n`) —
+    /// the sound correction for Caveat (a), now robust past lag 1. Takes the actual streams (not just `n`)
+    /// because the `ρ_k` depend on the outcomes. `n_eff < 1` ⇒ `infinity` (no effective power ⇒ never convict).
     let chshMarginAutocorr (delta: float) (a: ChshRound list) (b: ChshRound list) : float =
         let series = outcomeProductSeries a b
         let n = List.length series
-        let nEff = effectiveSampleSize n (lag1Autocorr series)
+        let nEff = effectiveSampleSizeHAC series (neweyWestBandwidth n)
         if nEff < 1.0 || delta <= 0.0 || delta >= 1.0 then infinity
         else sqrt (32.0 * log (1.0 / delta) / nEff)
 
@@ -281,6 +327,32 @@ module AntiSybil =
             let m2 = x.[h..] |> Array.average
             abs (m1 - m2) <= tol
 
+    /// **Multi-block stationarity gate** — strengthens `isApproxStationary`. Splits into `blocks`
+    /// contiguous blocks and requires **both** the spread of block MEANS and the spread of block
+    /// (population) VARIANCES to be `≤ tol`. The two-halves check saw only a first-moment difference of two
+    /// coarse blocks, so a within-half regime change or symmetric drift whose half-means cancel slipped
+    /// through (Soraya's defeat witness `[+1×10, −1×10, +1×10, −1×10]`, both half-means 0 yet a step
+    /// function). More blocks + a variance check catch it. A series shorter than `blocks` counts as
+    /// stationary (too little to compare). `blocks` is floored at 2.
+    let isApproxStationaryMultiBlock (tol: float) (blocks: int) (series: float list) : bool =
+        let x = List.toArray series
+        let n = x.Length
+        let b = max 2 blocks
+        if n < b then
+            true
+        else
+            let sz = n / b
+            let stats =
+                [ for i in 0 .. b - 1 ->
+                      let lo = i * sz
+                      let hi = if i = b - 1 then n - 1 else lo + sz - 1
+                      let seg = x.[lo..hi]
+                      let m = Array.average seg
+                      let v = seg |> Array.averageBy (fun z -> (z - m) * (z - m))
+                      m, v ]
+            let spread xs = List.max xs - List.min xs
+            spread (stats |> List.map fst) <= tol && spread (stats |> List.map snd) <= tol
+
     /// The CALIBRATED CHSH identity oracle: conviction at `2 + ε(n, δ)` with the
     /// pair's own run length, so an honestly-local pair at the bound is falsely
     /// convicted with probability ≤ δ (per pair) — the sound default the bare-
@@ -332,8 +404,9 @@ module AntiSybil =
         for i in 0 .. k - 1 do
             for j in i + 1 .. k - 1 do
                 let series = outcomeProductSeries arr.[i] arr.[j]
-                // Stationarity gate first: a non-stationary window is Unmeasured, never convicting.
-                if isApproxStationary stationarityTol series
+                // Stationarity gate first (multi-block: catches within-half drift the two-halves check
+                // missed): a non-stationary window is Unmeasured, never convicting.
+                if isApproxStationaryMultiBlock stationarityTol 4 series
                    && abs (chshS arr.[i] arr.[j]) > 2.0 + chshMarginAutocorr delta arr.[i] arr.[j] then
                     union i j
 

tests/Tests.FSharp/AntiSybil.Tests.fs

@@ -291,27 +291,62 @@ let ``effectiveSampleSize: n at rho<=0, strictly below n for rho>0, monotone dec
 // THE soundness property: the corrected margin is NEVER smaller than the i.i.d. margin, and STRICTLY
 // larger on a positively-autocorrelated outcome-product stream (n_eff < n). This is the Caveat-(a) fix.
 [<Fact>]
-let ``chshMarginAutocorr: >= the iid margin always, and STRICTLY larger on autocorrelated products`` () =
+let ``chshMarginAutocorr: >= iid always; STRICTLY larger on lag-1 AND lag-2 autocorrelation (HAC)`` () =
     let n = 64
-    // Runs product: a all +1, b in runs of 8 (+1.. then -1..) ⇒ product = b ⇒ positively autocorrelated.
     let a = [ for _ in 1..n -> { Setting = 0; Outcome = 1 } ]
+    // Runs of 8 ⇒ lag-1 positive ⇒ strictly larger.
     let bRuns = [ for i in 0 .. n - 1 -> { Setting = 0; Outcome = (if (i / 8) % 2 = 0 then 1 else -1) } ]
-    Assert.True(chshMarginAutocorr 0.05 a bRuns > chshMargin 0.05 n) // autocorrelated ⇒ strictly larger
-    // Alternating product (rho1 < 0 ⇒ clamped to 0 ⇒ n_eff = n) ⇒ EQUALS the i.i.d. margin.
+    Assert.True(chshMarginAutocorr 0.05 a bRuns > chshMargin 0.05 n)
+    // Alternating ⇒ lag-1 NEGATIVE but lag-2 = +1. The old lag-1-only model saw n_eff = n (equal margin);
+    // the HAC estimator catches the lag-2 structure ⇒ STRICTLY larger. This is Soraya's 3b fix in action.
     let bAlt = [ for i in 0 .. n - 1 -> { Setting = 0; Outcome = (if i % 2 = 0 then 1 else -1) } ]
-    Assert.Equal(chshMargin 0.05 n, chshMarginAutocorr 0.05 a bAlt, 9)
+    Assert.True(chshMarginAutocorr 0.05 a bAlt > chshMargin 0.05 n)
+    // A decorrelated pseudo-random product ⇒ never smaller than the i.i.d. margin (obligation b).
+    let bRand = [ for i in 0 .. n - 1 -> { Setting = 0; Outcome = (if (bits 5 n).[i] = 1 then 1 else -1) } ]
+    Assert.True(chshMarginAutocorr 0.05 a bRand >= chshMargin 0.05 n)
 
 [<Fact>]
-let ``isApproxStationary: stable series passes, a mean-shift fails`` () =
-    Assert.True(isApproxStationary 0.1 [ 1.0; -1.0; 1.0; -1.0 ]) // both halves mean 0
-    Assert.False(isApproxStationary 0.5 [ 1.0; 1.0; 1.0; 1.0; -1.0; -1.0; -1.0; -1.0 ]) // mean 1 → -1
+let ``lagKAutocorr: alternating has rho1 < 0 but rho2 = +1 (the lag-2 structure lag-1 misses)`` () =
+    let alt = [ for i in 0..19 -> if i % 2 = 0 then 1.0 else -1.0 ]
+    Assert.True(lagKAutocorr alt 1 < 0.0) // adjacent products anti-correlate
+    // biased estimator divides by full n ⇒ ρ₂ = (n−2)/n = 0.9 here (→ 1 as n grows); strongly positive.
+    Assert.True(lagKAutocorr alt 2 > 0.8)
+    Assert.True(lagKAutocorr alt 2 > lagKAutocorr alt 1) // lag-2 ≫ lag-1: the structure lag-1 misses
+
+// Soraya's lag-2 hole, fixed: HAC shrinks n_eff on a stream with weak rho1 but strong rho2, where the
+// AR(1) lag-1 estimator barely moves.
+[<Fact>]
+let ``effectiveSampleSizeHAC: catches lag-2 dependence the AR(1) lag-1 estimator misses`` () =
+    // +1 on evens (persistent), pseudo-random on odds ⇒ weak rho1, strong rho2.
+    let series = [ for i in 0..199 -> if i % 2 = 0 then 1.0 else (if (bits 3 200).[i] = 1 then 1.0 else -1.0) ]
+    let arNeff = effectiveSampleSize 200 (lag1Autocorr series) // lag-1 only ⇒ barely shrinks
+    let hacNeff = effectiveSampleSizeHAC series (neweyWestBandwidth 200) // sums lags ⇒ shrinks meaningfully
+    Assert.True(hacNeff <= 200.0) // obligation (a) still holds under HAC
+    Assert.True(arNeff > 180.0) // lag-1 barely moved (the hole)
+    Assert.True(hacNeff < arNeff) // HAC sees the lag-2 dependence AR(1)-lag1 missed
 
-// The oracle is STRICTLY more conservative: it can only keep MORE identities distinct (drop false
-// collapses), never fewer, than the i.i.d.-calibrated oracle — for any batch, any delta.
 [<Fact>]
-let ``chshSybilAutocorrCalibrated: never collapses MORE than the iid-calibrated oracle`` () =
-    let mk seed n = [ for i in 0 .. n - 1 -> { Setting = (bits seed n).[i]; Outcome = (if (bits (seed + 7) n).[i] = 1 then 1 else -1) } ]
-    let streams = [ mk 1 128; mk 2 128; mk 3 128; mk 1 128 ] // last is a replay of the first
-    let iid = chshSybilCalibrated 0.05 streams
-    let corrected = chshSybilAutocorrCalibrated 0.05 0.5 streams
-    Assert.True(corrected.DistinctCount >= iid.DistinctCount) // more conservative ⇒ >= distinct sources
+let ``isApproxStationary(+MultiBlock): the multi-block gate catches a step-function the two-halves check passes`` () =
+    Assert.True(isApproxStationary 0.1 [ 1.0; -1.0; 1.0; -1.0 ]) // stable
+    Assert.False(isApproxStationary 0.5 [ 1.0; 1.0; 1.0; 1.0; -1.0; -1.0; -1.0; -1.0 ]) // level shift caught
+    // Soraya's defeat witness: two-halves means both 0 ⇒ FALSE-passes; 4-block means +1,-1,+1,-1 ⇒ caught.
+    let witness =
+        [ for _ in 1..10 -> 1.0 ] @ [ for _ in 1..10 -> -1.0 ] @ [ for _ in 1..10 -> 1.0 ] @ [ for _ in 1..10 -> -1.0 ]
+    Assert.True(isApproxStationary 0.1 witness) // two-halves is fooled
+    Assert.False(isApproxStationaryMultiBlock 0.1 4 witness) // multi-block is not
+
+// Machine-check of obligation (c): the autocorr oracle's same-source equivalence is a REFINEMENT (subset)
+// of the i.i.d. oracle's ⇒ DistinctCount_autocorr >= DistinctCount_iid, over many generated batches
+// (including replays and autocorrelated streams). Strictly-more-conservative, never a new false collapse.
+[<Fact>]
+let ``chshSybilAutocorrCalibrated: conviction set SUBSET of the iid oracle over many batches`` () =
+    let mk s n = [ for i in 0 .. n - 1 -> { Setting = (bits s n).[i]; Outcome = (if (bits (s + 7) n).[i] = 1 then 1 else -1) } ]
+    let runs s n = [ for i in 0 .. n - 1 -> { Setting = (bits s n).[i]; Outcome = (if (i / 8) % 2 = 0 then 1 else -1) } ]
+    for seed in 1..40 do
+        let streams = [ mk seed 96; mk (seed + 1) 96; mk seed 96; runs (seed + 2) 96; runs (seed + 2) 96 ]
+        let iid = chshSybilCalibrated 0.05 streams
+        let corrected = chshSybilAutocorrCalibrated 0.05 0.5 streams
+        Assert.True(
+            corrected.DistinctCount >= iid.DistinctCount,
+            sprintf "seed %d: autocorr distinct %d < iid distinct %d" seed corrected.DistinctCount iid.DistinctCount
+        )