Repository: neurodata/hyppo (language: python)

Issue:
Friedman Rafsky PR
@sampan501 

closes #239 


Discussion / hints:

Fix patch:
diff --git a/benchmarks/indep_power_dimension.py b/benchmarks/indep_power_dimension.py
index 8a4bf3986..239552fed 100644
--- a/benchmarks/indep_power_dimension.py
+++ b/benchmarks/indep_power_dimension.py
@@ -56,7 +56,7 @@
 ]
 
 # haven't run for MaxMargin yet
-remove = ["maxmargin"]
+remove = ["maxmargin", "friedman_rafsky"]
 INDEP_TESTS = dict([(k, v) for k, v in INDEP_TESTS.items() if k not in remove])
 
 
diff --git a/benchmarks/indep_power_sampsize.py b/benchmarks/indep_power_sampsize.py
index 0984806ac..ef0ab0a53 100644
--- a/benchmarks/indep_power_sampsize.py
+++ b/benchmarks/indep_power_sampsize.py
@@ -58,7 +58,7 @@
 ]
 
 # these tests only make sense for > 1 dimension data
-remove = ["maxmargin", "kmerf"]
+remove = ["maxmargin", "kmerf", "friedman_rafsky"]
 INDEP_TESTS = dict([(k, v) for k, v in INDEP_TESTS.items() if k not in remove])
 
 
diff --git a/docs/api/index.rst b/docs/api/index.rst
index 3a98e6923..4695a89e9 100644
--- a/docs/api/index.rst
+++ b/docs/api/index.rst
@@ -21,6 +21,7 @@ Independence
    HHG
    CCA
    RV
+   FriedmanRafsky
 
 
 
diff --git a/docs/refs.bib b/docs/refs.bib
index 4d032f678..38f5ee17f 100644
--- a/docs/refs.bib
+++ b/docs/refs.bib
@@ -29,22 +29,7 @@ @article{hardoonCanonicalCorrelationAnalysis2004
   number = {12}
 }
 
