Repository: xarray-contrib/xskillscore

Issue:
metrics to add
Keep this comment in this issue updated with metrics to add to xskillscore. Consolidate from other issues as well as comments that appear below. The format for inputting will be:

 - METRIC_API_NAME (LONG_METRIC_NAME) [RELATED ISSUE (IF EXISTS)] {METRIC SOURCE/EQUATION}

The full list of metrics current in xskillscore can be found [here](https://github.com/xarray-contrib/xskillscore/blob/master/docs/source/api.rst). Remove issues from here once they are added.

Correlation Metrics
-------------------
 - `pearson_r_auto` (Pearson R Autocorrelation) [#205] {https://github.com/bradyrx/esmtools/blob/master/esmtools/stats.py#L171 }

Distance Metrics
-----------------
 - `medape` (Median Absolute Percentage Error)
 - `rmspe` (Root Mean Square Percentage Error) [#46] {https://www.kaggle.com/c/rossmann-store-sales/overview/evaluation }
 - `msle` (Mean Squared Log Error) [#47] {https://scikit-learn.org/stable/modules/generated/sklearn.metrics.mean_squared_log_error.html }
 - `crmse` (Centered Root Mean Square Error) {https://solarforecastarbiter.org/metrics/#crmse}
 - `maape` (Mean Arctangent Absolute Percentage Error) [#85] {https://github.com/xarray-contrib/xskillscore/issues/85#issuecomment-653135691 }
 - `explained_var` (Explained Variance Score) [#213] {https://scikit-learn.org/stable/modules/model_evaluation.html#explained-variance-score }
 - `mean_pinball` (Mean Pinball Loss) [#274] {https://scikit-learn.org/dev/modules/generated/sklearn.metrics.mean_pinball_loss.html}

Probabilistic Metrics
--------------------
- `brier_skill_score` (Brier Skill Score) [#49] {https://github.com/xarray-contrib/xskillscore/issues/49#issue-526878892 }
- Add `fair` arg to `crps_ensemble` [#260] {? }
- `crpss` (Continuous Ranked Probability Skill Score) [#49] {https://github.com/pangeo-data/climpred/blob/main/climpred/metrics.py#L2176}
- `rpss` (Ranked Probability Skill Score) [#49] {https://github.com/pangeo-data/climpred/blob/main/climpred/tests/test_probabilistic.py#L252}

Dichotomous-Only (yes/no) Metrics
------------------------------------
- `f1` (f1 Score) [#138] {https://scikit-learn.org/stable/modules/generated/sklearn.metrics.f1_score.html#sklearn.metrics.f1_score }
- `tpr` (True Positive Rate/Recall Score) [#138] {https://scikit-learn.org/stable/modules/generated/sklearn.metrics.recall_score.html#sklearn.metrics.recall_score }
- `precision` (Positive Predictive Value/Precision Score) {https://scikit-learn.org/stable/modules/generated/sklearn.metrics.precision_score.html#sklearn.metrics.precision_score }
- `rel_value` (Relative Value Score) [#229] {https://www.ecmwf.int/sites/default/files/elibrary/2007/15489-verification-probability-forecasts.pdf }

Multi-Category Metrics
-----------------------
- `mc_threat_score` (Multi-Category Threat Score) [#187] {?}

Comparative
-------------
- `ttest_ind` (T-test for the means of two independent samples of scores) [#175] {https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.ttest_ind.html }

Resampling
------------


Metric glossaries:
 - https://solarforecastarbiter.org/metrics/

Discussion:

Patch:
diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 9b2d148..e4e726a 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -12,6 +12,7 @@ Features
   observations are probability distributions ``p`` or cumulative
   distributionss ``c``. See :py:func:`~xskillscore.rps` docstrings and doctests for
   examples. (:pr:`300`) `Aaron Spring`_
+- Added slope of linear fit :py:func:`~xskillscore.linslope`. `Ray Bell`_
 
 Internal Changes
 ~~~~~~~~~~~~~~~~
diff --git a/ci/docs_notebooks.yml b/ci/docs_notebooks.yml
index 540ad2b..6e4b78d 100644
--- a/ci/docs_notebooks.yml
+++ b/ci/docs_notebooks.yml
@@ -25,6 +25,6 @@ dependencies:
   - sphinx_rtd_theme
   - pip
   - pip:
-    - sphinx_autosummary_accessors
-    # Install latest version of xskillscore.
-    - -e ..
+      - sphinx_autosummary_accessors
+      # Install latest version of xskillscore.
+      - -e ..
diff --git a/docs/source/api.rst b/docs/source/api.rst
index a59cbc4..a351298 100644
--- a/docs/source/api.rst
+++ b/docs/source/api.rst
@@ -24,6 +24,7 @@ Correlation Metrics
     spearman_r_eff_p_value
     effective_sample_size
     r2
+    linslope
 
 Distance Metrics
 ~~~~~~~~~~~~~~~~
diff --git a/docs/source/api/xskillscore.linslope.rst b/docs/source/api/xskillscore.linslope.rst
new file mode 100644
index 0000000..f714a49
--- /dev/null
+++ b/docs/source/api/xskillscore.linslope.rst
@@ -0,0 +1,6 @@
+xskillscore.linslope
+====================
+
+.. currentmodule:: xskillscore
+
+.. autofunction:: linslope
diff --git a/docs/source/quick-start.ipynb b/docs/source/quick-start.ipynb
index d31375e..d0e10e6 100644
--- a/docs/source/quick-start.ipynb
+++ b/docs/source/quick-start.ipynb
@@ -66,6 +66,7 @@
     "* Spearman Correlation p value (`spearman_r_p_value`)\n",
     "* Spearman Correlation effective p value (`spearman_r_eff_p_value`)\n",
     "* Effective Sample Size (`effective_sample_size`)\n",
+    "* Slope of Linear Fit (`linslope`)\n",
     "* Coefficient of Determination (`r2`)\n",
     "\n",
     "### Distance-Based\n",
diff --git a/xskillscore/__init__.py b/xskillscore/__init__.py
index ff9018c..9366a13 100644
--- a/xskillscore/__init__.py
+++ b/xskillscore/__init__.py
@@ -7,6 +7,7 @@ from .core.comparative import mae_test, sign_test
 from .core.contingency import Contingency
 from .core.deterministic import (
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
diff --git a/xskillscore/core/accessor.py b/xskillscore/core/accessor.py
index 9bfc41b..31a70f5 100644
--- a/xskillscore/core/accessor.py
+++ b/xskillscore/core/accessor.py
@@ -2,6 +2,7 @@ import xarray as xr
 
 from .deterministic import (
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
@@ -46,6 +47,11 @@ class XSkillScoreAccessor(object):
         else:
             return self._obj[x]
 
+    def linslope(self, a, b, *args, **kwargs):
+        a = self._in_ds(a)
+        b = self._in_ds(b)
+        return linslope(a, b, *args, **kwargs)
+
     def pearson_r(self, a, b, *args, **kwargs):
         a = self._in_ds(a)
         b = self._in_ds(b)
diff --git a/xskillscore/core/deterministic.py b/xskillscore/core/deterministic.py
index 201f866..9bae56d 100644
--- a/xskillscore/core/deterministic.py
+++ b/xskillscore/core/deterministic.py
@@ -4,6 +4,7 @@ import xarray as xr
 
 from .np_deterministic import (
     _effective_sample_size,
+    _linslope,
     _mae,
     _mape,
     _me,
@@ -27,21 +28,22 @@ from .utils import (
 )
 
 __all__ = [
+    "effective_sample_size",
+    "linslope",
+    "mae",
+    "mape",
+    "me",
+    "median_absolute_error",
+    "mse",
     "pearson_r",
-    "pearson_r_p_value",
     "pearson_r_eff_p_value",
-    "me",
+    "pearson_r_p_value",
+    "r2",
     "rmse",
-    "mse",
-    "mae",
-    "median_absolute_error",
     "smape",
-    "mape",
     "spearman_r",
-    "spearman_r_p_value",
     "spearman_r_eff_p_value",
-    "effective_sample_size",
-    "r2",
+    "spearman_r_p_value",
 ]
 
 
@@ -71,6 +73,75 @@ def _determine_input_core_dims(dim, weights):
     return input_core_dims
 
 
+def linslope(a, b, dim=None, weights=None, skipna=False, keep_attrs=False):
+    """Slope of linear fit.
+
+    .. math::
+        s_{ab} = \\frac{ \\sum_{i=i}^{n} (a_{i} - \\bar{a}) (b_{i} - \\bar{b}) }
+                 { \\sum_{i=1}^{n} (a_{i} - \\bar{a})^{2} }
+
+    Parameters
+    ----------
+    a : xarray.Dataset or xarray.DataArray
+        Labeled array(s) over which to apply the function.
+    b : xarray.Dataset or xarray.DataArray
+        Labeled array(s) over which to apply the function.
+    dim : str, list
+        The dimension(s) to apply the correlation along. Note that this dimension will
+        be reduced as a result. Defaults to None reducing all dimensions.
+    weights : xarray.Dataset or xarray.DataArray or None
+        Weights matching dimensions of ``dim`` to apply during the function.
+    skipna : bool
+        If True, skip NaNs when computing function.
+    keep_attrs : bool
+        If True, the attributes (attrs) will be copied
+        from the first input to the new one.
+        If False (default), the new object will
+        be returned without attributes.
+
+    Returns
+    -------
+    xarray.DataArray or xarray.Dataset
+        Slope of linear fit.
+
+    See Also
+    --------
+    scipy.stats.linregress
+
+    Examples
+    --------
+    >>> a = xr.DataArray(np.random.rand(5, 3, 3),
+    ...                  dims=['time', 'x', 'y'])
+    >>> b = xr.DataArray(np.random.rand(5, 3, 3),
+    ...                  dims=['time', 'x', 'y'])
+    >>> xs.linslope(a, b, dim='time')
+    <xarray.DataArray (x: 3, y: 3)>
+    array([[-0.30948771, -0.21562529, -0.63141304],
+           [ 0.31446077,  2.23858011,  0.44743617],
+           [-0.22243944,  0.47034784,  1.08512859]])
+    Dimensions without coordinates: x, y
+    """
+    _fail_if_dim_empty(dim)
+    dim, _ = _preprocess_dims(dim, a)
+    a, b = xr.broadcast(a, b, exclude=dim)
+    a, b, new_dim, weights = _stack_input_if_needed(a, b, dim, weights)
+    weights = _preprocess_weights(a, dim, new_dim, weights)
+
+    input_core_dims = _determine_input_core_dims(new_dim, weights)
+
+    return xr.apply_ufunc(
+        _linslope,
+        a,
+        b,
+        weights,
+        input_core_dims=input_core_dims,
+        kwargs={"axis": -1, "skipna": skipna},
+        dask="parallelized",
+        output_dtypes=[float],
+        keep_attrs=keep_attrs,
+    )
+
+
 def pearson_r(a, b, dim=None, weights=None, skipna=False, keep_attrs=False):
     """Pearson's correlation coefficient.
 
diff --git a/xskillscore/core/np_deterministic.py b/xskillscore/core/np_deterministic.py
index 180cbd3..aada365 100644
--- a/xskillscore/core/np_deterministic.py
+++ b/xskillscore/core/np_deterministic.py
@@ -6,21 +6,22 @@ from scipy.stats import distributions
 from .utils import suppress_warnings
 
 __all__ = [
+    "_effective_sample_size",
+    "_linslope",
+    "_mae",
+    "_mape",
+    "_me",
+    "_median_absolute_error",
+    "_mse",
     "_pearson_r",
-    "_pearson_r_p_value",
     "_pearson_r_eff_p_value",
-    "_me",
+    "_pearson_r_p_value",
+    "_r2",
     "_rmse",
-    "_mse",
-    "_mae",
-    "_median_absolute_error",
     "_smape",
-    "_mape",
     "_spearman_r",
-    "_spearman_r_p_value",
     "_spearman_r_eff_p_value",
-    "_effective_sample_size",
-    "_r2",
+    "_spearman_r_p_value",
 ]
 
 
@@ -153,6 +154,55 @@ def _effective_sample_size(a, b, axis, skipna):
     return n_eff
 
 
+def _linslope(a, b, weights, axis, skipna):
+    """ndarray implementation of scipy.stats.linregress[slope].
+
+    Parameters
+    ----------
+    a : ndarray
+        Input array.
+    b : ndarray
+        Input array.
+    axis : int
+        The axis to apply the linear slope along.
+    weights : ndarray
+        Input array of weights for a and b.
+    skipna : bool
+        If True, skip NaNs when computing function.
+
+    Returns
+    -------
+    res : ndarray
+        slope of linear fit.
+
+    See Also
+    --------
+    scipy.stats.linregress
+    """
+    sumfunc, meanfunc = _get_numpy_funcs(skipna)
+    if skipna:
+        a, b, weights = _match_nans(a, b, weights)
+    weights = _check_weights(weights)
+    a = np.rollaxis(a, axis)
+    b = np.rollaxis(b, axis)
+    if weights is not None:
+        weights = np.rollaxis(weights, axis)
+
+    am, bm = __compute_anomalies(a, b, weights=weights, axis=0, skipna=skipna)
+
+    if weights is not None:
+        s_num = sumfunc(weights * am * bm, axis=0)
+        s_den = sumfunc(weights * am * am, axis=0)
+    else:
+        s_num = sumfunc(am * bm, axis=0)
+        s_den = sumfunc(am * am, axis=0)
+
+    with suppress_warnings("invalid value encountered in true_divide"):
+        with suppress_warnings("invalid value encountered in double_scalars"):
+            res = s_num / s_den
+    return res
+
+
 def _r2(a, b, weights, axis, skipna):
     """ndarray implementation of sklearn.metrics.r2_score.
 


Test patch:
diff --git a/ci/minimum-tests.yml b/ci/minimum-tests.yml
index ce70b61..69f4717 100644
--- a/ci/minimum-tests.yml
+++ b/ci/minimum-tests.yml
@@ -23,4 +23,4 @@ dependencies:
   - pytest-xdist
   - pip
   - pip:
-    - -e ..
+      - -e ..
diff --git a/xskillscore/tests/test_accessor_deterministic.py b/xskillscore/tests/test_accessor_deterministic.py
index 1dd0b7b..d36271b 100644
--- a/xskillscore/tests/test_accessor_deterministic.py
+++ b/xskillscore/tests/test_accessor_deterministic.py
@@ -4,6 +4,7 @@ from xarray.tests import assert_allclose
 
 from xskillscore.core.deterministic import (
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
@@ -21,6 +22,7 @@ from xskillscore.core.deterministic import (
 )
 
 correlation_metrics = [
+    linslope,
     pearson_r,
     r2,
     pearson_r_p_value,
diff --git a/xskillscore/tests/test_deterministic.py b/xskillscore/tests/test_deterministic.py
index 0becc45..bf97a1a 100644
--- a/xskillscore/tests/test_deterministic.py
+++ b/xskillscore/tests/test_deterministic.py
@@ -7,6 +7,7 @@ from xskillscore.core.deterministic import (
     _preprocess_dims,
     _preprocess_weights,
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
@@ -24,6 +25,7 @@ from xskillscore.core.deterministic import (
 )
 from xskillscore.core.np_deterministic import (
     _effective_sample_size,
+    _linslope,
     _mae,
     _mape,
     _me,
@@ -41,6 +43,7 @@ from xskillscore.core.np_deterministic import (
 )
 
 correlation_metrics = [
+    (linslope, _linslope),
     (pearson_r, _pearson_r),
     (r2, _r2),
     (pearson_r_p_value, _pearson_r_p_value),
@@ -51,6 +54,7 @@ correlation_metrics = [
     (effective_sample_size, _effective_sample_size),
 ]
 correlation_metrics_names = [
+    "linslope",
     "pearson_r",
     "r2",
     "pearson_r_p_value",
diff --git a/xskillscore/tests/test_gridded_metrics.py b/xskillscore/tests/test_gridded_metrics.py
index 180abe1..cf9009f 100644
--- a/xskillscore/tests/test_gridded_metrics.py
+++ b/xskillscore/tests/test_gridded_metrics.py
@@ -2,6 +2,7 @@ import numpy as np
 import pytest
 
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -17,6 +18,7 @@ from xskillscore.core.deterministic import (
 )
 
 METRICS = [
+    linslope,
     mae,
     mse,
     median_absolute_error,
diff --git a/xskillscore/tests/test_mask_skipna.py b/xskillscore/tests/test_mask_skipna.py
index 879e3b9..7f038e4 100644
--- a/xskillscore/tests/test_mask_skipna.py
+++ b/xskillscore/tests/test_mask_skipna.py
@@ -2,6 +2,7 @@ import numpy as np
 import pytest
 
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     median_absolute_error,
diff --git a/xskillscore/tests/test_metric_results_accurate.py b/xskillscore/tests/test_metric_results_accurate.py
index 8899540..8b40209 100644
--- a/xskillscore/tests/test_metric_results_accurate.py
+++ b/xskillscore/tests/test_metric_results_accurate.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 import sklearn.metrics
-from scipy.stats import pearsonr, spearmanr
+from scipy.stats import linregress, pearsonr, spearmanr
 from sklearn.metrics import (
     mean_absolute_error,
     mean_absolute_percentage_error,
@@ -11,6 +11,7 @@ from sklearn.metrics import (
 
 import xskillscore as xs
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -36,6 +37,7 @@ xs_skl_metrics_with_zeros = [
 ]
 
 xs_scipy_metrics = [
+    (linslope, linregress, 0),
     (pearson_r, pearsonr, 0),
     (spearman_r, spearmanr, 0),
     (pearson_r_p_value, pearsonr, 1),
diff --git a/xskillscore/tests/test_skipna_functionality.py b/xskillscore/tests/test_skipna_functionality.py
index 50adb72..7024cf9 100644
--- a/xskillscore/tests/test_skipna_functionality.py
+++ b/xskillscore/tests/test_skipna_functionality.py
@@ -4,6 +4,7 @@ import pytest
 from xarray.tests import CountingScheduler, assert_allclose, raise_if_dask_computes
 
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -19,6 +20,7 @@ from xskillscore.core.deterministic import (
 )
 
 WEIGHTED_METRICS = [
+    linslope,
     pearson_r,
     pearson_r_p_value,
     spearman_r,
diff --git a/xskillscore/tests/test_weighted_metric_results_accurate.py b/xskillscore/tests/test_weighted_metric_results_accurate.py
index 5ea8e63..6de2338 100644
--- a/xskillscore/tests/test_weighted_metric_results_accurate.py
+++ b/xskillscore/tests/test_weighted_metric_results_accurate.py
@@ -1,7 +1,5 @@
 import numpy as np
 import pytest
-import sklearn.metrics
-from scipy.stats import pearsonr
 from sklearn.metrics import (
     mean_absolute_error,
     mean_absolute_percentage_error,
@@ -11,6 +9,7 @@ from sklearn.metrics import (
 
 import xskillscore as xs
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -44,7 +43,15 @@ def weighted_pearsonr(x, y, w):
     return r_num / r_den
 
 
-xs_scipy_metrics = [(pearson_r, weighted_pearsonr)]
+def weighted_linslope(x, y, w):
+    xm = x - (np.sum(x * w) / np.sum(w))
+    ym = y - (np.sum(y * w) / np.sum(w))
+    s_num = np.sum(w * xm * ym)
+    s_den = np.sum(w * xm * xm)
+    return s_num / s_den
+
+
+xs_scipy_metrics = [(pearson_r, weighted_pearsonr), (linslope, weighted_linslope)]
 
 
 xs_np_metrics = [