Merge pull request #6360 from jonathanunderwood/legfit_restrict_terms

ENH: Allow specification of terms to fit in polynomial fitting functions

15 files changed, 402 insertions, 35 deletions

diff --git a/doc/release/1.11.0-notes.rst b/doc/release/1.11.0-notes.rst
index e541a4739..ea0e41694 100644
--- a/doc/release/1.11.0-notes.rst
+++ b/doc/release/1.11.0-notes.rst
@@ -240,6 +240,15 @@ arguments cannot be cast to a common type, it could have raised a ``TypeError``
 or ``ValueError`` depending on their order. Now, ``np.dot`` will now always
 raise a ``TypeError``.
 
+numpy.polynomial.*fit now supports restricted fitting
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+The ``deg`` parameter was extended to allow restricted fitting of
+specified terms in the polynomial expansion for all polynomial
+types. This change is backward compatible and it is still possible to
+specify ``deg`` as a single integer to specify the maximum
+order/degree of polynomial used in the fit, but it is now possible for
+``deg`` to be a list specifying which terms in the series to fit.
+
 Deprecations
 ============
 
diff --git a/numpy/polynomial/_polybase.py b/numpy/polynomial/_polybase.py
index 41bdf0484..37eb59f03 100644
--- a/numpy/polynomial/_polybase.py
+++ b/numpy/polynomial/_polybase.py
@@ -739,8 +739,14 @@ class ABCPolyBase(object):
             y-coordinates of the sample points. Several data sets of sample
             points sharing the same x-coordinates can be fitted at once by
             passing in a 2D-array that contains one dataset per column.
-        deg : int
-            Degree of the fitting polynomial.
+        deg : int or array_like
+            Degree of the fitting polynomial. If `deg` is a single integer
+            all terms up to and including the `deg`'th term are included.
+            `deg` may alternatively be a list or array specifying which
+            terms in the Legendre expansion to include in the fit.
+
+            .. versionchanged:: 1.11.0
+            `deg` may be a list specifying which terms to fit
         domain : {None, [beg, end], []}, optional
             Domain to use for the returned series. If ``None``,
             then a minimal domain that covers the points `x` is chosen.  If
diff --git a/numpy/polynomial/chebyshev.py b/numpy/polynomial/chebyshev.py
index de5cbb734..9db613b78 100644
--- a/numpy/polynomial/chebyshev.py
+++ b/numpy/polynomial/chebyshev.py
@@ -1617,8 +1617,14 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
         y-coordinates of the sample points. Several data sets of sample
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
-    deg : int
-        Degree of the fitting series
+    deg : int or array_like
+        Degree of the fitting series. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included.
+        `deg` may alternatively be a list or array specifying which
+        terms in the Legendre expansion to include in the fit.
+
+        .. versionchanged:: 1.11.0
+        `deg` may be a list specifying which terms to fit
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1710,12 +1716,14 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
     --------
 
     """
-    order = int(deg) + 1
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
+    deg = np.asarray([deg,], dtype=int).flatten()
 
     # check arguments.
-    if deg < 0:
+    if deg.size < 1:
+        raise TypeError("expected deg to be one or more integers")
+    if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
         raise TypeError("expected 1D vector for x")
@@ -1726,8 +1734,20 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
+    if deg.size == 1:
+        restricted_fit = False
+        lmax = deg[0]
+        order = lmax + 1
+    else:
+        restricted_fit = True
+        lmax = deg.max()
+        order = deg.size
+
     # set up the least squares matrices in transposed form
-    lhs = chebvander(x, deg).T
+    van = chebvander(x, lmax)
+    if restricted_fit:
+        van = van[:, deg]
+    lhs = van.T
     rhs = y.T
     if w is not None:
         w = np.asarray(w) + 0.0
@@ -1755,6 +1775,15 @@ def chebfit(x, y, deg, rcond=None, full=False, w=None):
     c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
     c = (c.T/scl).T
 
+    # Expand c to include non-fitted coefficients which are set to zero
+    if restricted_fit:
+        if c.ndim == 2:
+            cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
+        else:
+            cc = np.zeros(lmax+1, dtype=c.dtype)
+        cc[deg] = c
+        c = cc
+
     # warn on rank reduction
     if rank != order and not full:
         msg = "The fit may be poorly conditioned"
diff --git a/numpy/polynomial/hermite.py b/numpy/polynomial/hermite.py
index 2ce1d97a8..5d4b357fe 100644
--- a/numpy/polynomial/hermite.py
+++ b/numpy/polynomial/hermite.py
@@ -1388,8 +1388,14 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
         y-coordinates of the sample points. Several data sets of sample
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
-    deg : int
-        Degree of the fitting polynomial
+    deg : int or array_like
+        Degree of the fitting polynomial. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included.
+        `deg` may alternatively be a list or array specifying which
+        terms in the Legendre expansion to include in the fit.
+
+        .. versionchanged:: 1.11.0
+        `deg` may be a list specifying which terms to fit
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1486,12 +1492,14 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
     array([ 0.97902637,  1.99849131,  3.00006   ])
 
     """
