Merge pull request #307 from raybellwaves/add-linslope

Author: raybellwaves

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
 ~~~~~~~~~~~~~~~~

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 ..

ci/minimum-tests.yml

@@ -23,4 +23,4 @@ dependencies:
   - pytest-xdist
   - pip
   - pip:
-    - -e ..
+      - -e ..

docs/source/api.rst

@@ -24,6 +24,7 @@ Correlation Metrics
     spearman_r_eff_p_value
     effective_sample_size
     r2
+    linslope
 
 Distance Metrics
 ~~~~~~~~~~~~~~~~

docs/source/api/xskillscore.linslope.rst

@@ -0,0 +1,6 @@
+xskillscore.linslope
+====================
+
+.. currentmodule:: xskillscore
+
+.. autofunction:: linslope

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",

xskillscore/__init__.py

@@ -7,6 +7,7 @@
 from .core.contingency import Contingency
 from .core.deterministic import (
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,

xskillscore/core/accessor.py

@@ -2,6 +2,7 @@
 
 from .deterministic import (
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
@@ -46,6 +47,11 @@ def _in_ds(self, x):
         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)

xskillscore/core/deterministic.py

@@ -4,6 +4,7 @@
 
 from .np_deterministic import (
     _effective_sample_size,
+    _linslope,
     _mae,
     _mape,
     _me,
@@ -27,21 +28,22 @@
 )
 
 __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.
 

xskillscore/core/np_deterministic.py

@@ -6,21 +6,22 @@
 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.
 

xskillscore/tests/test_accessor_deterministic.py

@@ -4,6 +4,7 @@
 
 from xskillscore.core.deterministic import (
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
@@ -21,6 +22,7 @@
 )
 
 correlation_metrics = [
+    linslope,
     pearson_r,
     r2,
     pearson_r_p_value,

xskillscore/tests/test_deterministic.py

@@ -7,6 +7,7 @@
     _preprocess_dims,
     _preprocess_weights,
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
@@ -24,6 +25,7 @@
 )
 from xskillscore.core.np_deterministic import (
     _effective_sample_size,
+    _linslope,
     _mae,
     _mape,
     _me,
@@ -41,6 +43,7 @@
 )
 
 correlation_metrics = [
+    (linslope, _linslope),
     (pearson_r, _pearson_r),
     (r2, _r2),
     (pearson_r_p_value, _pearson_r_p_value),
@@ -51,6 +54,7 @@
     (effective_sample_size, _effective_sample_size),
 ]
 correlation_metrics_names = [
+    "linslope",
     "pearson_r",
     "r2",
     "pearson_r_p_value",

xskillscore/tests/test_gridded_metrics.py

@@ -2,6 +2,7 @@
 import pytest
 
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -17,6 +18,7 @@
 )
 
 METRICS = [
+    linslope,
     mae,
     mse,
     median_absolute_error,

xskillscore/tests/test_mask_skipna.py

@@ -2,6 +2,7 @@
 import pytest
 
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     median_absolute_error,

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 @@
 
 import xskillscore as xs
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -36,6 +37,7 @@
 ]
 
 xs_scipy_metrics = [
+    (linslope, linregress, 0),
     (pearson_r, pearsonr, 0),
     (spearman_r, spearmanr, 0),
     (pearson_r_p_value, pearsonr, 1),

xskillscore/tests/test_skipna_functionality.py

@@ -4,6 +4,7 @@
 from xarray.tests import CountingScheduler, assert_allclose, raise_if_dask_computes
 
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -19,6 +20,7 @@
 )
 
 WEIGHTED_METRICS = [
+    linslope,
     pearson_r,
     pearson_r_p_value,
     spearman_r,