Repository: vnmabus/dcor

Issue:
Counting the distance from a point to itself
Hi, I've hit a bit of a problem. I was trying to work out why I was getting different results from the `ecp` R package versus `dcor`. After some intense investigation, I think the cause seems to be at the point of taking the mean of each within-sample distance. Note, this is before we apply the coefficient or consider the between-sample distances. Precisely, I'm referring to the mean taken here: https://github.com/vnmabus/dcor/blob/161a6f5928ec0f30ce89fcfd5e90e6ed9315e383/dcor/_energy.py#L41-L42

In all the Székely and Rizzo papers (e.g. Székely & Rizzo, 2004), this mean is defined as the arithmetic mean, and the same as you have used in `dcor`:
![image](https://user-images.githubusercontent.com/5019367/109004266-a52da400-76fc-11eb-8a4c-23254fbef42c.png)

However in the Matteson and James papers I have been looking at (e.g. Matteson & James, 2014; James et al., 2016), they seem to define it as follows:
 
![image](https://user-images.githubusercontent.com/5019367/109005122-b4612180-76fd-11eb-81c2-d90d78e6a277.png)

What they seem to be doing here is summing the lower triangle of the matrix, *excluding* the diagonal, and then divided by the combination n choose 2. So if we had a sample with 5 items, the full distance matrix would be 5 x 5 = 25 items, but the lower triangle would only have 10 items in it. They would sum these distances and divide by 5 choose 2, which is 10. So this is also taking the mean, but it's the mean *excluding* the diagonal, which is of course always 0 in a within-sample distance matrix. The ultimate outcome is that their "mean" is actually ![\frac{n}{n-1} \mu](https://render.githubusercontent.com/render/math?math=%5Cdisplaystyle+%5Cfrac%7Bn%7D%7Bn-1%7D+%5Cmu), which is larger than it should be, as it isn't counting the 0s on the diagonal.

Note that this is also visible in the implementation of their work, in the `ecp` package. Here, they sum the matrix but then divide by ![n \times n - 1](https://render.githubusercontent.com/render/math?math=%5Cdisplaystyle+n+%5Ctimes+n+-+1), which is equivalent to the above, but not equivalent to the true mean:
https://github.com/zwenyu/ecp/blob/65a9bb56308d25ce3c6be4d6388137f428118248/src/energyChangePoint.cpp#L112

My question is this: are they simply wrong? If no, is there any theory supporting this alternative formula? If there is, should this be something supported in `dcor`? Fortunately it kind of already is thanks to my customizable average feature. But it could be called out specifically. I appreciate your input here as you likely understand this domain better than I do.

Discussion:

Patch:
diff --git a/dcor/__init__.py b/dcor/__init__.py
index de2d8b6..8df9903 100644
--- a/dcor/__init__.py
+++ b/dcor/__init__.py
@@ -25,7 +25,7 @@ from ._dcor_internals import (double_centered, u_centered,  # noqa
                               mean_product, u_product,
                               u_projection,
                               u_complementary_projection)
-from ._energy import energy_distance  # noqa
+from ._energy import energy_distance, EstimationStatistic  # noqa
 from ._partial_dcor import (partial_distance_covariance,  # noqa
                             partial_distance_correlation)
 from ._rowwise import rowwise, RowwiseMode
diff --git a/dcor/_energy.py b/dcor/_energy.py
index f91cc39..548d33a 100644
--- a/dcor/_energy.py
+++ b/dcor/_energy.py
@@ -4,10 +4,45 @@ import warnings
 
 import numpy as np
 
+from enum import Enum, auto
+
 from . import distances
 from ._utils import _transform_to_2d
 
 
+class EstimationStatistic(Enum):
+    """
+    A type of estimation statistic used for calculating energy distance.
+    """
+    @classmethod
+    def from_string(cls, item):
+        """
+        Allows EstimationStatistic.from_string('u'),
+        EstimationStatistic.from_string('V'),
+        EstimationStatistic.from_string('V_STATISTIC'),
+        EstimationStatistic.from_string('u_statistic') etc
+        """
+        upper = item.upper()
+        if upper == 'U':
+            return cls.U_STATISTIC
+        elif upper == 'V':
+            return cls.V_STATISTIC
+        else:
+            return cls[upper]
+
+    U_STATISTIC = auto()
+    """
+    Hoeffding's unbiased U-statistics
+    (does not include the distance from each point to itself)
+    """
+
+    V_STATISTIC = auto()
+    """
+    von Mises's biased V-statistics
+    (does include the distance from each point to itself)
+    """
+
+
 def _check_valid_energy_exponent(exponent):
     if not 0 < exponent < 2:
         warning_msg = ('The energy distance is not guaranteed to be '
@@ -19,7 +54,8 @@ def _check_valid_energy_exponent(exponent):
 
 
 def _energy_distance_from_distance_matrices(
-        distance_xx, distance_yy, distance_xy, average=None):
+        distance_xx, distance_yy, distance_xy, average=None,
+        estimation_stat=EstimationStatistic.V_STATISTIC):
     """
     Compute energy distance with precalculated distance matrices.
 
@@ -28,11 +64,25 @@ def _energy_distance_from_distance_matrices(
     average: Callable[[ArrayLike], float]
         A function that will be used to calculate an average of distances.
         This defaults to np.mean.
-
+    estimation_stat: Union[str, EstimationStatistic]
+        If EstimationStatistic.U_STATISTIC, calculate energy distance using
+        Hoeffding's unbiased U-statistics. Otherwise, use von Mises's biased
+        V-statistics.
+        If this is provided as a string, it will first be converted to
+        an EstimationStatistic enum instance.
     """
+    if isinstance(estimation_stat, str):
+        estimation_stat = EstimationStatistic.from_string(estimation_stat)
+
     if average is None:
         average = np.mean
 
+    if estimation_stat == EstimationStatistic.U_STATISTIC:
+        # If using u-statistics, we exclude the central diagonal of 0s for the
+        # within-sample distances
+        distance_xx = distance_xx[np.triu_indices_from(distance_xx, k=1)]
+        distance_yy = distance_yy[np.triu_indices_from(distance_yy, k=1)]
+
     return (
         2 * average(distance_xy) -
         average(distance_xx) -
@@ -40,7 +90,8 @@ def _energy_distance_from_distance_matrices(
     )
 
 
-def energy_distance(x, y, *, average=None, exponent=1):
+def energy_distance(x, y, *, average=None, exponent=1,
+                    estimation_stat=EstimationStatistic.V_STATISTIC):
     """
     Estimator for energy distance.
 
@@ -61,6 +112,12 @@ def energy_distance(x, y, *, average=None, exponent=1):
     average: Callable[[ArrayLike], float]
         A function that will be used to calculate an average of distances.
         This defaults to np.mean.
+    estimation_stat: Union[str, EstimationStatistic]
+        If EstimationStatistic.U_STATISTIC, calculate energy distance using
+        Hoeffding's unbiased U-statistics. Otherwise, use von Mises's biased
+        V-statistics.
+        If this is provided as a string, it will first be converted to
+        an EstimationStatistic enum instance.
 
     Returns
     -------
@@ -111,4 +168,5 @@ def energy_distance(x, y, *, average=None, exponent=1):
         distance_yy=distance_yy,
         distance_xy=distance_xy,
         average=average,
+        estimation_stat=estimation_stat
     )
diff --git a/dcor/homogeneity.py b/dcor/homogeneity.py
index 1095489..5c0fd36 100644
--- a/dcor/homogeneity.py
+++ b/dcor/homogeneity.py
@@ -11,6 +11,7 @@ import numpy as _np
 from . import _energy, _hypothesis
 from . import distances as _distances
 from ._utils import _transform_to_2d
+from ._energy import EstimationStatistic
 
 
 def _energy_test_statistic_coefficient(n, m):
@@ -19,17 +20,26 @@ def _energy_test_statistic_coefficient(n, m):
 
 
 def _energy_test_statistic_from_distance_matrices(
-        distance_xx, distance_yy, distance_xy, n, m, average=None):
+        distance_xx, distance_yy, distance_xy, n, m, average=None,
+        estimation_stat=_energy.EstimationStatistic.V_STATISTIC):
     """Test statistic with precomputed distance matrices."""
     energy_distance = _energy._energy_distance_from_distance_matrices(
         distance_xx=distance_xx, distance_yy=distance_yy,
-        distance_xy=distance_xy, average=average
+        distance_xy=distance_xy, average=average,
+        estimation_stat=estimation_stat
     )
 
     return _energy_test_statistic_coefficient(n, m) * energy_distance
 
 
-def energy_test_statistic(x, y, *, exponent=1, average=None):
+def energy_test_statistic(
+        x,
+        y,
+        *,
+        exponent=1,
+        average=None,
+        estimation_stat=EstimationStatistic.V_STATISTIC
+):
     """
     Homogeneity statistic.
 
@@ -49,6 +59,12 @@ def energy_test_statistic(x, y, *, exponent=1, average=None):
     average: Callable[[ArrayLike], float]
         A function that will be used to calculate an average of distances.
         This defaults to np.mean.
+    estimation_stat: Union[str, EstimationStatistic]
+        If EstimationStatistic.U_STATISTIC, calculate energy distance using
+        Hoeffding's unbiased U-statistics. Otherwise, use von Mises's biased
+        V-statistics.
+        If this is provided as a string, it will first be converted to
+        an EstimationStatistic enum instance.
 
     Returns
     -------
@@ -91,11 +107,12 @@ def energy_test_statistic(x, y, *, exponent=1, average=None):
         y,
         exponent=exponent,
         average=average,
+        estimation_stat=estimation_stat
     )
 
 
 def _energy_test_statistic_multivariate_from_distance_matrix(
-        distance, indexes, sizes, average=None):
+        distance, indexes, sizes, average=None, estimation_stat='v'):
     """Statistic for several random vectors given the distance matrix."""
     energy = 0.0
 
@@ -113,7 +130,9 @@ def _energy_test_statistic_multivariate_from_distance_matrix(
 
             pairwise_energy = _energy_test_statistic_from_distance_matrices(
                 distance_xx=distance_xx, distance_yy=distance_yy,
-                distance_xy=distance_xy, n=n, m=m, average=average)
+                distance_xy=distance_xy, n=n, m=m, average=average,
+                estimation_stat=estimation_stat
+                )
 
             energy += pairwise_energy
 
@@ -126,6 +145,7 @@ def energy_test(
     exponent=1,
     random_state=None,
     average=None,
+    estimation_stat=_energy.EstimationStatistic.V_STATISTIC
 ):
     """
     Test of homogeneity based on the energy distance.
@@ -150,6 +170,12 @@ def energy_test(
     average: Callable[[ArrayLike], float]
         A function that will be used to calculate an average of distances.
         This defaults to np.mean.
+    estimation_stat: Union[str, EstimationStatistic]
+        If EstimationStatistic.U_STATISTIC, calculate energy distance using
+        Hoeffding's unbiased U-statistics. Otherwise, use von Mises's biased
+        V-statistics.
+        If this is provided as a string, it will first be converted to
+        an EstimationStatistic enum instance.
 
     Returns
     -------
@@ -225,7 +251,9 @@ def energy_test(
             distance=distance_matrix,
             indexes=sample_indexes,
             sizes=sample_sizes,
-            average=average)
+            average=average,
+            estimation_stat=estimation_stat
+            )
 
     return _hypothesis._permutation_test_with_sym_matrix(
         sample_distances,


Test patch:
diff --git a/dcor/tests/test_energy.py b/dcor/tests/test_energy.py
new file mode 100644
index 0000000..e7d9571
--- /dev/null
+++ b/dcor/tests/test_energy.py
@@ -0,0 +1,65 @@
+import numpy as np
+from dcor._energy import _energy_distance_from_distance_matrices
+from dcor import EstimationStatistic
+import unittest
+
+
+class TestEnergyTest(unittest.TestCase):
+    """Tests for energy statistics"""
+
+    def test_u_v_statistics(self):
+        """
+        Test that we are correctly applying the provided type of estimation
+        statistics
+
+        Also tests that we can choose a statistic type using either a string
+        or an enum instance
+        """
+        distance_within = np.array([
+            [0, 1, 1],
+            [1, 0, 1],
+            [1, 1, 0],
+        ])
+        distance_between = np.array([
+            [1, 1, 1],
+            [1, 1, 1],
+            [1, 1, 1],
+        ])
+
+        u_dist = _energy_distance_from_distance_matrices(
+            distance_xx=distance_within,
+            distance_yy=distance_within,
+            distance_xy=distance_between,
+            estimation_stat=EstimationStatistic.U_STATISTIC
+        )
+        u_dist_2 = _energy_distance_from_distance_matrices(
+            distance_xx=distance_within,
+            distance_yy=distance_within,
+            distance_xy=distance_between,
+            estimation_stat='U_STATISTIC'
+        )
+        u_dist_3 = _energy_distance_from_distance_matrices(
+            distance_xx=distance_within,
+            distance_yy=distance_within,
+            distance_xy=distance_between,
+            estimation_stat='u_statistic'
+        )
+        v_dist = _energy_distance_from_distance_matrices(
+            distance_xx=distance_within,
+            distance_yy=distance_within,
+            distance_xy=distance_between,
+            estimation_stat='v'
+        )
+
+        # The first 3 invocations should be identical
+        self.assertEqual(u_dist, u_dist_2)
+        self.assertEqual(u_dist, u_dist_3)
+
+        # V-statistics count the 0 terms on the diagonal, so the impact of the
+        # within-sample distance will be smaller, making the overall distance
+        # larger.
+        self.assertGreater(v_dist, u_dist)
+
+        # Also test for exact values
+        self.assertEqual(u_dist, 0)
+        self.assertAlmostEqual(v_dist, 2/3)