-    order = int(deg) + 1
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
+    deg = np.asarray([deg,], dtype=int).flatten()
 
     # check arguments.
-    if deg < 0:
+    if deg.size < 1:
+        raise TypeError("expected deg to be one or more integers")
+    if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
         raise TypeError("expected 1D vector for x")
@@ -1502,8 +1510,20 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
+    if deg.size == 1:
+        restricted_fit = False
+        lmax = deg[0]
+        order = lmax + 1
+    else:
+        restricted_fit = True
+        lmax = deg.max()
+        order = deg.size
+
     # set up the least squares matrices in transposed form
-    lhs = hermvander(x, deg).T
+    van = hermvander(x, lmax)
+    if restricted_fit:
+        van = van[:, deg]
+    lhs = van.T
     rhs = y.T
     if w is not None:
         w = np.asarray(w) + 0.0
@@ -1531,6 +1551,15 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
     c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
     c = (c.T/scl).T
 
+    # Expand c to include non-fitted coefficients which are set to zero
+    if restricted_fit:
+        if c.ndim == 2:
+            cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
+        else:
+            cc = np.zeros(lmax+1, dtype=c.dtype)
+        cc[deg] = c
+        c = cc
+
     # warn on rank reduction
     if rank != order and not full:
         msg = "The fit may be poorly conditioned"
diff --git a/numpy/polynomial/hermite_e.py b/numpy/polynomial/hermite_e.py
index 394f256a5..da441af83 100644
--- a/numpy/polynomial/hermite_e.py
+++ b/numpy/polynomial/hermite_e.py
@@ -1385,8 +1385,14 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
         y-coordinates of the sample points. Several data sets of sample
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
-    deg : int
-        Degree of the fitting polynomial
+    deg : int or array_like
+        Degree of the fitting polynomial. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included.
+        `deg` may alternatively be a list or array specifying which
+        terms in the Legendre expansion to include in the fit.
+
+        .. versionchanged:: 1.11.0
+        `deg` may be a list specifying which terms to fit
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1483,12 +1489,14 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
     array([ 1.01690445,  1.99951418,  2.99948696])
 
     """
-    order = int(deg) + 1
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
+    deg = np.asarray([deg,], dtype=int).flatten()
 
     # check arguments.
-    if deg < 0:
+    if deg.size < 1:
+        raise TypeError("expected deg to be one or more integers")
+    if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
         raise TypeError("expected 1D vector for x")
@@ -1499,8 +1507,20 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
+    if deg.size == 1:
+        restricted_fit = False
+        lmax = deg[0]
+        order = lmax + 1
+    else:
+        restricted_fit = True
+        lmax = deg.max()
+        order = deg.size
+
     # set up the least squares matrices in transposed form
-    lhs = hermevander(x, deg).T
+    van = hermevander(x, lmax)
+    if restricted_fit:
+        van = van[:, deg]
+    lhs = van.T
     rhs = y.T
     if w is not None:
         w = np.asarray(w) + 0.0
@@ -1528,6 +1548,15 @@ def hermefit(x, y, deg, rcond=None, full=False, w=None):
     c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
     c = (c.T/scl).T
 
+    # Expand c to include non-fitted coefficients which are set to zero
+    if restricted_fit:
+        if c.ndim == 2:
+            cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
+        else:
+            cc = np.zeros(lmax+1, dtype=c.dtype)
+        cc[deg] = c
+        c = cc
+
     # warn on rank reduction
     if rank != order and not full:
         msg = "The fit may be poorly conditioned"
diff --git a/numpy/polynomial/laguerre.py b/numpy/polynomial/laguerre.py
index fffe9e6b6..280e28159 100644
--- a/numpy/polynomial/laguerre.py
+++ b/numpy/polynomial/laguerre.py
@@ -1387,8 +1387,14 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
         y-coordinates of the sample points. Several data sets of sample
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
-    deg : int
-        Degree of the fitting polynomial
+    deg : int or array_like
+        Degree of the fitting polynomial. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included.
+        `deg` may alternatively be a list or array specifying which
+        terms in the Legendre expansion to include in the fit.
+
+        .. versionchanged:: 1.11.0
+        `deg` may be a list specifying which terms to fit
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1485,12 +1491,14 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
     array([ 0.96971004,  2.00193749,  3.00288744])
 
     """
-    order = int(deg) + 1
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
+    deg = np.asarray([deg,], dtype=int).flatten()
 
     # check arguments.
-    if deg < 0:
+    if deg.size < 1:
+        raise TypeError("expected deg to be one or more integers")
+    if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
         raise TypeError("expected 1D vector for x")
@@ -1501,8 +1509,20 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
+    if deg.size == 1:
+        restricted_fit = False
+        lmax = deg[0]
+        order = lmax + 1
+    else:
+        restricted_fit = True
+        lmax = deg.max()
+        order = deg.size
+
     # set up the least squares matrices in transposed form
