Allow returning cov

See numpy#7880

numpy/polynomial/hermite.py

@@ -1272,7 +1272,7 @@ def hermvander3d(x, y, z, deg):
     return pu._vander_nd_flat((hermvander, hermvander, hermvander), (x, y, z), deg)
 
 
-def hermfit(x, y, deg, rcond=None, full=False, w=None):
+def hermfit(x, y, deg, rcond=None, full=False, w=None, cov=False):
     """
     Least squares fit of Hermite series to data.
 
@@ -1315,7 +1315,10 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
         chosen so that the errors of the products ``w[i]*y[i]`` all have the
         same variance.  When using inverse-variance weighting, use
         ``w[i] = 1/sigma(y[i])``.  The default value is None.
-
+   cov : bool, optional
+        Return the estimate and the covariance matrix of the estimate
+        If full is True, then cov is not returned.
+     
     Returns
     -------
     coef : ndarray, shape (M,) or (M, K)
@@ -1330,7 +1333,13 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
         - rank -- the numerical rank of the scaled Vandermonde matrix
         - singular_values -- singular values of the scaled Vandermonde matrix
         - rcond -- value of `rcond`.
-
+    V : ndarray, shape (M,M) or (M,M,K)
+          Present only if `full` = False and `cov`=True.  The covariance
+          matrix of the polynomial coefficient estimates.  The diagonal of
+          this matrix are the variance estimates for each coefficient.  If y
+          is a 2-D array, then the covariance matrix for the `k`-th data set
+          are in ``V[:,:,k]``
+   
         For more details, see `numpy.linalg.lstsq`.
 
     Warns
@@ -1401,7 +1410,7 @@ def hermfit(x, y, deg, rcond=None, full=False, w=None):
     array([1.0218, 1.9986, 2.9999]) # may vary
 
     """
-    return pu._fit(hermvander, x, y, deg, rcond, full, w)
+    return pu._fit(hermvander, x, y, deg, rcond, full, w, cov)
 
 
 def hermcompanion(c):