Repository: statsmodels/statsmodels
Pull request #6401 (2020-01-10 07:45:40)

ENH: Distance dependence measures
- [x] tests added / passed. 
- [x] code/documentation is well formatted.  
- [x] properly formatted commit message. See 
      [NumPy's guide](https://docs.scipy.org/doc/numpy-1.15.1/dev/gitwash/development_workflow.html#writing-the-commit-message). 

An implementation of Distance dependence measures and the Distance Covariance (dCov) test, as introduced by Szekely, G.J., Rizzo, M.L., and Bakirov, N.K. (2007) “Measuring and testing dependence by correlation of distances”. Annals of Statistics, Vol. 35 No. 6, pp. 2769-2794.

Any and all feedback is most welcome and I'd be happy to make any changes if need be.




pep8speaks: Hello @roncho12! Thanks for updating this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:



* In the file [`statsmodels/stats/tests/test_dist_dependant_measures.py`](https://github.com/statsmodels/statsmodels/blob/e25965fe04d084d1e9fb4dca6f7924573bb4e84a/statsmodels/stats/tests/test_dist_dependant_measures.py):

> [Line 2:1](https://github.com/statsmodels/statsmodels/blob/e25965fe04d084d1e9fb4dca6f7924573bb4e84a/statsmodels/stats/tests/test_dist_dependant_measures.py#L2): [F401](https://duckduckgo.com/?q=pep8%20F401) 'numpy.testing.assert_allclose' imported but unused


codecov[bot]: # [Codecov](https://codecov.io/gh/statsmodels/statsmodels/pull/6401?src=pr&el=h1) Report
> Merging [#6401](https://codecov.io/gh/statsmodels/statsmodels/pull/6401?src=pr&el=desc) into [master](https://codecov.io/gh/statsmodels/statsmodels/commit/da43833e9a66abc0c7c74577c88cdd4152aed354?src=pr&el=desc) will **increase** coverage by `0.01%`.
> The diff coverage is `97.64%`.

[![Impacted file tree graph](https://codecov.io/gh/statsmodels/statsmodels/pull/6401/graphs/tree.svg?width=650&token=Og706qJV5S&height=150&src=pr)](https://codecov.io/gh/statsmodels/statsmodels/pull/6401?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master    #6401      +/-   ##
==========================================
+ Coverage   85.04%   85.05%   +0.01%     
==========================================
  Files         640      642       +2     
  Lines      103193   103363     +170     
  Branches    11277    11284       +7     
==========================================
+ Hits        87759    87916     +157     
- Misses      12957    12962       +5     
- Partials     2477     2485       +8
```


| [Impacted Files](https://codecov.io/gh/statsmodels/statsmodels/pull/6401?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [...models/stats/tests/test\_dist\_dependant\_measures.py](https://codecov.io/gh/statsmodels/statsmodels/pull/6401/diff?src=pr&el=tree#diff-c3RhdHNtb2RlbHMvc3RhdHMvdGVzdHMvdGVzdF9kaXN0X2RlcGVuZGFudF9tZWFzdXJlcy5weQ==) | `100% <100%> (ø)` | |
| [statsmodels/stats/dist\_dependence\_measures.py](https://codecov.io/gh/statsmodels/statsmodels/pull/6401/diff?src=pr&el=tree#diff-c3RhdHNtb2RlbHMvc3RhdHMvZGlzdF9kZXBlbmRlbmNlX21lYXN1cmVzLnB5) | `94.8% <94.8%> (ø)` | |
| [statsmodels/tests/test\_package.py](https://codecov.io/gh/statsmodels/statsmodels/pull/6401/diff?src=pr&el=tree#diff-c3RhdHNtb2RlbHMvdGVzdHMvdGVzdF9wYWNrYWdlLnB5) | `70.58% <0%> (-29.42%)` | :arrow_down: |
| [statsmodels/tools/docstring.py](https://codecov.io/gh/statsmodels/statsmodels/pull/6401/diff?src=pr&el=tree#diff-c3RhdHNtb2RlbHMvdG9vbHMvZG9jc3RyaW5nLnB5) | `77.03% <0%> (-0.96%)` | :arrow_down: |
| [statsmodels/regression/mixed\_linear\_model.py](https://codecov.io/gh/statsmodels/statsmodels/pull/6401/diff?src=pr&el=tree#diff-c3RhdHNtb2RlbHMvcmVncmVzc2lvbi9taXhlZF9saW5lYXJfbW9kZWwucHk=) | `89.84% <0%> (-0.76%)` | :arrow_down: |
| [statsmodels/tsa/stattools.py](https://codecov.io/gh/statsmodels/statsmodels/pull/6401/diff?src=pr&el=tree#diff-c3RhdHNtb2RlbHMvdHNhL3N0YXR0b29scy5weQ==) | `90.05% <0%> (-0.75%)` | :arrow_down: |
| [statsmodels/tools/grouputils.py](https://codecov.io/gh/statsmodels/statsmodels/pull/6401/diff?src=pr&el=tree#diff-c3RhdHNtb2RlbHMvdG9vbHMvZ3JvdXB1dGlscy5weQ==) | `69.46% <0%> (-0.45%)` | :arrow_down: |
| [statsmodels/datasets/utils.py](https://codecov.io/gh/statsmodels/statsmodels/pull/6401/diff?src=pr&el=tree#diff-c3RhdHNtb2RlbHMvZGF0YXNldHMvdXRpbHMucHk=) | `84.66% <0%> (+2.45%)` | :arrow_up: |
| ... and [1 more](https://codecov.io/gh/statsmodels/statsmodels/pull/6401/diff?src=pr&el=tree-more) | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/statsmodels/statsmodels/pull/6401?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/statsmodels/statsmodels/pull/6401?src=pr&el=footer). Last update [da43833...fa01ccd](https://codecov.io/gh/statsmodels/statsmodels/pull/6401?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

coveralls: 
[![Coverage Status](https://coveralls.io/builds/28025308/badge)](https://coveralls.io/builds/28025308)

Coverage increased (+0.02%) to 87.46% when pulling **c48e9751b0f08d1bd72459a7e76635bf61ffb04f on roncho12:distance_dependence_measurements** into **da43833e9a66abc0c7c74577c88cdd4152aed354 on statsmodels:master**.

roncho12: @bashtage Thanks!
I've changed the names.
You're correct, i used the same naming as the paper, but it does make sense to make them more explicit.

Patch:
diff --git a/codecov.yml b/.codecov.yml
similarity index 54%
rename from codecov.yml
rename to .codecov.yml
index 7d38087e5..86c981126 100644
--- a/codecov.yml
+++ b/.codecov.yml
@@ -1,16 +1,11 @@
 codecov:
-  branch: master
   notify:
     require_ci_to_pass: no
     after_n_builds: 1
 
-comment: off
-
 coverage:
   status:
     project:
       default:
-        target: '82'
-    patch:
-      default:
-        target: '50'
+        # basic
+        threshold: 0.05
diff --git a/docs/source/contrasts.rst b/docs/source/contrasts.rst
index cd6ec5705..c349f5460 100644
--- a/docs/source/contrasts.rst
+++ b/docs/source/contrasts.rst
@@ -3,7 +3,7 @@
 Patsy: Contrast Coding Systems for categorical variables
 ===========================================================
 
-.. note:: This document is based on `this excellent resource from UCLA <https://stats.idre.ucla.edu/r/library/r-library-contrast-coding-systems-for-categorical-variables/>`__.
+.. note:: This document is based heavily on `this excellent resource from UCLA <http://www.ats.ucla.edu/stat/r/library/contrast_coding.htm>`__.
 
 A categorical variable of K categories, or levels, usually enters a regression as a sequence of K-1 dummy variables. This amounts to a linear hypothesis on the level means. That is, each test statistic for these variables amounts to testing whether the mean for that level is statistically significantly different from the mean of the base category. This dummy coding is called Treatment coding in R parlance, and we will follow this convention. There are, however, different coding methods that amount to different sets of linear hypotheses.
 
diff --git a/docs/source/release/version0.11.rst b/docs/source/release/version0.11.rst
index 9c1000346..a284fa90a 100644
--- a/docs/source/release/version0.11.rst
+++ b/docs/source/release/version0.11.rst
@@ -42,6 +42,10 @@ Statistics
 Add the Oaxaca-Blinder decomposition (:class:`~statsmodels.stats.oaxaca.OaxacaBlinder`) that
 decomposes the difference in group means into with and between components.
 
+Add the Distance dependence measures statistics
+(:func:`~statsmodels.stats.dist_dependence_measures.distance_statistics`) and the Distance Covariance
+test (:func:`~statsmodels.stats.dist_dependence_measures.distance_covariance_test`).
+
 Statespace Models
 -----------------
 
diff --git a/docs/source/stats.rst b/docs/source/stats.rst
index 89dc74b07..f374c8e40 100644
--- a/docs/source/stats.rst
+++ b/docs/source/stats.rst
@@ -516,4 +516,25 @@ This function attempts to port the functionality of the oaxaca command in STATA
    :toctree: generated/
 
    OaxacaBlinder
-   OaxacaResults
\ No newline at end of file
+   OaxacaResults
+
+
+Distance Dependence Measures
+----------------------------
+
+Distance dependence measures and the Distance Covariance (dCov) test.
+
+.. module:: statsmodels.stats.dist_dependence_measures
+   :synopsis: Distance Dependence Measures
+
+.. currentmodule:: statsmodels.stats.dist_dependence_measures
+
+.. autosummary::
+   :toctree: generated/
+
+   distance_covariance_test
+   distance_statistics
+   distance_correlation
+   distance_covariance
+   distance_variance
+
diff --git a/statsmodels/stats/dist_dependence_measures.py b/statsmodels/stats/dist_dependence_measures.py
new file mode 100644
index 000000000..3a9b40579
--- /dev/null
+++ b/statsmodels/stats/dist_dependence_measures.py
@@ -0,0 +1,507 @@
+# -*- coding: utf-8 -*-
+""" Distance dependence measure and the dCov test.
+
+Implementation of Székely et al. (2007) calculation of distance
+dependence statistics, including the Distance covariance (dCov) test
+for independence of random vectors of arbitrary length.
+
+Author: Ron Itzikovitch
+
+References
+----------
+.. Szekely, G.J., Rizzo, M.L., and Bakirov, N.K. (2007)
+   "Measuring and testing dependence by correlation of distances".
+   Annals of Statistics, Vol. 35 No. 6, pp. 2769-2794.
+
+"""
+import numpy as np
+import warnings
+from collections import namedtuple
+
+from scipy.spatial.distance import pdist, squareform
+from scipy.stats import norm
+
+
+DistDependStat = namedtuple(
+    "DistDependStat",
+    ["test_statistic", "distance_correlation",
+     "distance_covariance", "dvar_x", "dvar_y", "S"],
+)
+
+
+def distance_covariance_test(x, y, B=None, method="auto"):
+    r"""The Distance Covariance (dCov) test
+
+    Apply the Distance Covariance (dCov) test of independence to `x` and `y`.
+    This test was introduced in [1]_, and is based on the distance covariance
+    statistic. The test is applicable to random vectors of arbitrary length
+    (see the notes section for more details).
+
+    Parameters
+    ----------
+    x : array_like, 1-D or 2-D
+        If `x` is 1-D than it is assumed to be a vector of observations of a
+        single random variable. If `x` is 2-D than the rows should be
+        observations and the columns are treated as the components of a
+        random vector, i.e., each column represents a different component of
+        the random vector `x`.
+    y : array_like, 1-D or 2-D
+        Same as `x`, but only the number of observation has to match that of
+        `x`. If `y` is 2-D note that the number of columns of `y` (i.e., the
+        number of components in the random vector) does not need to match
+        the number of columns in `x`.
+    B : int, optional, default=`None`
+        The number of iterations to perform when evaluating the null
+        distribution of the test statistic when the `emp` method is
+        applied (see below). if `B` is `None` than as in [1]_ we set
+        `B` to be ``B = 200 + 5000/n``, where `n` is the number of
+        observations.
+    method : {'auto', 'emp', 'asym'}, optional, default=auto
+        The method by which to obtain the p-value for the test.
+
+        - `auto` : Default method. The number of observations will be used to
+          determine the method.
+        - `emp` : Empirical evaluation of the p-value using permutations of
+          the rows of `y` to obtain the null distribution.
+        - `asym` : An asymptotic approximation of the distribution of the test
+          statistic is used to find the p-value.
+
+    Returns
+    -------
+    test_statistic : float
+        The value of the test statistic used in the test.
+    pval : float
+        The p-value.
+    chosen_method : str
+        The method that was used to obtain the p-value. Mostly relevant when
+        the function is called with `method='auto'`.
+
+    Notes
+    -----
+    The test applies to random vectors of arbitrary dimensions, i.e., `x`
+    can be a 1-D vector of observations for a single random variable while
+    `y` can be a `k` by `n` 2-D array (where `k > 1`). In other words, it
+    is also possible for `x` and `y` to both be 2-D arrays and have the
+    same number of rows (observations) while differing in the number of
+    columns.
+
+    As noted in [1]_ the statistics are sensitive to all types of departures
+    from independence, including nonlinear or nonmonotone dependence
+    structure.
+
+    References
+    ----------
+    .. [1] Szekely, G.J., Rizzo, M.L., and Bakirov, N.K. (2007)
+       "Measuring and testing by correlation of distances".
+       Annals of Statistics, Vol. 35 No. 6, pp. 2769-2794.
+
+    Examples
+    --------
+    >>> from statsmodels.stats.dist_dependence_measures import
+    ... distance_covariance_test
+    >>> data = np.random.rand(1000, 10)
+    >>> x, y = data[:, :3], data[:, 3:]
+    >>> x.shape
+    (1000, 3)
+    >>> y.shape
+    (1000, 7)
+    >>> distance_covariance_test(x, y)
+    (1.0426404792714983, 0.2971148340813543, 'asym')
+    # (test_statistic, pval, chosen_method)
+
+    """
+    x, y = _validate_and_tranform_x_and_y(x, y)
+
+    n = x.shape[0]
+    stats = distance_statistics(x, y)
+
+    if method == "auto" and n <= 500 or method == "emp":
+        chosen_method = "emp"
+        test_statistic, pval = _empirical_pvalue(x, y, B, n, stats)
+
+    elif method == "auto" and n > 500 or method == "asym":
+        chosen_method = "asym"
+        test_statistic, pval = _asymptotic_pvalue(stats)
+
+    else:
+        raise ValueError("Unknown 'method' parameter: {}".format(method))
+
+    # In case we got an extreme p-value (0 or 1) when using the empirical
+    # distribution of the test statistic under the null, we fall back
+    # to the asymptotic approximation.
+    if chosen_method == "emp" and pval in [0, 1]:
+        msg = (
+            "p-value was {} when using the empirical method. ".format(pval)
+            + "The asymptotic approximation will be used instead"
+        )
+        warnings.warn(msg)
+        _, pval = _asymptotic_pvalue(stats)
+
+    return test_statistic, pval, chosen_method
+
+
+def _validate_and_tranform_x_and_y(x, y):
+    r"""Ensure `x` and `y` have proper shape and transform/reshape them if
+    required.
+
+    Parameters
+    ----------
+    x : array_like, 1-D or 2-D
+        If `x` is 1-D than it is assumed to be a vector of observations of a
+        single random variable. If `x` is 2-D than the rows should be
+        observations and the columns are treated as the components of a
+        random vector, i.e., each column represents a different component of
+        the random vector `x`.
+    y : array_like, 1-D or 2-D
+        Same as `x`, but only the number of observation has to match that of
+        `x`. If `y` is 2-D note that the number of columns of `y` (i.e., the
+        number of components in the random vector) does not need to match
+        the number of columns in `x`.
+
+    Returns
+    -------
+    x : array_like, 1-D or 2-D
+    y : array_like, 1-D or 2-D
+
+    Raises
+    ------
+    ValueError
+        If `x` and `y` have a different number of observations.
+
+    """
+    x = np.asanyarray(x)
+    y = np.asanyarray(y)
+
+    if x.shape[0] != y.shape[0]:
+        raise ValueError(
+            "x and y must have the same number of observations (rows)."
+        )
+
+    if len(x.shape) == 1:
+        x = x.reshape((x.shape[0], 1))
+
+    if len(y.shape) == 1:
+        y = y.reshape((y.shape[0], 1))
+
+    return x, y
+
+
+def _empirical_pvalue(x, y, B, n, stats):
+    r"""Calculate the empirical p-value based on permutations of `y`'s rows
+
+    Parameters
+    ----------
+    x : array_like, 1-D or 2-D
+        If `x` is 1-D than it is assumed to be a vector of observations of a
+        single random variable. If `x` is 2-D than the rows should be
+        observations and the columns are treated as the components of a
+        random vector, i.e., each column represents a different component of
+        the random vector `x`.
+    y : array_like, 1-D or 2-D
+        Same as `x`, but only the number of observation has to match that of
+        `x`. If `y` is 2-D note that the number of columns of `y` (i.e., the
+        number of components in the random vector) does not need to match
+        the number of columns in `x`.
+    B : int
+        The number of iterations when evaluating the null distribution.
+    n : Number of observations found in each of `x` and `y`.
+    stats: namedtuple
+        The result obtained from calling ``distance_statistics(x, y)``.
+
+    Returns
+    -------
+    test_statistic : float
+        The empirical test statistic.
+    pval: float
+        The empirical p-value.
+
+    """
+    B = int(B) if B else int(np.floor(200 + 5000 / n))
+    empirical_dist = _get_test_statistic_distribution(x, y, B)
+    pval = 1 - np.searchsorted(
+        sorted(empirical_dist), stats.test_statistic
+    ) / len(empirical_dist)
+    test_statistic = stats.test_statistic
+
+    return test_statistic, pval
+
+
+def _asymptotic_pvalue(stats):
+    r"""Calculate the p-value based on an approximation of the distribution of
+    the test statistic under the null.
+
+    Parameters
+    ----------
+    stats: namedtuple
+        The result obtained from calling ``distance_statistics(x, y)``.
+
+    Returns
+    -------
+    test_statistic : float
+        The test statistic.
+    pval: float
+        The asymptotic p-value.
+
+    """
+    test_statistic = np.sqrt(stats.test_statistic / stats.S)
+    pval = (1 - norm.cdf(test_statistic)) * 2
+
+    return test_statistic, pval
+
+
+def _get_test_statistic_distribution(x, y, B):
+    r"""
+    Parameters
+    ----------
+    x : array_like, 1-D or 2-D
+        If `x` is 1-D than it is assumed to be a vector of observations of a
+        single random variable. If `x` is 2-D than the rows should be
+        observations and the columns are treated as the components of a
+        random vector, i.e., each column represents a different component of
+        the random vector `x`.
+    y : array_like, 1-D or 2-D
+        Same as `x`, but only the number of observation has to match that of
+        `x`. If `y` is 2-D note that the number of columns of `y` (i.e., the
+        number of components in the random vector) does not need to match
+        the number of columns in `x`.
+    B : int
+        The number of iterations to perform when evaluating the null
+        distribution.
+
+    Returns
+    -------
+    emp_dist : array_like
+        The empirical distribution of the test statistic.
+
+    """
+    y = y.copy()
+    emp_dist = np.zeros(B)
+    x_dist = squareform(pdist(x, "euclidean"))
+
+    for i in range(B):
+        np.random.shuffle(y)
+        emp_dist[i] = distance_statistics(x, y, x_dist=x_dist).test_statistic
+
+    return emp_dist
+
+
+def distance_statistics(x, y, x_dist=None, y_dist=None):
+    r"""Calculate various distance dependence statistics.
+
+    Calculate several distance dependence statistics as described in [1]_.
+
+    Parameters
+    ----------
+    x : array_like, 1-D or 2-D
+        If `x` is 1-D than it is assumed to be a vector of observations of a
+        single random variable. If `x` is 2-D than the rows should be
+        observations and the columns are treated as the components of a
+        random vector, i.e., each column represents a different component of
+        the random vector `x`.
+    y : array_like, 1-D or 2-D
+        Same as `x`, but only the number of observation has to match that of
+        `x`. If `y` is 2-D note that the number of columns of `y` (i.e., the
+        number of components in the random vector) does not need to match
+        the number of columns in `x`.
+    x_dist : array_like, 2-D, optional
+        A square 2-D array_like object whose values are the euclidean
+        distances between `x`'s rows.
+    y_dist : array_like, 2-D, optional
+        A square 2-D array_like object whose values are the euclidean
+        distances between `y`'s rows.
+
+    Returns
+    -------
+    namedtuple
+        A named tuple of distance dependance statistics (DistDependStat) with
+        the following values:
+
+        - test_statistic : float - The "basic" test statistic (i.e., the one
+          used when the `emp` method is chosen when calling
+          ``distance_covariance_test()``
+        - distance_correlation : float - The distance correlation
+          between `x` and `y`.
+        - distance_covariance : float - The distance covariance of
+          `x` and `y`.
+        - dvar_x : float - The distance variance of `x`.
+        - dvar_y : float - The distance variance of `y`.
+        - S : float - The mean of the euclidean distances in `x` multiplied
+          by those of `y`. Mostly used internally.
+
+    References
+    ----------
+    .. [1] Szekely, G.J., Rizzo, M.L., and Bakirov, N.K. (2007)
+       "Measuring and testing dependence by correlation of distances".
+       Annals of Statistics, Vol. 35 No. 6, pp. 2769-2794.
+
+    Examples
+    --------
+
+    >>> from statsmodels.stats.dist_dependence_measures import
+    ... distance_statistics
+    >>> distance_statistics(np.random.random(1000), np.random.random(1000))
+    DistDependStat(test_statistic=0.07948284320205831,
+    distance_correlation=0.04269511890990793,
+    distance_covariance=0.008915315092696293,
+    dvar_x=0.20719027438266704, dvar_y=0.21044934264957588,
+    S=0.10892061635588891)
+
+    """
+    x, y = _validate_and_tranform_x_and_y(x, y)
+
+    n = x.shape[0]
+
+    a = x_dist if x_dist is not None else squareform(pdist(x, "euclidean"))
+    b = y_dist if y_dist is not None else squareform(pdist(y, "euclidean"))
+
+    a_row_means = a.mean(axis=0, keepdims=True)
+    b_row_means = b.mean(axis=0, keepdims=True)
+    a_col_means = a.mean(axis=1, keepdims=True)
+    b_col_means = b.mean(axis=1, keepdims=True)
+    a_mean = a.mean()
+    b_mean = b.mean()
+
+    A = a - a_row_means - a_col_means + a_mean
+    B = b - b_row_means - b_col_means + b_mean
+
+    S = a_mean * b_mean
+    dcov = np.sqrt(np.multiply(A, B).mean())
+    dvar_x = np.sqrt(np.multiply(A, A).mean())
+    dvar_y = np.sqrt(np.multiply(B, B).mean())
+    dcor = dcov / np.sqrt(dvar_x * dvar_y)
+
+    test_statistic = n * dcov ** 2
+
+    return DistDependStat(
+        test_statistic=test_statistic,
+        distance_correlation=dcor,
+        distance_covariance=dcov,
+        dvar_x=dvar_x,
+        dvar_y=dvar_y,
+        S=S,
+    )
+
+
+def distance_covariance(x, y):
+    r"""Distance covariance.
+
+    Calculate the empirical distance covariance as described in [1]_.
+
+    Parameters
+    ----------
+    x : array_like, 1-D or 2-D
+        If `x` is 1-D than it is assumed to be a vector of observations of a
+        single random variable. If `x` is 2-D than the rows should be
+        observations and the columns are treated as the components of a
+        random vector, i.e., each column represents a different component of
+        the random vector `x`.
+    y : array_like, 1-D or 2-D
+        Same as `x`, but only the number of observation has to match that of
+        `x`. If `y` is 2-D note that the number of columns of `y` (i.e., the
+        number of components in the random vector) does not need to match
+        the number of columns in `x`.
+
+    Returns
+    -------
+    float
+        The empirical distance covariance between `x` and `y`.
+
+    References
+    ----------
+    .. [1] Szekely, G.J., Rizzo, M.L., and Bakirov, N.K. (2007)
+       "Measuring and testing dependence by correlation of distances".
+       Annals of Statistics, Vol. 35 No. 6, pp. 2769-2794.
+
+    Examples
+    --------
+
+    >>> from statsmodels.stats.dist_dependence_measures import
+    ... distance_covariance
+    >>> distance_covariance(np.random.random(1000), np.random.random(1000))
+    0.007575063951951362
+
+    """
+    return distance_statistics(x, y).distance_covariance
+
+
+def distance_variance(x):
+    r"""Distance variance.
+
+    Calculate the empirical distance variance as described in [1]_.
+
+    Parameters
+    ----------
+    x : array_like, 1-D or 2-D
+        If `x` is 1-D than it is assumed to be a vector of observations of a
+        single random variable. If `x` is 2-D than the rows should be
+        observations and the columns are treated as the components of a
+        random vector, i.e., each column represents a different component of
+        the random vector `x`.
+
+    Returns
+    -------
+    float
+        The empirical distance variance of `x`.
+
+    References
+    ----------
+    .. [1] Szekely, G.J., Rizzo, M.L., and Bakirov, N.K. (2007)
+       "Measuring and testing dependence by correlation of distances".
+       Annals of Statistics, Vol. 35 No. 6, pp. 2769-2794.
+
+    Examples
+    --------
+
+    >>> from statsmodels.stats.dist_dependence_measures import
+    ... distance_variance
+    >>> distance_variance(np.random.random(1000))
+    0.21732609190659702
+
+    """
+    return distance_covariance(x, x)
+
+
+def distance_correlation(x, y):
+    r"""Distance correlation.
+
+    Calculate the empirical distance correlation as described in [1]_.
+    This statistic is analogous to product-moment correlation and describes
+    the dependence between `x` and `y`, which are random vectors of
+    arbitrary length. The statistics' values range between 0 (implies
+    independence) and 1 (implies complete dependence).
+
+    Parameters
+    ----------
+    x : array_like, 1-D or 2-D
+        If `x` is 1-D than it is assumed to be a vector of observations of a
+        single random variable. If `x` is 2-D than the rows should be
+        observations and the columns are treated as the components of a
+        random vector, i.e., each column represents a different component of
+        the random vector `x`.
+    y : array_like, 1-D or 2-D
+        Same as `x`, but only the number of observation has to match that of
+        `x`. If `y` is 2-D note that the number of columns of `y` (i.e., the
+        number of components in the random vector) does not need to match
+        the number of columns in `x`.
+
+    Returns
+    -------
+    float
+        The empirical distance correlation between `x` and `y`.
+
+    References
+    ----------
+    .. [1] Szekely, G.J., Rizzo, M.L., and Bakirov, N.K. (2007)
+       "Measuring and testing dependence by correlation of distances".
+       Annals of Statistics, Vol. 35 No. 6, pp. 2769-2794.
+
+    Examples
+    --------
+
+    >>> from statsmodels.stats.dist_dependence_measures import
+    ... distance_correlation
+    >>> distance_correlation(np.random.random(1000), np.random.random(1000))
+    0.04060497840149489
+
+    """
+    return distance_statistics(x, y).distance_correlation


Test patch:
diff --git a/statsmodels/stats/tests/test_dist_dependant_measures.py b/statsmodels/stats/tests/test_dist_dependant_measures.py
new file mode 100644
index 000000000..1e5fbcb06
--- /dev/null
+++ b/statsmodels/stats/tests/test_dist_dependant_measures.py
@@ -0,0 +1,198 @@
+import numpy as np
+from numpy.testing import assert_almost_equal
+from pytest import raises as assert_raises, warns as assert_warns
+
+import statsmodels.stats.dist_dependence_measures as ddm
+from statsmodels.datasets import get_rdataset
+
+
+class TestDistDependenceMeasures(object):
+    @classmethod
+    def setup_class(cls):
+        """
+        Values were obtained via the R `energy` package.
+
+        R code:
+        ------
+        > dcov.test(x, y, R=200)
+
+            dCov independence test (permutation test)
+
+        data:  index 1, replicates 200
+        nV^2 = 45829, p-value = 0.004975
+        sample estimates:
+            dCov
+        47.86925
+
+        > DCOR(x, y)
+        $dCov
+        [1] 47.86925
+
+        $dCor
+        [1] 0.9999704
+
+        $dVarX
+        [1] 47.28702
+
+        $dVarY
+        [1] 48.46151
+        """
+        np.random.seed(3)
+
+        cls.x = np.array(range(1, 101)).reshape((20, 5))
+        cls.y = cls.x + np.log(cls.x)
+
+        cls.dcor_exp = 0.9999704
+        cls.dcov_exp = 47.86925
+        cls.dvar_x_exp = 47.28702
+        cls.dvar_y_exp = 48.46151
+        cls.pval_emp_exp = 0.004975
+        cls.test_stat_emp_exp = 45829
+
+        # The values above are functions of the following values, and
+        # therefore when the above group of variables is computed correctly
+        # it means this group of variables was also correctly calculated.
+        cls.S_exp = 5686.03162
+        cls.test_stat_asym_exp = 2.8390102
+        cls.pval_asym_exp = 0.00452
+
+    def test_input_validation_nobs(self):
+        with assert_raises(ValueError, match="same number of observations"):
+            ddm.distance_covariance_test(self.x[:2, :], self.y)
+
+    def test_input_validation_unknown_method(self):
+        with assert_raises(ValueError, match="Unknown 'method' parameter"):
+            ddm.distance_covariance_test(self.x, self.y, method="wrong_name")
+
+    def test_statistic_value_asym_method(self):
+        statistic, pval, method = ddm.distance_covariance_test(
+            self.x, self.y, method="asym")
+
+        assert method == "asym"
+        assert_almost_equal(statistic, self.test_stat_asym_exp, 4)
+        assert_almost_equal(pval, self.pval_asym_exp, 3)
+
+    def test_statistic_value_emp_method(self):
+        statistic, pval, method = ddm.distance_covariance_test(
+            self.x, self.y, method="emp")
+
+        assert method == "emp"
+        assert_almost_equal(statistic, self.test_stat_emp_exp, 0)
+        assert_almost_equal(pval, self.pval_emp_exp, 3)
+
+    def test_fallback_to_asym_method(self):
+        match_text = "The asymptotic approximation will be used"
+        with assert_warns(UserWarning, match=match_text):
+            statistic, pval, _ = ddm.distance_covariance_test(
+                self.x, self.y, method="emp", B=200
+            )
+            assert_almost_equal(statistic, self.test_stat_emp_exp, 0)
+            assert_almost_equal(pval, self.pval_asym_exp, 3)
+
+    def test_statistics_for_2d_input(self):
+        stats = ddm.distance_statistics(self.x, self.y)
+
+        assert_almost_equal(stats.test_statistic, self.test_stat_emp_exp, 0)
+        assert_almost_equal(stats.distance_correlation, self.dcor_exp, 4)
+        assert_almost_equal(stats.distance_covariance, self.dcov_exp, 4)
+        assert_almost_equal(stats.dvar_x, self.dvar_x_exp, 4)
+        assert_almost_equal(stats.dvar_y, self.dvar_y_exp, 4)
+        assert_almost_equal(stats.S, self.S_exp, 4)
+
+    def test_statistics_for_1d_input(self):
+        x = np.array(range(1, 21))
+        y = x + np.log(x)
+
+        stats = ddm.distance_statistics(x, y)
+
+        # Values were obtained using the R `energy` package
+        assert_almost_equal(stats.test_statistic, 398.94623, 5)
+        assert_almost_equal(stats.distance_correlation, 0.9996107, 4)
+        assert_almost_equal(stats.distance_covariance, 4.4662414, 4)
+        assert_almost_equal(stats.dvar_x, 4.2294799, 4)
+        assert_almost_equal(stats.dvar_y, 4.7199304, 4)
+        assert_almost_equal(stats.S, 49.880200, 4)
+
+    def test_results_on_the_iris_dataset(self):
+        """
+        R code example from the `energy` package documentation for
+        `energy::distance_covariance.test`:
+
+        > x <- iris[1:50, 1:4]
+        > y <- iris[51:100, 1:4]
+        > set.seed(1)
+        > dcov.test(x, y, R=200)
+
+            dCov independence test (permutation test)
+
+        data:  index 1, replicates 200
+        nV^2 = 0.5254, p-value = 0.9552
+        sample estimates:
+             dCov
+        0.1025087
+        """
+        iris = get_rdataset("iris").data.values[:, :4]
+        x = iris[:50]
+        y = iris[50:100]
+
+        stats = ddm.distance_statistics(x, y)
+
+        assert_almost_equal(stats.test_statistic, 0.5254, 4)
+        assert_almost_equal(stats.distance_correlation, 0.3060479, 4)
+        assert_almost_equal(stats.distance_covariance, 0.1025087, 4)
+        assert_almost_equal(stats.dvar_x, 0.2712927, 4)
+        assert_almost_equal(stats.dvar_y, 0.4135274, 4)
+        assert_almost_equal(stats.S, 0.667456, 4)
+
+        test_statistic, _, method = ddm.distance_covariance_test(x, y, B=199)
+
+        assert_almost_equal(test_statistic, 0.5254, 4)
+        assert method == "emp"
+
+    def test_results_on_the_quakes_dataset(self):
+        """
+        R code:
+        ------
+
+        > data("quakes")
+        > x = quakes[1:50, 1:3]
+        > y = quakes[51:100, 1:3]
+        > dcov.test(x, y, R=200)
+
+            dCov independence test (permutation test)
+
+        data:  index 1, replicates 200
+        nV^2 = 45046, p-value = 0.4577
+        sample estimates:
+            dCov
+        30.01526
+        """
+        quakes = get_rdataset("quakes").data.values[:, :3]
+        x = quakes[:50]
+        y = quakes[50:100]
+
+        stats = ddm.distance_statistics(x, y)
+
+        assert_almost_equal(np.round(stats.test_statistic), 45046, 0)
+        assert_almost_equal(stats.distance_correlation, 0.1894193, 4)
+        assert_almost_equal(stats.distance_covariance, 30.01526, 4)
+        assert_almost_equal(stats.dvar_x, 170.1702, 4)
+        assert_almost_equal(stats.dvar_y, 147.5545, 4)
+        assert_almost_equal(stats.S, 52265, 0)
+
+        test_statistic, _, method = ddm.distance_covariance_test(x, y, B=199)
+
+        assert_almost_equal(np.round(test_statistic), 45046, 0)
+        assert method == "emp"
+
+    def test_dcor(self):
+        assert_almost_equal(ddm.distance_correlation(self.x, self.y),
+                            self.dcor_exp, 4)
+
+    def test_dcov(self):
+        assert_almost_equal(ddm.distance_covariance(self.x, self.y),
+                            self.dcov_exp, 4)
+
+    def test_dvar(self):
+        assert_almost_equal(ddm.distance_variance(self.x),
+                            self.dvar_x_exp, 4)