-    lhs = lagvander(x, deg).T
+    van = lagvander(x, lmax)
+    if restricted_fit:
+        van = van[:, deg]
+    lhs = van.T
     rhs = y.T
     if w is not None:
         w = np.asarray(w) + 0.0
@@ -1530,6 +1550,15 @@ def lagfit(x, y, deg, rcond=None, full=False, w=None):
     c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
     c = (c.T/scl).T
 
+    # Expand c to include non-fitted coefficients which are set to zero
+    if restricted_fit:
+        if c.ndim == 2:
+            cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
+        else:
+            cc = np.zeros(lmax+1, dtype=c.dtype)
+        cc[deg] = c
+        c = cc
+
     # warn on rank reduction
     if rank != order and not full:
         msg = "The fit may be poorly conditioned"
diff --git a/numpy/polynomial/legendre.py b/numpy/polynomial/legendre.py
index c91cb72ec..2035ba6e9 100644
--- a/numpy/polynomial/legendre.py
+++ b/numpy/polynomial/legendre.py
@@ -1418,8 +1418,14 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
         y-coordinates of the sample points. Several data sets of sample
         points sharing the same x-coordinates can be fitted at once by
         passing in a 2D-array that contains one dataset per column.
-    deg : int
-        Degree of the fitting polynomial
+    deg : int or array_like
+        Degree of the fitting polynomial. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included.
+        `deg` may alternatively be a list or array specifying which
+        terms in the Legendre expansion to include in the fit.
+
+        .. versionchanged:: 1.11.0
+        `deg` may be a list specifying which terms to fit
     rcond : float, optional
         Relative condition number of the fit. Singular values smaller than
         this relative to the largest singular value will be ignored. The
@@ -1440,9 +1446,11 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
     Returns
     -------
     coef : ndarray, shape (M,) or (M, K)
-        Legendre coefficients ordered from low to high. If `y` was 2-D,
-        the coefficients for the data in column k  of `y` are in column
-        `k`.
+        Legendre coefficients ordered from low to high. If `y` was
+        2-D, the coefficients for the data in column k of `y` are in
+        column `k`. If `deg` is specified as a list, coefficients for
+        terms not included in the fit are set equal to zero in the
+        returned `coef`.
 
     [residuals, rank, singular_values, rcond] : list
         These values are only returned if `full` = True
@@ -1511,12 +1519,14 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
     --------
 
     """
-    order = int(deg) + 1
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
+    deg = np.asarray([deg,], dtype=int).flatten()
 
     # check arguments.
-    if deg < 0:
+    if deg.size < 1:
+        raise TypeError("expected deg to be one or more integers")
+    if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
         raise TypeError("expected 1D vector for x")
@@ -1527,8 +1537,20 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
+    if deg.size == 1:
+        restricted_fit = False
+        lmax = deg[0]
+        order = lmax + 1
+    else:
+        restricted_fit = True
+        lmax = deg.max()
+        order = deg.size
+
     # set up the least squares matrices in transposed form
-    lhs = legvander(x, deg).T
+    van = legvander(x, lmax)
+    if restricted_fit:
+        van = van[:, deg]
+    lhs = van.T
     rhs = y.T
     if w is not None:
         w = np.asarray(w) + 0.0
@@ -1556,6 +1578,15 @@ def legfit(x, y, deg, rcond=None, full=False, w=None):
     c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
     c = (c.T/scl).T
 
+    # Expand c to include non-fitted coefficients which are set to zero
+    if restricted_fit:
+        if c.ndim == 2:
+            cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
+        else:
+            cc = np.zeros(lmax+1, dtype=c.dtype)
+        cc[deg] = c
+        c = cc
+
     # warn on rank reduction
     if rank != order and not full:
         msg = "The fit may be poorly conditioned"
diff --git a/numpy/polynomial/polynomial.py b/numpy/polynomial/polynomial.py
index 60e339a1d..7c922c11b 100644
--- a/numpy/polynomial/polynomial.py
+++ b/numpy/polynomial/polynomial.py
@@ -1217,8 +1217,14 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
         sharing the same x-coordinates can be (independently) fit with one
         call to `polyfit` by passing in for `y` a 2-D array that contains
         one data set per column.
-    deg : int
-        Degree of the polynomial(s) to be fit.
+    deg : int or array_like
+        Degree of the fitting polynomial. If `deg` is a single integer
+        all terms up to and including the `deg`'th term are included.
+        `deg` may alternatively be a list or array specifying which
+        terms in the Legendre expansion to include in the fit.
+
+        .. versionchanged:: 1.11.0
+        `deg` may be a list specifying which terms to fit
     rcond : float, optional
         Relative condition number of the fit.  Singular values smaller
         than `rcond`, relative to the largest singular value, will be
@@ -1332,12 +1338,14 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
     0.50443316,  0.28853036]), 1.1324274851176597e-014]
 
     """
-    order = int(deg) + 1
     x = np.asarray(x) + 0.0
     y = np.asarray(y) + 0.0
