fix(core): normalize Q statistics by residual count

src/pykrige/core.py

@@ -838,12 +838,12 @@ def _find_statistics(
 
 def calcQ1(epsilon):
     """Returns the Q1 statistic for the variogram fit (see [1])."""
-    return abs(np.sum(epsilon) / (epsilon.shape[0] - 1))
+    return abs(np.sum(epsilon) / epsilon.shape[0])
 
 
 def calcQ2(epsilon):
     """Returns the Q2 statistic for the variogram fit (see [1])."""
-    return np.sum(epsilon**2) / (epsilon.shape[0] - 1)
+    return np.sum(epsilon**2) / epsilon.shape[0]
 
 
 def calc_cR(Q2, sigma):

tests/test_core.py

@@ -375,6 +375,18 @@ def test_core_calculate_variogram_model():
     assert_allclose(res, np.array([3.0, 2.0, 1.0]), 0.001, 0.001)
 
 
+def test_q1_q2_use_number_of_epsilon_residuals():
+    epsilon = np.array([1.0, 1.0])
+    assert core.calcQ1(epsilon) == 1.0
+    assert core.calcQ2(epsilon) == 1.0
+
+
+def test_q1_q2_single_residual_are_finite():
+    epsilon = np.array([0.5])
+    assert core.calcQ1(epsilon) == 0.5
+    assert core.calcQ2(epsilon) == 0.25
+
+
 def test_core_krige():
     # Example 3.2 from Kitanidis
     data = np.array([[9.7, 47.6, 1.22], [43.8, 24.6, 2.822]])