title: Multivariate polynomial regression with numpy
body: <p>I have many samples <code>(y_i, (a_i, b_i, c_i))</code> where <code>y</code> is presumed to vary as a polynomial in <code>a,b,c</code>  up to a certain degree. For example for a given set of data and degree 2 I might produce the model </p>
  
  <p><code>y = a^2 + 2ab - 3cb + c^2 +.5ac</code></p>
  
  <p>This can be done using least squares and is a slight extension of numpy's polyfit routine. Is there a standard implementation somewhere in the Python ecosystem?</p>