+    deg = np.asarray([deg,], dtype=int).flatten()
 
     # check arguments.
-    if deg < 0:
+    if deg.size < 1:
+        raise TypeError("expected deg to be one or more integers")
+    if deg.min() < 0:
         raise ValueError("expected deg >= 0")
     if x.ndim != 1:
         raise TypeError("expected 1D vector for x")
@@ -1348,8 +1356,20 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
     if len(x) != len(y):
         raise TypeError("expected x and y to have same length")
 
+    if deg.size == 1:
+        restricted_fit = False
+        lmax = deg[0]
+        order = lmax + 1
+    else:
+        restricted_fit = True
+        lmax = deg.max()
+        order = deg.size
+
     # set up the least squares matrices in transposed form
-    lhs = polyvander(x, deg).T
+    van = polyvander(x, lmax)
+    if restricted_fit:
+        van = van[:, deg]
+    lhs = van.T
     rhs = y.T
     if w is not None:
         w = np.asarray(w) + 0.0
@@ -1377,6 +1397,15 @@ def polyfit(x, y, deg, rcond=None, full=False, w=None):
     c, resids, rank, s = la.lstsq(lhs.T/scl, rhs.T, rcond)
     c = (c.T/scl).T
 
+    # Expand c to include non-fitted coefficients which are set to zero
+    if restricted_fit:
+        if c.ndim == 2:
+            cc = np.zeros((lmax+1, c.shape[1]), dtype=c.dtype)
+        else:
+            cc = np.zeros(lmax+1, dtype=c.dtype)
+        cc[deg] = c
+        c = cc
+
     # warn on rank reduction
     if rank != order and not full:
         msg = "The fit may be poorly conditioned"
diff --git a/numpy/polynomial/tests/test_chebyshev.py b/numpy/polynomial/tests/test_chebyshev.py
index a596905f6..8d992c4f0 100644
--- a/numpy/polynomial/tests/test_chebyshev.py
+++ b/numpy/polynomial/tests/test_chebyshev.py
@@ -399,6 +399,9 @@ class TestFitting(TestCase):
         def f(x):
             return x*(x - 1)*(x - 2)
 
+        def f2(x):
+            return x**4 + x**2 + 1
+
         # Test exceptions
         assert_raises(ValueError, cheb.chebfit, [1], [1], -1)
         assert_raises(TypeError, cheb.chebfit, [[1]], [1], 0)
@@ -408,6 +411,9 @@ class TestFitting(TestCase):
         assert_raises(TypeError, cheb.chebfit, [1], [1, 2], 0)
         assert_raises(TypeError, cheb.chebfit, [1], [1], 0, w=[[1]])
         assert_raises(TypeError, cheb.chebfit, [1], [1], 0, w=[1, 1])
+        assert_raises(ValueError, cheb.chebfit, [1], [1], [-1,])
+        assert_raises(ValueError, cheb.chebfit, [1], [1], [2, -1, 6])
+        assert_raises(TypeError, cheb.chebfit, [1], [1], [])
 
         # Test fit
         x = np.linspace(0, 2)
@@ -416,13 +422,25 @@ class TestFitting(TestCase):
         coef3 = cheb.chebfit(x, y, 3)
         assert_equal(len(coef3), 4)
         assert_almost_equal(cheb.chebval(x, coef3), y)
+        coef3 = cheb.chebfit(x, y, [0, 1, 2, 3])
+        assert_equal(len(coef3), 4)
+        assert_almost_equal(cheb.chebval(x, coef3), y)
         #
         coef4 = cheb.chebfit(x, y, 4)
         assert_equal(len(coef4), 5)
         assert_almost_equal(cheb.chebval(x, coef4), y)
+        coef4 = cheb.chebfit(x, y, [0, 1, 2, 3, 4])
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(cheb.chebval(x, coef4), y)
+        # check things still work if deg is not in strict increasing
+        coef4 = cheb.chebfit(x, y, [2, 3, 4, 1, 0])
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(cheb.chebval(x, coef4), y)
         #
         coef2d = cheb.chebfit(x, np.array([y, y]).T, 3)
         assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
+        coef2d = cheb.chebfit(x, np.array([y, y]).T, [0, 1, 2, 3])
+        assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
         # test weighting
         w = np.zeros_like(x)
         yw = y.copy()
@@ -430,13 +448,26 @@ class TestFitting(TestCase):
         y[0::2] = 0
         wcoef3 = cheb.chebfit(x, yw, 3, w=w)
         assert_almost_equal(wcoef3, coef3)
+        wcoef3 = cheb.chebfit(x, yw, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef3, coef3)
         #
         wcoef2d = cheb.chebfit(x, np.array([yw, yw]).T, 3, w=w)
         assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
+        wcoef2d = cheb.chebfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
         # test scaling with complex values x points whose square
         # is zero when summed.
         x = [1, 1j, -1, -1j]
         assert_almost_equal(cheb.chebfit(x, x, 1), [0, 1])
