diff --git a/swarmmetrics.py b/swarmmetrics.py
index 98e9438..b6c7118 100644
--- a/swarmmetrics.py
+++ b/swarmmetrics.py
@@ -460,7 +460,7 @@ def shuffle_test(messages: list[dict], n_shuffles: int = 100,
         if std_s > 0:
             z = (obs - mean_s) / std_s
         else:
-            z = 0.0 if obs == mean_s else float('inf')
+            z = 0.0 if obs == mean_s else (10.0 if obs > mean_s else -10.0)
         result["z_scores"][node] = round(z, 4)
         result["significant"][node] = abs(z) > 2.0
 
@@ -504,12 +504,16 @@ def ablation_sensitivity(messages: list[dict], target: str,
     for node in set(z_full) | set(z_abl):
         if node == target:
             continue
-        deltas[node] = round(z_abl.get(node, 0) - z_full.get(node, 0), 4)
+        d = z_abl.get(node, 0) - z_full.get(node, 0)
+        # Clamp to avoid ±inf from zero-variance singletons
+        d = max(-100.0, min(100.0, d))
+        deltas[node] = round(d, 4)
 
     gainers = {k: v for k, v in deltas.items() if v > 2.0}
     losers = {k: v for k, v in deltas.items() if v < -2.0}
 
-    mean_delta = sum(deltas.values()) / max(len(deltas), 1)
+    finite_deltas = [v for v in deltas.values() if math.isfinite(v)]
+    mean_delta = sum(finite_deltas) / max(len(finite_deltas), 1)
     if mean_delta > 0.5:
         verdict = "suppressor"
     elif mean_delta < -0.5:
@@ -842,7 +846,9 @@ def phi_accrual(messages: list, t_now: float = None,
         # Compute inter-arrival intervals
         intervals = [times[i+1] - times[i] for i in range(len(times) - 1)]
         mean_ival = sum(intervals) / len(intervals)
-        variance = sum((x - mean_ival) ** 2 for x in intervals) / len(intervals)
+        # Sample variance (N-1) — population variance underestimates for small N
+        n_iv = len(intervals)
+        variance = sum((x - mean_ival) ** 2 for x in intervals) / max(n_iv - 1, 1)
         std_ival = math.sqrt(variance) if variance > 0 else mean_ival * 0.1
 
         liveness.mean_interval = mean_ival
@@ -859,15 +865,10 @@ def phi_accrual(messages: list, t_now: float = None,
             liveness.phi = 10.0 if t_diff > mean_ival * 1.5 else 0.0
             liveness.state = 'gray' if liveness.phi >= 3.0 else 'green'
         else:
-            # Normal CDF approximation (error function)
+            # Normal CDF via error function (exact, not logistic approx)
             # P(X <= t_diff) where X ~ N(mean, std)
             z = (t_diff - mean_ival) / std_ival
-            # Approximate CDF using logistic approximation
-            # F(z) ≈ 1 / (1 + exp(-1.7 * z))
-            try:
-                cdf = 1.0 / (1.0 + math.exp(-1.7 * z))
-            except OverflowError:
-                cdf = 1.0 if z > 0 else 0.0
+            cdf = 0.5 * (1.0 + math.erf(z / math.sqrt(2.0)))
 
             # φ = -log10(1 - F(t_diff))
             if cdf >= 1.0 - 1e-15: