Alternative solution for the original gist for fitting plane/surface to a set of data points in Python with optional using of weights on data points. Polynomial features are generated with PolynomialFeature form scikit-learn and must not be coded manually. Additional information are in my answer to the original gist