Return polyfit coefficients (or SymPy expression)

[sigvaldm]

Today polyfit() doesn't return fitting coefficient like NumPy's does, but NumPy's is only 1D. Now that localreg is multivariate, it would be nice to make it more similar to NumPy's for replacement. The ordering of the coefficients are more complicated for the multivariate case, but perhaps there could also be a class or something that could return a SymPy expression of the fit?