-@incollection{hotellingRelationsTwoSets1992,
-  title = {Relations {{Between Two Sets}} of {{Variates}}},
-  booktitle = {Breakthroughs in {{Statistics}}: {{Methodology}} and {{Distribution}}},
-  author = {Hotelling, Harold},
-  editor = {Kotz, Samuel and Johnson, Norman L.},
-  year = {1992},
-  pages = {162--190},
-  publisher = {{Springer}},
-  address = {{New York, NY}},
-  doi = {10.1007/978-1-4612-4380-9_14},
-  abstract = {Concepts of correlation and regression may be applied not only to ordinary one-dimensional variates but also to variates of two or more dimensions. Marksmen side by side firing simultaneous shots at targets, so that the deviations are in part due to independent individual errors and in part to common causes such as wind, provide a familiar introduction to the theory of correlation; but only the correlation of the horizontal components is ordinarily discussed, whereas the complex consisting of horizontal and vertical deviations may be even more interesting. The wind at two places may be compared, using both components of the velocity in each place. A fluctuating vector is thus matched at each moment with another fluctuating vector. The study of individual differences in mental and physical traits calls for a detailed study of the relations between sets of correlated variates. For example the scores on a number of mental tests may be compared with physical measurements on the same persons. The questions then arise of determining the number and nature of the independent relations of mind and body shown by these data to exist, and of extracting from the multiplicity of correlations in the system suitable characterizations of these independent relations. As another example, the inheritance of intelligence in rats might be studied by applying not one but s different mental tests to N mothers and to a daughter of each},
-  isbn = {978-1-4612-4380-9},
-  keywords = {Canonical Correlation,Canonical Variate,Flat Space,Sample Covariance,Vector Correlation},
-  language = {en},
-  series = {Springer {{Series}} in {{Statistics}}}
-}
+./.
 
 @article{knappCanonicalCorrelationAnalysis1978,
   title = {Canonical Correlation Analysis: {{A}} General Parametric Significance-Testing System},
@@ -497,6 +482,22 @@ @article{grettonKernelJointIndependence2016
   archiveprefix = {arXiv},
   keywords = {Math - Statistics Theory, Statistics - Machine Learning},
  }
+ 
+ @article{friedmanMultivariateGeneralizationsoftheWaldWolfowitzandSmirnovTwoSampleTests1979,
+  title = {Multivariate Generalizations of the Wald-Wolfowitz and Smirnov Two-Sample Tests},
+  author = {Friedman, Jerome and Rafsky, Lawrence},
+  year = {1979},
+  month = jul,
+  journal = {Ann. Statist.},
+  volume = {7},
+  number = {4},
+  pages = {697-717},
+  doi = {10.1214/aos/1176344722},
+  abstract = {Multivariate generalizations of the Wald-Wolfowitz runs statistic and the Smirnov maximum deviation statistic for the two-sample problem are presented. They are based on the minimal spanning tree of the pooled sample points. Some null distribution results are derived and a simulation study of power is reported.},
+  copyright = {All rights reserved},
+  keywords = {Wald-Wolfowitz Runs Test,randomness,minimum spanning tree},
+  annotation = {\_eprint: https://doi.org/10.1214/aos/1176344722},
+ }
 
 @inproceedings{hellerMultivariateTestsOfAssociation2016,
   title = {Multivariate Tests of Association Based on Univariate Tests},
@@ -515,4 +516,4 @@ @misc{sasHoeffdingDependenceCoefficient
   author = {SAS},
   howpublished = "\url{https://support.sas.com/documentation/cdl/en/procstat/63104/HTML/default/viewer.htm#procstat_corr_sect016.htm}",
   note = {Accessed: 2021-12-17},
-}
\ No newline at end of file
+}
diff --git a/hyppo/independence/__init__.py b/hyppo/independence/__init__.py
index fe96d0a06..6b97a1749 100644
--- a/hyppo/independence/__init__.py
+++ b/hyppo/independence/__init__.py
@@ -7,6 +7,7 @@
 from .max_margin import MaxMargin
 from .mgc import MGC
 from .rv import RV
+from .friedman_rafsky import FriedmanRafsky
 
 
 __all__ = [s for s in dir()]  # add imported tests to __all__
@@ -21,4 +22,5 @@
     "mgc": MGC,
     "kmerf": KMERF,
     "maxmargin": MaxMargin,
+    "friedman_rafsky": FriedmanRafsky,
 }
diff --git a/hyppo/independence/friedman_rafsky.py b/hyppo/independence/friedman_rafsky.py
new file mode 100644
index 000000000..c9d491c37
--- /dev/null
+++ b/hyppo/independence/friedman_rafsky.py
@@ -0,0 +1,223 @@
+import random
+from numba import jit
+import numpy as np
+
+from .base import IndependenceTest, IndependenceTestOutput
+from ..tools import perm_test
+
+
+class FriedmanRafsky(IndependenceTest):
+    r"""
+    Friedman-Rafksy (FR) test statistic and p-value.
+    This is a multivariate extension of the Wald-Wolfowitz
+    runs test for randomness. The normal concept of a 'run'
+    is replaced by a minimum spanning tree (MST) calculated between
+    the points in respective data sets with edge weights defined
+    as the Euclidean distance between two such points. After MST
+    has been determined, all edges such that both corresponding
+    nodes do not belong to the same class are severed and the
+    number of independent resulting trees is counted. This test is
+    consistent against similar tests.
+
+    Notes
+    -----
+    The statistic can be derived as follows
+    :footcite:p:`friedmanMultivariateGeneralizationsoftheWaldWolfowitzandSmirnovTwoSampleTests1979`
+
+    Let :math:`x` be a combined sample of :math:`(n, p)` and :math:`(m, p)`
+    samples of random variables :math:`X` and let :math:`y` be a :math:`(n+m, 1)`
+    array of labels :math:`Y`. We can then create a graph such that each point in
+    :math:`X` is connected to each other point in :math:`X` by an edge weighted by
+    the euclidean distance inbetween those points. The minimum spanning tree is then
+    calculated and all edges such that the labels in :math:`Y` are not from the same
+    class are removed. The number of independent graphs is then summed to determine
+    the uncorrected statistic for the test.
+
+    The p-value and null distribution for the corrected statistic are calculated via
+    a permutation test using :meth:`hyppo.tools.perm_test`.
+    """
+
+    def __init__(self, **kwargs):
+
+        IndependenceTest.__init__(self, **kwargs)
+
+    def _num_runs(self, labels, MST_connections):
+        r"""
+        Helper function to determine number of independent
+        'runs' from MST connections.
+
+        Parameters
+        ----------
+        labels : ndarry of float
+            Lables corresponding to respective classes of samples.
+        MST_connections: list of int
+            List containing pairs of points connected in final MST.
+
+        Returns
+        -------
+        run_count : int
+            Number of runs after severing all such edges with nodes of
+            differing class labels.
+        """
+        run_count = 1
+
+        for x in MST_connections:
+            if labels[x[0]] != labels[x[1]]:
+                run_count += 1
+
+        return run_count
+
+    def statistic(self, x, y):
+        r"""
+        Helper function that calculates the Friedman Rafksy test statistic.
+
+        Parameters
+        ----------
+        x,y : ndarray of float
+            Input data matrices. ``x`` and ``y`` must have the same number of
+            rows. That is, the shapes must be ``(n, p)`` and ``(n, 1)`` where
+            `n` is the number of combined samples and `p` is the number of
+            dimensions. ``y`` is the array of labels corresponding to the two
+            samples, respectively.
+
+        Returns
+        -------
+        stat : float
+            The computed Friedman Rafsky statistic.
+        """
+        x = np.transpose(x)
+        labels = np.transpose(y)
+
+        MST_connections = MST(x, labels)
+        stat = self._num_runs(labels, MST_connections)
+
+        return stat
+
+    def test(
+        self,
+        x,
+        y,
+        reps=1000,
+        workers=1,
+        random_state=None,
+    ):
+        r"""
+        Calculates the Friedman Rafsky test statistic and p-value.
+
+        Parameters
+        ----------
+        x,y : ndarray of float
+            Input data matrices. ``x`` and ``y`` must have the same number of
+            rows. That is, the shapes must be ``(n, p)`` and ``(n, 1)`` where
+            `n` is the number of combined samples and `p` is the number of
+            dimensions. ``y`` is the array of labels corresponding to the two
+            samples, respectively.
+        reps : int, default: 1000
+            The number of replications used to estimate the null distribution
+            when using the permutation test used to calculate the p-value.
+        workers : int, default: 1
+            The number of cores to parallelize the p-value computation over.
+            Supply ``-1`` to use all cores available to the Process.
+        random_state : int, default: None
+            The random_state for permutation testing to be fixed for
+            reproducibility.
+
+        Returns
+        -------
+        stat : float
+            The computed Friedman Rafsky statistic.
+        pvalue : float
+            The computed Friedman Rafsky p-value.
+        """
+
+        stat, pvalue, null_dist = perm_test(
+            self.statistic,
+            x,
+            y,
+            reps,
+            workers,
+            is_distsim=False,
+            random_state=random_state,
+        )
+        stat = (stat - np.mean(null_dist)) / np.std(null_dist)
+        self.stat = stat
+
+        return IndependenceTestOutput(stat, pvalue)
+
+
+@jit(nopython=True, cache=True)  # pragma: no cover
+def prim(weight_mat, labels):
+    r"""
+    Helper function to read weighted matrix input and compute minimum
+    spanning tree via Prim's algorithm.
+
+    Parameters
+    ----------
+    weight_mat : ndarry of float
+        Weighted connection matrix.
+    labels : ndarry of int
+        Lables corresponding to respective classes of samples.
+
+    Returns
+    -------
+    MST_connections : list of int
+        List of pairs of nodes connected in final MST.
+    """
+    INF = 9999999
+    V = len(labels)
+    selected = np.zeros(len(labels))
+    no_edge = 0
+    selected[0] = True
+    MST_connections = []
+
+    while no_edge < V - 1:
+        minimum = INF
+        x = 0
+        y = 0
+        for i in range(V):
+            if selected[i]:
+                for j in range(V):
+                    if (not selected[j]) and weight_mat[i][j]:
+                        # not in selected and there is an edge
+                        if minimum > weight_mat[i][j]:
+                            minimum = weight_mat[i][j]
+                            x = i
+                            y = j
+        MST_connections.append([x, y])
+        selected[y] = True
+        no_edge += 1
+
+    return MST_connections
+
+
+@jit(nopython=True, cache=True)
+def MST(x, labels):  # pragma: no cover
+    r"""
+    Helper function to read input data and calculate Euclidean distance
+    between each possible pair of points before finding MST.
+
+    Parameters
+    ----------
+    x : ndarry of float
+        Dataset such that each column corresponds to a point of data.
+    labels : ndarry of int
+        Lables corresponding to respective classes of samples.
+
+    Returns
+    -------
+    MST_connections : list
+        List of pairs of nodes connected in final MST.
+    """
+
+    G = np.zeros((len(x[0]), len(x[0])))
+
+    for i in range(len(x[0])):
+
+        for j in range(i + 1, len(x[0])):
+            weight = np.linalg.norm(x[:, i] - x[:, j])
+            G[i][j] = weight
+            G[j][i] = weight
+
+    MST_connections = prim(G, labels)
+
+    return MST_connections
diff --git a/tutorials/images/fr_full_tree.png b/tutorials/images/fr_full_tree.png
new file mode 100644
index 000000000..6313abdfb
Binary files /dev/null and b/tutorials/images/fr_full_tree.png differ
diff --git a/tutorials/images/fr_runs_tree.png b/tutorials/images/fr_runs_tree.png
new file mode 100644
index 000000000..4164c0e29
Binary files /dev/null and b/tutorials/images/fr_runs_tree.png differ
diff --git a/tutorials/independence.py b/tutorials/independence.py
index ce7089e53..1f0e04bbf 100644
--- a/tutorials/independence.py
+++ b/tutorials/independence.py
@@ -78,7 +78,7 @@
 #    :Cons: - Very slow (computationally complex)
 #           - Test statistic not very interpretable, not between (-1, 1)
 #
-# This test runs like :ref:`any other test<general indep>` and can be implemented 
+# This test runs like :ref:`any other test<general indep>` and can be implemented
 # as below:
 
 import timeit
@@ -95,7 +95,7 @@
 t_fast_hhg = timeit.Timer(stmt="HHG().test(x, y, auto=True)", setup=setup_code)
 
 hhg_time = np.array(t_hhg.timeit(number=1))  # original HHG
-fast_hhg_time  = np.array(t_fast_hhg.timeit(number=5)) / 5  # fast HHG
+fast_hhg_time = np.array(t_fast_hhg.timeit(number=5)) / 5  # fast HHG
 
 print(u"Original HHG time: {0:.3g}s".format(hhg_time))
 print(u"Fast HHG time: {0:.3g}s".format(fast_hhg_time))
@@ -334,5 +334,44 @@
 # the relevant section of reference documentation in :mod:`hyppo.independence` for more
 # information.
 # These tests runs like :ref:`any other test<general indep>`.
+
+########################################################################################
+# Friedman Rafsky Test for Randomness
+# --------------------------------------------
+#
+# This notebook will introduce the usage of the Friedman Rafsky test, a multivariate
+# extension of the Wald-Wolfowitz runs test to test for randomness between two multivariate
+# samples. More specifically, the function tests whether two multivariate samples were
+# independently drawn from the same distribution.
+#
+# The question proposed in 'Multivariate Generalizations of the Wald-Wolfowitz and Smirnov Two-Sample Tests'
+# is that of how to extend the univariate Wald-Wolfowitz runs test to a multivariate setting.
+#
+# The univariate Wald-Wolfowitz runs test is a non-parametric statistical test that checks a randomness hypothesis
+# for a two-valued data sequence. More specifically, it can be used to test the hypothesis that the elements of a sequence
+# are mutually independent. For a data sequence with identifiers of two groups,
+# say: :math:`X , Y` we begin by sorting the combined data set
+# :math: `W` in numerical ascending order. The number of runs is then defined by the number
+# of maximal, non-empty segments of the sequence consisting of adjacent and equal elements.
+# So if we designate every :math:`X = +,  Y = -` an example assortment of the 15 element long sequence
+# of both sets could be given as :math:`+++----++---+++` which contains 5 runs, 3 of which positive and 2 of which negative.
+# By randomly permuting the labels of our data a large number of times, we can compare our true number of runs to that of the
+# random permutations to determine the test statistic and p-value.
+#
+# For extension to the multivariate case, we determine a neighbor not by proximity in numerical order,
+# but instead by euclidean distance between each point. The data is then 'sorted' via a calculation of
+# the minimum spanning tree in which each point is connected to each other point with edge weight equal
+# to that of the euclidean distance between the two points. The number of runs is then determined by severing
+# each edge of the MST for which the points connected do not belong to the same family. Similarly to the univariate
+# case, the labels associated with each node are then permuted randomly a large number of times, we can
+# compare our number of true runs to the random distribution of runs to determine the multivariate test
+# statistic and p-value.
+#
+# As such, we see that this set of data contains 3 such runs, and again by randomizing the
+# labels of each data point we can determine the test statistic and p-value for the hypothesis
+# that to determine the randomness of our combined sample.
+# It's worth noting that X and Y have the same number of samples, just that they posess the same number of multivariate
+# features. Lastly, labels for X and Y need not be 0 and 1, they just need be consistent across samples.
+# These tests runs like :ref:`any other test<general indep>`.
 #
 # .. _[1]: https://link.springer.com/article/10.1007/s10182-020-00378-1


Test patch:
diff --git a/hyppo/independence/tests/test_friedman_rafsky.py b/hyppo/independence/tests/test_friedman_rafsky.py
new file mode 100644
index 000000000..567339f72
--- /dev/null
+++ b/hyppo/independence/tests/test_friedman_rafsky.py
@@ -0,0 +1,37 @@
+import numpy as np
+import pytest
+from numpy.testing import assert_almost_equal, assert_raises, assert_warns
+
+from ...tools import linear, power
+from .. import FriedmanRafsky
+
+
+class TestFriedmanRafskyStat:
+    @pytest.mark.parametrize("n", [100])
+    @pytest.mark.parametrize("num_runs", [48])
+    @pytest.mark.parametrize("obs_stat", [-0.527])
+    @pytest.mark.parametrize("obs_pvalue", [0.748])
+    def test_linear_oned(self, n, num_runs, obs_stat, obs_pvalue):
+        np.random.seed(123456789)
+        x, y = linear(n, 1)
+        num_rows, num_cols = x.shape
+        y = np.random.choice([0, 1], num_rows, p=[0.5, 0.5])
+        y = np.transpose(y)
+        stat1, pvalue1 = FriedmanRafsky().test(x, y)
+        stat2 = FriedmanRafsky().statistic(x, y)
+
+        assert_almost_equal(stat1, obs_stat, decimal=2)
+        assert_almost_equal(stat2, num_runs, decimal=2)
+        assert_almost_equal(pvalue1, obs_pvalue, decimal=2)
+
+    @pytest.mark.parametrize("n", [100, 200])
+    def test_rep(self, n):
+        x, y = linear(n, 1)
+        num_rows, num_cols = x.shape
+        y = np.random.choice([0, 1], num_rows, p=[0.5, 0.5])
+        y = np.transpose(y)
+        stat1, pvalue1 = FriedmanRafsky().test(x, y, random_state=2)
+        stat2, pvalue2 = FriedmanRafsky().test(x, y, random_state=2)
+
+        assert stat1 == stat2
+        assert pvalue1 == pvalue2