+        assert_almost_equal(cheb.chebfit(x, x, [0, 1]), [0, 1])
+        # test fitting only even polynomials
+        x = np.linspace(-1, 1)
+        y = f2(x)
+        coef1 = cheb.chebfit(x, y, 4)
+        assert_almost_equal(cheb.chebval(x, coef1), y)
+        coef2 = cheb.chebfit(x, y, [0, 2, 4])
+        assert_almost_equal(cheb.chebval(x, coef2), y)
+        assert_almost_equal(coef1, coef2)
 
 
 class TestCompanion(TestCase):
diff --git a/numpy/polynomial/tests/test_classes.py b/numpy/polynomial/tests/test_classes.py
index cd5a54687..a7cf7209c 100644
--- a/numpy/polynomial/tests/test_classes.py
+++ b/numpy/polynomial/tests/test_classes.py
@@ -173,11 +173,18 @@ def check_fit(Poly):
     assert_almost_equal(p(x), y)
     assert_almost_equal(p.domain, d)
     assert_almost_equal(p.window, w)
+    p = Poly.fit(x, y, [0, 1, 2, 3], domain=d, window=w)
+    assert_almost_equal(p(x), y)
+    assert_almost_equal(p.domain, d)
+    assert_almost_equal(p.window, w)
 
     # check with class domain default
     p = Poly.fit(x, y, 3, [])
     assert_equal(p.domain, Poly.domain)
     assert_equal(p.window, Poly.window)
+    p = Poly.fit(x, y, [0, 1, 2, 3], [])
+    assert_equal(p.domain, Poly.domain)
+    assert_equal(p.window, Poly.window)
 
     # check that fit accepts weights.
     w = np.zeros_like(x)
@@ -185,7 +192,9 @@ def check_fit(Poly):
     w[::2] = 1
     p1 = Poly.fit(x[::2], z[::2], 3)
     p2 = Poly.fit(x, z, 3, w=w)
+    p3 = Poly.fit(x, z, [0, 1, 2, 3], w=w)
     assert_almost_equal(p1(x), p2(x))
+    assert_almost_equal(p2(x), p3(x))
 
 
 def check_equal(Poly):
diff --git a/numpy/polynomial/tests/test_hermite.py b/numpy/polynomial/tests/test_hermite.py
index e67625a88..04da72b26 100644
--- a/numpy/polynomial/tests/test_hermite.py
+++ b/numpy/polynomial/tests/test_hermite.py
@@ -387,6 +387,9 @@ class TestFitting(TestCase):
         def f(x):
             return x*(x - 1)*(x - 2)
 
+        def f2(x):
+            return x**4 + x**2 + 1
+
         # Test exceptions
         assert_raises(ValueError, herm.hermfit, [1], [1], -1)
         assert_raises(TypeError, herm.hermfit, [[1]], [1], 0)
@@ -396,6 +399,9 @@ class TestFitting(TestCase):
         assert_raises(TypeError, herm.hermfit, [1], [1, 2], 0)
         assert_raises(TypeError, herm.hermfit, [1], [1], 0, w=[[1]])
         assert_raises(TypeError, herm.hermfit, [1], [1], 0, w=[1, 1])
+        assert_raises(ValueError, herm.hermfit, [1], [1], [-1,])
+        assert_raises(ValueError, herm.hermfit, [1], [1], [2, -1, 6])
+        assert_raises(TypeError, herm.hermfit, [1], [1], [])
 
         # Test fit
         x = np.linspace(0, 2)
@@ -404,13 +410,25 @@ class TestFitting(TestCase):
         coef3 = herm.hermfit(x, y, 3)
         assert_equal(len(coef3), 4)
         assert_almost_equal(herm.hermval(x, coef3), y)
+        coef3 = herm.hermfit(x, y, [0, 1, 2, 3])
+        assert_equal(len(coef3), 4)
+        assert_almost_equal(herm.hermval(x, coef3), y)
         #
         coef4 = herm.hermfit(x, y, 4)
         assert_equal(len(coef4), 5)
         assert_almost_equal(herm.hermval(x, coef4), y)
+        coef4 = herm.hermfit(x, y, [0, 1, 2, 3, 4])
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(herm.hermval(x, coef4), y)
+        # check things still work if deg is not in strict increasing
+        coef4 = herm.hermfit(x, y, [2, 3, 4, 1, 0])
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(herm.hermval(x, coef4), y)
         #
         coef2d = herm.hermfit(x, np.array([y, y]).T, 3)
         assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
+        coef2d = herm.hermfit(x, np.array([y, y]).T, [0, 1, 2, 3])
+        assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
         # test weighting
         w = np.zeros_like(x)
         yw = y.copy()
@@ -418,13 +436,26 @@ class TestFitting(TestCase):
         y[0::2] = 0
         wcoef3 = herm.hermfit(x, yw, 3, w=w)
         assert_almost_equal(wcoef3, coef3)
+        wcoef3 = herm.hermfit(x, yw, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef3, coef3)
         #
         wcoef2d = herm.hermfit(x, np.array([yw, yw]).T, 3, w=w)
         assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
