From 2462aa25074df95245081e9fcbd939c5d8b70c84 Mon Sep 17 00:00:00 2001 From: Anton Soloviev Date: Thu, 15 Oct 2026 00:23:06 -0700 Subject: [PATCH] Preserve correlation for valid small-angle observations --- pycircstat2/correlation.py | 14 +++++++++--- tests/test_correlation_small_angles.py | 30 ++++++++++++++++++++++++++ 2 files changed, 41 insertions(+), 3 deletions(-) create mode 100644 tests/test_correlation_small_angles.py diff --git a/pycircstat2/correlation.py b/pycircstat2/correlation.py index 535fb87..4d410bb 100644 --- a/pycircstat2/correlation.py +++ b/pycircstat2/correlation.py @@ -107,6 +107,14 @@ def _coerce_angles( return array, None +def _has_zero_pairwise_variation(angles: np.ndarray) -> bool: + """Detect constant or antipodal directions without an absolute variance floor.""" + differences = angles - angles[0] + # sin(k*pi) has rounding residuals proportional to the argument's magnitude. + tolerance = np.finfo(float).eps * np.abs(differences) + return bool(np.all(np.abs(np.sin(differences)) <= tolerance)) + + def _circ_corrcc_fl( a: Union[Circular, np.ndarray, Sequence[float]], b: Union[Circular, np.ndarray, Sequence[float]], @@ -146,7 +154,7 @@ def _corr(a, b): diff_b = np.sin(np.subtract.outer(b, b)[np.triu_indices(len(b), k=1)]) num = np.sum(diff_a * diff_b) den = np.sqrt(np.sum(diff_a**2) * np.sum(diff_b**2)) - if np.isclose(den, 0.0): + if den == 0.0 or _has_zero_pairwise_variation(a) or _has_zero_pairwise_variation(b): raise ValueError("Degenerate data produced zero variance in pairwise differences.") return num / den @@ -223,7 +231,7 @@ def _circ_corrcc_js( Sb = np.sin(bbar) num = np.sum(Sa * Sb) den = np.sqrt(np.sum(Sa**2) * np.sum(Sb**2)) - if np.isclose(den, 0.0): + if den == 0.0 or _has_zero_pairwise_variation(a_alpha) or _has_zero_pairwise_variation(b_alpha): raise ValueError("Degenerate data produced zero variance around the mean direction.") r = num / den @@ -232,7 +240,7 @@ def _circ_corrcc_js( l20 = np.mean(Sa**2) l02 = np.mean(Sb**2) l22 = np.mean(Sa**2 * Sb**2) - if np.isclose(l22, 0.0): + if l22 == 0.0: raise ValueError("Degenerate data caused division by zero in variance term.") test_stat = np.sqrt(n) * np.sqrt(l20 * l02 / l22) * r p_val = 2 * (1 - norm.cdf(np.abs(test_stat))) diff --git a/tests/test_correlation_small_angles.py b/tests/test_correlation_small_angles.py new file mode 100644 index 0000000..01e80a7 --- /dev/null +++ b/tests/test_correlation_small_angles.py @@ -0,0 +1,30 @@ +import numpy as np +import pytest + +from pycircstat2.correlation import circ_corrcc + + +@pytest.mark.parametrize("method", ["fl", "js"]) +@pytest.mark.parametrize("test", [False, True]) +@pytest.mark.parametrize("scale", [1e-3, 1e-5, 1e-20]) +def test_nonconstant_small_angles_have_unit_self_correlation(method, test, scale): + angles = scale * np.array([0.0, 1.0, 4.0, 2.0, 3.0]) + result = circ_corrcc(angles, angles, method=method, test=test) + r = result.r if test else result + assert r == pytest.approx(1.0) + if method == "js" and test: + assert np.isfinite(result.test_stat) + assert 0.0 < result.p_value < 1.0 + + +@pytest.mark.parametrize("method", ["fl", "js"]) +@pytest.mark.parametrize("test", [False, True]) +@pytest.mark.parametrize("angles", [ + np.zeros(5), + np.full(5, 0.7), + np.array([0.0, np.pi, 0.0, np.pi, 0.0]), + np.array([0.7, 0.7 + 2 * np.pi, 0.7, 0.7 + 2 * np.pi, 0.7]), +]) +def test_constant_or_antipodal_angles_remain_degenerate(method, test, angles): + with pytest.raises(ValueError, match="zero variance"): + circ_corrcc(angles, angles, method=method, test=test)