To follow up on these comments, an enhancement proposal (https://github.com/numpy/numpy/pull/20889), which would enable numpy.polynomial.polynomial.polyfit to return the covariance matrix, appears to be moving forward. This is great news, and addresses item (1) of my previous message. After thinking a bit more about item (2), is it too late to change the default setting to "[]" for the domain parameter in the convenience class `fit` methods? A user who is aware of the scaling and wants to preserve the rescaled domain information would still be able to do that. The rest of us could remain blissfully unaware of what is happening under the hood. This option would also facilitate adding a covariance matrix output to these fit methods, which I would also recommend. The covariance matrix will also depend on the scaling, and I assume that most users would want to know its element values for the original scale. Making these changes to the API would enable users to replace numpy.polynomial.polynomial.polyfit and numpy.polynomial.polynomial.polval (and their equivalents for other polynomial classes) with the more sophisticated object-oriented interface more readily. I like many of the features of the convenience classes, but the lack of parameter uncertainty output is a deal-breaker. The overhead of managing the rescaling up front is also a barrier for introductory data analysis applications. Either way, as the example in my previous email demonstrates, the current behavior produces output that is unnecessarily confusing. Even if the API remains unchanged, the way polynomial instances are represented to the user needs improvement.