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