From 5017fba6cb52b8cb2422f066ab4b4ab37859938d Mon Sep 17 00:00:00 2001
From: Anton Soloviev <anton@praviel.com>
Date: Thu, 15 Oct 2026 00:12:58 -0700
Subject: [PATCH] Use jackknife standard errors for Fisher-Lee inference

---
 pycircstat2/correlation.py |  6 ++++--
 tests/test_correlation.py  | 21 +++++++++++++++++++++
 2 files changed, 25 insertions(+), 2 deletions(-)

diff --git a/pycircstat2/correlation.py b/pycircstat2/correlation.py
index 535fb87..2a01f6a 100644
--- a/pycircstat2/correlation.py
+++ b/pycircstat2/correlation.py
@@ -162,8 +162,10 @@ def _corr(a, b):
         m_raas = np.mean(raas)
         s2_raas = np.var(raas, ddof=1)
         z = norm.ppf(0.975)
-        lb = n * r - (n - 1) * m_raas - z * np.sqrt(s2_raas / n)
-        ub = n * r - (n - 1) * m_raas + z * np.sqrt(s2_raas / n)
+        # The pseudovalues are n*r - (n-1)*raas, so their SD is (n-1) times larger.
+        standard_error = (n - 1) * np.sqrt(s2_raas / n)
+        lb = n * r - (n - 1) * m_raas - z * standard_error
+        ub = n * r - (n - 1) * m_raas + z * standard_error
 
         reject = ~(lb <= 0 <= ub)
 
diff --git a/tests/test_correlation.py b/tests/test_correlation.py
index 69f10e0..89cb345 100644
--- a/tests/test_correlation.py
+++ b/tests/test_correlation.py
@@ -1,9 +1,30 @@
 import numpy as np
+from scipy.stats import norm
 
 from pycircstat2 import Circular, load_data
 from pycircstat2.correlation import circ_corrcc, circ_corrcl
 
 
+def test_fisher_lee_jackknife_uses_pseudovalue_standard_error():
+    data = load_data("milwaukee", source="jammalamadaka")
+    a = np.deg2rad(data["theta"].to_numpy())
+    b = np.deg2rad(data["psi"].to_numpy())
+    n = len(a)
+    r = circ_corrcc(a, b, method="fl")
+    leave_one_out = np.array([
+        circ_corrcc(np.delete(a, i), np.delete(b, i), method="fl")
+        for i in range(n)
+    ])
+    pseudovalues = n * r - (n - 1) * leave_one_out
+    half_width = norm.ppf(0.975) * pseudovalues.std(ddof=1) / np.sqrt(n)
+    lower = pseudovalues.mean() - half_width
+    upper = pseudovalues.mean() + half_width
+    assert lower < 0 < upper
+    result = circ_corrcc(a, b, method="fl", test=True)
+    assert result.r == r
+    assert not result.reject_null
+
+
 def test_circ_corrcc():
 
     # Example 27.20 (Zar, 2010)