+        wcoef2d = herm.hermfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
         # test scaling with complex values x points whose square
         # is zero when summed.
         x = [1, 1j, -1, -1j]
         assert_almost_equal(herm.hermfit(x, x, 1), [0, .5])
+        assert_almost_equal(herm.hermfit(x, x, [0, 1]), [0, .5])
+        # test fitting only even Legendre polynomials
+        x = np.linspace(-1, 1)
+        y = f2(x)
+        coef1 = herm.hermfit(x, y, 4)
+        assert_almost_equal(herm.hermval(x, coef1), y)
+        coef2 = herm.hermfit(x, y, [0, 2, 4])
+        assert_almost_equal(herm.hermval(x, coef2), y)
+        assert_almost_equal(coef1, coef2)
 
 
 class TestCompanion(TestCase):
diff --git a/numpy/polynomial/tests/test_hermite_e.py b/numpy/polynomial/tests/test_hermite_e.py
index f8601a828..1162502dc 100644
--- a/numpy/polynomial/tests/test_hermite_e.py
+++ b/numpy/polynomial/tests/test_hermite_e.py
@@ -388,6 +388,9 @@ class TestFitting(TestCase):
         def f(x):
             return x*(x - 1)*(x - 2)
 
+        def f2(x):
+            return x**4 + x**2 + 1
+
         # Test exceptions
         assert_raises(ValueError, herme.hermefit, [1], [1], -1)
         assert_raises(TypeError, herme.hermefit, [[1]], [1], 0)
@@ -397,6 +400,9 @@ class TestFitting(TestCase):
         assert_raises(TypeError, herme.hermefit, [1], [1, 2], 0)
         assert_raises(TypeError, herme.hermefit, [1], [1], 0, w=[[1]])
         assert_raises(TypeError, herme.hermefit, [1], [1], 0, w=[1, 1])
+        assert_raises(ValueError, herme.hermefit, [1], [1], [-1,])
+        assert_raises(ValueError, herme.hermefit, [1], [1], [2, -1, 6])
+        assert_raises(TypeError, herme.hermefit, [1], [1], [])
 
         # Test fit
         x = np.linspace(0, 2)
@@ -405,13 +411,25 @@ class TestFitting(TestCase):
         coef3 = herme.hermefit(x, y, 3)
         assert_equal(len(coef3), 4)
         assert_almost_equal(herme.hermeval(x, coef3), y)
+        coef3 = herme.hermefit(x, y, [0, 1, 2, 3])
+        assert_equal(len(coef3), 4)
+        assert_almost_equal(herme.hermeval(x, coef3), y)
         #
         coef4 = herme.hermefit(x, y, 4)
         assert_equal(len(coef4), 5)
         assert_almost_equal(herme.hermeval(x, coef4), y)
+        coef4 = herme.hermefit(x, y, [0, 1, 2, 3, 4])
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(herme.hermeval(x, coef4), y)
+        # check things still work if deg is not in strict increasing
+        coef4 = herme.hermefit(x, y, [2, 3, 4, 1, 0])
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(herme.hermeval(x, coef4), y)
         #
         coef2d = herme.hermefit(x, np.array([y, y]).T, 3)
         assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
+        coef2d = herme.hermefit(x, np.array([y, y]).T, [0, 1, 2, 3])
+        assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
         # test weighting
         w = np.zeros_like(x)
         yw = y.copy()
@@ -419,13 +437,26 @@ class TestFitting(TestCase):
         y[0::2] = 0
         wcoef3 = herme.hermefit(x, yw, 3, w=w)
         assert_almost_equal(wcoef3, coef3)
+        wcoef3 = herme.hermefit(x, yw, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef3, coef3)
         #
         wcoef2d = herme.hermefit(x, np.array([yw, yw]).T, 3, w=w)
         assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
+        wcoef2d = herme.hermefit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
         # test scaling with complex values x points whose square
         # is zero when summed.
         x = [1, 1j, -1, -1j]
         assert_almost_equal(herme.hermefit(x, x, 1), [0, 1])
+        assert_almost_equal(herme.hermefit(x, x, [0, 1]), [0, 1])
+        # test fitting only even Legendre polynomials
+        x = np.linspace(-1, 1)
+        y = f2(x)
+        coef1 = herme.hermefit(x, y, 4)
+        assert_almost_equal(herme.hermeval(x, coef1), y)
+        coef2 = herme.hermefit(x, y, [0, 2, 4])
+        assert_almost_equal(herme.hermeval(x, coef2), y)
+        assert_almost_equal(coef1, coef2)
 
 
 class TestCompanion(TestCase):
diff --git a/numpy/polynomial/tests/test_laguerre.py b/numpy/polynomial/tests/test_laguerre.py
index 1dc57a960..c25476088 100644
--- a/numpy/polynomial/tests/test_laguerre.py
+++ b/numpy/polynomial/tests/test_laguerre.py
@@ -393,6 +393,9 @@ class TestFitting(TestCase):
         assert_raises(TypeError, lag.lagfit, [1], [1, 2], 0)
         assert_raises(TypeError, lag.lagfit, [1], [1], 0, w=[[1]])
         assert_raises(TypeError, lag.lagfit, [1], [1], 0, w=[1, 1])
