From 43a90e5d2757007364afa6c0be954740a830713d Mon Sep 17 00:00:00 2001
From: Julien Chevallier <julien.chevallier1@univ-grenoble-alpes.fr>
Date: Sun, 26 Apr 2026 15:37:35 +0200
Subject: [PATCH 1/2] refactor: use prefix sums in scalar estimators

---
 src/estimation.jl | 11 +++++++++--
 1 file changed, 9 insertions(+), 2 deletions(-)

diff --git a/src/estimation.jl b/src/estimation.jl
index 7598723..5c5265b 100644
--- a/src/estimation.jl
+++ b/src/estimation.jl
@@ -22,18 +22,25 @@ function estimators(
 
     Z_T = sum(X; dims=2)
     ∑X = sum(X; dims=1)
+    cs_∑X = cumsum(vec(∑X))
     Z̄_T = mean(Z_T)
     m̂ = Z̄_T / T
     v̂ = (N) * (T + 1) * T^(-3) * (mean(Z_T .^ 2) - T / (T + 1) * (Z̄_T + Z̄_T^2))
     WΔ = 0.0
     for iter in 1:div(T, Δ)
-        WΔ += (N) / T * (sum(∑X[(1 + (iter - 1) * Δ):(iter * Δ)]) / (N) - Δ * m̂)^2
+        lo = 1 + (iter - 1) * Δ
+        hi = iter * Δ
+        block_sum = lo == 1 ? cs_∑X[hi] : cs_∑X[hi] - cs_∑X[lo - 1]
+        WΔ += (N) / T * (block_sum / (N) - Δ * m̂)^2
     end
     W2Δ = 0.0
     for iter in 1:div(T, 2 * Δ)
+        lo = 1 + (iter - 1) * 2 * Δ
+        hi = iter * 2 * Δ
+        block_sum = lo == 1 ? cs_∑X[hi] : cs_∑X[hi] - cs_∑X[lo - 1]
         W2Δ +=
             (N) / T *
-            (sum(∑X[(1 + (iter - 1) * 2 * Δ):(iter * 2 * Δ)]) / (N) - 2 * Δ * m̂)^2
+            (block_sum / (N) - 2 * Δ * m̂)^2
     end
     ŵ = 2 * W2Δ - WΔ
 

From f474639a969ff947f49bd25d706f98f886a02feb Mon Sep 17 00:00:00 2001
From: Julien Chevallier <julien.chevallier1@univ-grenoble-alpes.fr>
Date: Sun, 26 Apr 2026 15:38:58 +0200
Subject: [PATCH 2/2] refactor: use shared prefix sums in vector estimators

---
 src/estimation.jl | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/src/estimation.jl b/src/estimation.jl
index 5c5265b..dad75cd 100644
--- a/src/estimation.jl
+++ b/src/estimation.jl
@@ -38,9 +38,7 @@ function estimators(
         lo = 1 + (iter - 1) * 2 * Δ
         hi = iter * 2 * Δ
         block_sum = lo == 1 ? cs_∑X[hi] : cs_∑X[hi] - cs_∑X[lo - 1]
-        W2Δ +=
-            (N) / T *
-            (block_sum / (N) - 2 * Δ * m̂)^2
+        W2Δ += (N) / T * (block_sum / (N) - 2 * Δ * m̂)^2
     end
     ŵ = 2 * W2Δ - WΔ
 
@@ -59,6 +57,7 @@ function estimators(
 
     Z_T = sum(X; dims=2)
     ∑X = sum(X; dims=1)
+    cs_∑X = cumsum(vec(∑X))
     Z̄_T = mean(Z_T)
     m̂ = Z̄_T / T
     v̂ = (N) * (T + 1) * T^(-3) * (mean(Z_T .^ 2) - T / (T + 1) * (Z̄_T + Z̄_T^2))
@@ -69,13 +68,17 @@ function estimators(
         end
         WΔ = 0.0
         for iter in 1:div(T, Δ)
-            WΔ += (N) / T * (sum(∑X[(1 + (iter - 1) * Δ):(iter * Δ)]) / (N) - Δ * m̂)^2
+            lo = 1 + (iter - 1) * Δ
+            hi = iter * Δ
+            block_sum = lo == 1 ? cs_∑X[hi] : cs_∑X[hi] - cs_∑X[lo - 1]
+            WΔ += (N) / T * (block_sum / (N) - Δ * m̂)^2
         end
         W2Δ = 0.0
         for iter in 1:div(T, 2 * Δ)
-            W2Δ +=
-                (N) / T *
-                (sum(∑X[(1 + (iter - 1) * 2 * Δ):(iter * 2 * Δ)]) / (N) - 2 * Δ * m̂)^2
+            lo = 1 + (iter - 1) * 2 * Δ
+            hi = iter * 2 * Δ
+            block_sum = lo == 1 ? cs_∑X[hi] : cs_∑X[hi] - cs_∑X[lo - 1]
+            W2Δ += (N) / T * (block_sum / (N) - 2 * Δ * m̂)^2
         end
         push!(ŵ, 2 * W2Δ - WΔ)
     end