Repository: fatiando/verde
Pull request #162 (2018-12-19 22:05:50)

Allow zero degree polynomials in Trend (a mean)
There was no real reason to not allow them in the past. A zero degree
polynomial is a mean. This can be useful to test if a simple static
shift in the data is enough of a trend.


**Reminders**

- [ ] Run `make format` and `make check` to make sure the code follows the style guide.
- [ ] Add tests for new features or tests that would have caught the bug that you're fixing.
- [ ] Add new public functions/methods/classes to `doc/api/index.rst` and `verde/__init__.py`.
- [ ] Write detailed docstrings for all functions/classes/methods.
- [ ] If adding new functionality, add an example to the docstring, gallery, and/or tutorials.



Patch:
diff --git a/verde/trend.py b/verde/trend.py
index f51d8c2..c574fe5 100644
--- a/verde/trend.py
+++ b/verde/trend.py
@@ -9,20 +9,28 @@ from .coordinates import get_region
 
 
 class Trend(BaseGridder):
-    """
+    r"""
     Fit a 2D polynomial trend to spatial data.
 
-    The trend is estimated through weighted least-squares regression.
+    The polynomial of degree :math:`N` is defined as:
+
+    .. math::
+
+        f(e, n) = \sum\limits_{l=0}^{N}\sum\limits_{m=0}^{N - l} e^l n^m
 
-    The Jacobian (design, sensitivity, feature, etc) matrix for the regression
-    is normalized using :class:`sklearn.preprocessing.StandardScaler` without
-    centering the mean so that the transformation can be undone in the
-    estimated coefficients.
+    in which :math:`e` and :math:`n` are the easting and northing coordinates,
+    respectively.
+
+    The trend is estimated through weighted least-squares regression. The Jacobian
+    (design, sensitivity, feature, etc) matrix for the regression is normalized using
+    :class:`sklearn.preprocessing.StandardScaler` without centering the mean so that the
+    transformation can be undone in the estimated coefficients.
 
     Parameters
     ----------
     degree : int
-        The degree of the polynomial. Must be >= 1.
+        The degree of the polynomial. Must be >= 0 (a degree of zero would estimate the
+        mean of the data).
 
     Attributes
     ----------
@@ -37,22 +45,37 @@ class Trend(BaseGridder):
     --------
 
     >>> from verde import grid_coordinates
+    >>> import numpy as np
     >>> coordinates = grid_coordinates((1, 5, -5, -1), shape=(5, 5))
     >>> data = 10 + 2*coordinates[0] - 0.4*coordinates[1]
     >>> trend = Trend(degree=1).fit(coordinates, data)
-    >>> print(', '.join(['{:.1f}'.format(i) for i in trend.coef_]))
-    10.0, 2.0, -0.4
-    >>> import numpy as np
+    >>> print("Coefficients:", ', '.join(['{:.1f}'.format(i) for i in trend.coef_]))
+    Coefficients: 10.0, 2.0, -0.4
     >>> np.allclose(trend.predict(coordinates), data)
     True
-    >>> # Use weights to account for outliers
+
+    A zero degree polynomial estimates the mean of the data:
+
+    >>> mean = Trend(degree=0).fit(coordinates, data)
+    >>> np.allclose(mean.predict(coordinates), data.mean())
+    True
+    >>> print("Data mean:", '{:.2f}'.format(data.mean()))
+    Data mean: 17.20
+    >>> print("Coefficient:", '{:.2f}'.format(mean.coef_[0]))
+    Coefficient: 17.20
+
+    We can use weights to account for outliers or data points with variable
+    uncertainties (see :func:`verde.variance_to_weights`):
+
     >>> data_out = data.copy()
     >>> data_out[2, 2] += 500
     >>> weights = np.ones_like(data)
     >>> weights[2, 2] = 1e-10
     >>> trend_out = Trend(degree=1).fit(coordinates, data_out, weights)
-    >>> print(', '.join(['{:.1f}'.format(i) for i in trend_out.coef_]))
-    10.0, 2.0, -0.4
+    >>> # Still recover the coefficients even with the added outlier
+    >>> print("Coefficients:", ', '.join(['{:.1f}'.format(i) for i in trend_out.coef_]))
+    Coefficients: 10.0, 2.0, -0.4
+    >>> # The residual at the outlier location should be values we added to that point
     >>> residual = data_out - trend_out.predict(coordinates)
     >>> print('{:.2f}'.format(residual[2, 2]))
     500.00
@@ -205,9 +228,12 @@ def polynomial_power_combinations(degree):
     >>> # This is a long polynomial so split it in two lines
     >>> print(" ".join([str(c) for c in polynomial_power_combinations(3)]))
     (0, 0) (1, 0) (0, 1) (2, 0) (1, 1) (0, 2) (3, 0) (2, 1) (1, 2) (0, 3)
+    >>> # A degree zero polynomial would be just the mean
+    >>> print(polynomial_power_combinations(0))
+    ((0, 0),)
 
     """
-    if degree < 1:
-        raise ValueError("Invalid polynomial degree '{}'. Must be >= 1.".format(degree))
+    if degree < 0:
+        raise ValueError("Invalid polynomial degree '{}'. Must be >= 0.".format(degree))
     combinations = ((i, j) for j in range(degree + 1) for i in range(degree + 1 - j))
     return tuple(sorted(combinations, key=sum))


Test patch:
diff --git a/verde/tests/test_trend.py b/verde/tests/test_trend.py
index 3a21a3b..f8264e8 100644
--- a/verde/tests/test_trend.py
+++ b/verde/tests/test_trend.py
@@ -23,10 +23,20 @@ def test_trend(simple_model):
     "Test the trend estimation on a simple problem"
     coords, coefs, data = simple_model
     trend = Trend(degree=1).fit(coords, data)
+    assert trend.coef_.size == 3
     npt.assert_allclose(trend.coef_, coefs)
     npt.assert_allclose(trend.predict(coords), data)
 
 
+def test_trend_mean(simple_model):
+    "Test the trend estimation with degree zero (a mean) on a simple problem"
+    coords, _, data = simple_model
+    trend = Trend(degree=0).fit(coords, data)
+    assert trend.coef_.size == 1
+    npt.assert_allclose(trend.coef_, data.mean())
+    npt.assert_allclose(trend.predict(coords), data.mean())
+
+
 def test_trend_weights(simple_model):
     "Use weights to account for outliers"
     coords, coefs, data = simple_model
@@ -43,14 +53,13 @@ def test_trend_weights(simple_model):
 
 def test_polynomial_combinations_fails():
     "Test failing conditions for the combinations"
-    with pytest.raises(ValueError):
-        polynomial_power_combinations(degree=0)
     with pytest.raises(ValueError):
         polynomial_power_combinations(degree=-10)
 
 
 def test_trend_jacobian_fails():
     "Test failing conditions for the trend jacobian builder"
+    # Coordinates must have the same number of elements
     east, north = np.arange(50), np.arange(30)
     trend = Trend(degree=1)
     with pytest.raises(ValueError):