+        assert_raises(ValueError, lag.lagfit, [1], [1], [-1,])
+        assert_raises(ValueError, lag.lagfit, [1], [1], [2, -1, 6])
+        assert_raises(TypeError, lag.lagfit, [1], [1], [])
 
         # Test fit
         x = np.linspace(0, 2)
@@ -401,13 +404,21 @@ class TestFitting(TestCase):
         coef3 = lag.lagfit(x, y, 3)
         assert_equal(len(coef3), 4)
         assert_almost_equal(lag.lagval(x, coef3), y)
+        coef3 = lag.lagfit(x, y, [0, 1, 2, 3])
+        assert_equal(len(coef3), 4)
+        assert_almost_equal(lag.lagval(x, coef3), y)
         #
         coef4 = lag.lagfit(x, y, 4)
         assert_equal(len(coef4), 5)
         assert_almost_equal(lag.lagval(x, coef4), y)
+        coef4 = lag.lagfit(x, y, [0, 1, 2, 3, 4])
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(lag.lagval(x, coef4), y)
         #
         coef2d = lag.lagfit(x, np.array([y, y]).T, 3)
         assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
+        coef2d = lag.lagfit(x, np.array([y, y]).T, [0, 1, 2, 3])
+        assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
         # test weighting
         w = np.zeros_like(x)
         yw = y.copy()
@@ -415,13 +426,18 @@ class TestFitting(TestCase):
         y[0::2] = 0
         wcoef3 = lag.lagfit(x, yw, 3, w=w)
         assert_almost_equal(wcoef3, coef3)
+        wcoef3 = lag.lagfit(x, yw, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef3, coef3)
         #
         wcoef2d = lag.lagfit(x, np.array([yw, yw]).T, 3, w=w)
         assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
+        wcoef2d = lag.lagfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
         # test scaling with complex values x points whose square
         # is zero when summed.
         x = [1, 1j, -1, -1j]
         assert_almost_equal(lag.lagfit(x, x, 1), [1, -1])
+        assert_almost_equal(lag.lagfit(x, x, [0, 1]), [1, -1])
 
 
 class TestCompanion(TestCase):
diff --git a/numpy/polynomial/tests/test_legendre.py b/numpy/polynomial/tests/test_legendre.py
index 8ac1feb58..9c259d14c 100644
--- a/numpy/polynomial/tests/test_legendre.py
+++ b/numpy/polynomial/tests/test_legendre.py
@@ -388,6 +388,9 @@ class TestFitting(TestCase):
         def f(x):
             return x*(x - 1)*(x - 2)
 
+        def f2(x):
+            return x**4 + x**2 + 1
+
         # Test exceptions
         assert_raises(ValueError, leg.legfit, [1], [1], -1)
         assert_raises(TypeError, leg.legfit, [[1]], [1], 0)
@@ -397,6 +400,9 @@ class TestFitting(TestCase):
         assert_raises(TypeError, leg.legfit, [1], [1, 2], 0)
         assert_raises(TypeError, leg.legfit, [1], [1], 0, w=[[1]])
         assert_raises(TypeError, leg.legfit, [1], [1], 0, w=[1, 1])
+        assert_raises(ValueError, leg.legfit, [1], [1], [-1,])
+        assert_raises(ValueError, leg.legfit, [1], [1], [2, -1, 6])
+        assert_raises(TypeError, leg.legfit, [1], [1], [])
 
         # Test fit
         x = np.linspace(0, 2)
@@ -405,13 +411,25 @@ class TestFitting(TestCase):
         coef3 = leg.legfit(x, y, 3)
         assert_equal(len(coef3), 4)
         assert_almost_equal(leg.legval(x, coef3), y)
+        coef3 = leg.legfit(x, y, [0, 1, 2, 3])
+        assert_equal(len(coef3), 4)
+        assert_almost_equal(leg.legval(x, coef3), y)
         #
         coef4 = leg.legfit(x, y, 4)
         assert_equal(len(coef4), 5)
         assert_almost_equal(leg.legval(x, coef4), y)
+        coef4 = leg.legfit(x, y, [0, 1, 2, 3, 4])
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(leg.legval(x, coef4), y)
+        # check things still work if deg is not in strict increasing
+        coef4 = leg.legfit(x, y, [2, 3, 4, 1, 0])
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(leg.legval(x, coef4), y)
         #
         coef2d = leg.legfit(x, np.array([y, y]).T, 3)
         assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
+        coef2d = leg.legfit(x, np.array([y, y]).T, [0, 1, 2, 3])
+        assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
         # test weighting
         w = np.zeros_like(x)
         yw = y.copy()
@@ -419,13 +437,26 @@ class TestFitting(TestCase):
         y[0::2] = 0
         wcoef3 = leg.legfit(x, yw, 3, w=w)
         assert_almost_equal(wcoef3, coef3)
