See issue #7780: https://github.com/numpy/numpy/issues/7780 Would be good if the numpy.polynomial fit methods (e.g. numpy.polynomial.legendre.Legendre.fit) could return the covariance matrix of the fitted parameters, in the same way as numpy.polyfit. Let me know if there are any suggestions/concerns. Otherwise I’m happy to add the code and update the docs. Tom