+        wcoef3 = leg.legfit(x, yw, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef3, coef3)
         #
         wcoef2d = leg.legfit(x, np.array([yw, yw]).T, 3, w=w)
         assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
+        wcoef2d = leg.legfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
         # test scaling with complex values x points whose square
         # is zero when summed.
         x = [1, 1j, -1, -1j]
         assert_almost_equal(leg.legfit(x, x, 1), [0, 1])
+        assert_almost_equal(leg.legfit(x, x, [0, 1]), [0, 1])
+        # test fitting only even Legendre polynomials
+        x = np.linspace(-1, 1)
+        y = f2(x)
+        coef1 = leg.legfit(x, y, 4)
+        assert_almost_equal(leg.legval(x, coef1), y)
+        coef2 = leg.legfit(x, y, [0, 2, 4])
+        assert_almost_equal(leg.legval(x, coef2), y)
+        assert_almost_equal(coef1, coef2)
 
 
 class TestCompanion(TestCase):
diff --git a/numpy/polynomial/tests/test_polynomial.py b/numpy/polynomial/tests/test_polynomial.py
index c806a8497..00a52ebce 100644
--- a/numpy/polynomial/tests/test_polynomial.py
+++ b/numpy/polynomial/tests/test_polynomial.py
@@ -419,6 +419,9 @@ class TestMisc(TestCase):
         def f(x):
             return x*(x - 1)*(x - 2)
 
+        def f2(x):
+            return x**4 + x**2 + 1
+
         # Test exceptions
         assert_raises(ValueError, poly.polyfit, [1], [1], -1)
         assert_raises(TypeError, poly.polyfit, [[1]], [1], 0)
@@ -428,6 +431,9 @@ class TestMisc(TestCase):
         assert_raises(TypeError, poly.polyfit, [1], [1, 2], 0)
         assert_raises(TypeError, poly.polyfit, [1], [1], 0, w=[[1]])
         assert_raises(TypeError, poly.polyfit, [1], [1], 0, w=[1, 1])
+        assert_raises(ValueError, poly.polyfit, [1], [1], [-1,])
+        assert_raises(ValueError, poly.polyfit, [1], [1], [2, -1, 6])
+        assert_raises(TypeError, poly.polyfit, [1], [1], [])
 
         # Test fit
         x = np.linspace(0, 2)
@@ -436,13 +442,21 @@ class TestMisc(TestCase):
         coef3 = poly.polyfit(x, y, 3)
         assert_equal(len(coef3), 4)
         assert_almost_equal(poly.polyval(x, coef3), y)
+        coef3 = poly.polyfit(x, y, [0, 1, 2, 3])
+        assert_equal(len(coef3), 4)
+        assert_almost_equal(poly.polyval(x, coef3), y)
         #
         coef4 = poly.polyfit(x, y, 4)
         assert_equal(len(coef4), 5)
         assert_almost_equal(poly.polyval(x, coef4), y)
+        coef4 = poly.polyfit(x, y, [0, 1, 2, 3, 4])
+        assert_equal(len(coef4), 5)
+        assert_almost_equal(poly.polyval(x, coef4), y)
         #
         coef2d = poly.polyfit(x, np.array([y, y]).T, 3)
         assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
+        coef2d = poly.polyfit(x, np.array([y, y]).T, [0, 1, 2, 3])
+        assert_almost_equal(coef2d, np.array([coef3, coef3]).T)
         # test weighting
         w = np.zeros_like(x)
         yw = y.copy()
@@ -450,13 +464,26 @@ class TestMisc(TestCase):
         yw[0::2] = 0
         wcoef3 = poly.polyfit(x, yw, 3, w=w)
         assert_almost_equal(wcoef3, coef3)
+        wcoef3 = poly.polyfit(x, yw, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef3, coef3)
         #
         wcoef2d = poly.polyfit(x, np.array([yw, yw]).T, 3, w=w)
         assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
+        wcoef2d = poly.polyfit(x, np.array([yw, yw]).T, [0, 1, 2, 3], w=w)
+        assert_almost_equal(wcoef2d, np.array([coef3, coef3]).T)
         # test scaling with complex values x points whose square
         # is zero when summed.
         x = [1, 1j, -1, -1j]
         assert_almost_equal(poly.polyfit(x, x, 1), [0, 1])
+        assert_almost_equal(poly.polyfit(x, x, [0, 1]), [0, 1])
+        # test fitting only even Polyendre polynomials
+        x = np.linspace(-1, 1)
+        y = f2(x)
+        coef1 = poly.polyfit(x, y, 4)
+        assert_almost_equal(poly.polyval(x, coef1), y)
+        coef2 = poly.polyfit(x, y, [0, 2, 4])
+        assert_almost_equal(poly.polyval(x, coef2), y)
+        assert_almost_equal(coef1, coef2)
 
     def test_polytrim(self):
         coef = [2, -1, 1, 0]