Polynomial regression using scikit-learn

Question

I am trying to use scikit-learn for polynomial regression. From what I read polynomial regression is a special case of linear regression. I was hopping that maybe one of scikit's generalized linear models can be parameterised to fit higher order polynomials but I see no option for doing that.

I did manage to use a Support Vector Regressor with a poly kernel. That worked well with a subset of my data, but it takes much to long to fit larger data sets so I still need to find something faster (even if trading some precision).

Am I missing something obvious here?

Answer 1

Given data $\mathbf{x}$, a column vector, and $\mathbf{y}$, the target vector, you can perform polynomial regression by appending polynomials of $\mathbf{x}$. For example, consider if

$$ \mathbf{x} = \begin{bmatrix} 2 \\[0.3em] -1 \\[0.3em] \frac{1}{3} \end{bmatrix}$$

Using just this vector in linear regression implies the model:

$$ y = \alpha_1 x $$

We can add columns that are powers of the vector above, which represent adding polynomials to the regression. Below we show this for polynomials up to power 3: $$ X = \begin{bmatrix} 2 & 4 & 8 \\[0.3em] -1 & 1 & -1 \\[0.3em] \frac{1}{3} & \frac{1}{3^2} & \frac{1}{3^3} \end{bmatrix}$$

This is our new data matrix that we use in sklearn's linear regression, and it represents the model:

$$ y = \alpha_1 x + \alpha_2x^2 + \alpha_3x^3$$

---

Note that I did not add a constant vector of $1$'s, as sklearn will automatically include this.

Answer 2

Theory

Polynomial regression is a special case of linear regression. With the main idea of how do you select your features. Looking at the multivariate regression with 2 variables: x1 and x2. Linear regression will look like this: y = a1 * x1 + a2 * x2.

Now you want to have a polynomial regression (let's make 2 degree polynomial). We will create a few additional features: x1*x2, x1^2 and x2^2. So we will get your 'linear regression':

y = a1 * x1 + a2 * x2 + a3 * x1*x2 + a4 * x1^2 + a5 * x2^2

This nicely shows an important concept curse of dimensionality, because the number of new features grows much faster than linearly with the growth of degree of polynomial.

Practice with scikit-learn

You do not need to do all this in scikit. Polynomial regression is already available there (in 0.15 version. Check how to update it here).

from sklearn.preprocessing import PolynomialFeatures
from sklearn import linear_model

X = [[0.44, 0.68], [0.99, 0.23]]
vector = [109.85, 155.72]
predict= [[0.49, 0.18]]
#Edit: added second square bracket above to fix the ValueError problem

poly = PolynomialFeatures(degree=2)
X_ = poly.fit_transform(X)
predict_ = poly.fit_transform(predict)

clf = linear_model.LinearRegression()
clf.fit(X_, vector)
print clf.predict(predict_)

Answer 3

In case you are using a multivariate regression and not just a univariate regression, do not forget the cross terms. For instance if you have two variables $x_1$ and $x_2$, and you want polynomials up to power 2, you should use $y = a_1x_1 + a_2x_2 + a_3x_1^2 + a_4x_2^2 + a_5x_1x_2$ where the last term ($a_5x_1x_2$) is the one I am talking about.

Answer 4

Using a similar approach to @Cam.Davidson.Pilon, I wrote a couple functions to help demo this approach in Python. It can be expanded by adding more terms in the np.concatenate vectors. The output for the y_pred would not change, but getting the coefficients, regr.coef_[0][2], would need to be included.

from sklearn import linear_model
import numpy as np
from matplotlib import pyplot as plt
def regress_2nd_order(x, y):
    # The x values are transformed into the second order matrix
    X = np.concatenate([x, x**2], axis=1)
    # Use model to fit to the data 
    regr.fit(X, y)
    # Extract the values of interest for forming the equation...
    # y_pred = c0 + c1 * x + c2 * x^2 
    c0 = round(regr.intercept_[0],2)
    c1 = round(regr.coef_[0][0],2)
    c2 = round(regr.coef_[0][1],4)
    # c3 = round(regr.coef_[0][2],4)
    # Another way to get this is using the `regr.predict` function
    # Predict will calculate the values for you
    y_pred = regr.predict(X)
    # Print the intercept (c0) and the corresponding coefficients
    print('Scikit learn - \nEquation: %.2f + %.2f*T + %.5f*T^2' % (c0,c1,c2) )
    #  return y_pred, (c0, c1, c2, c3)
    return y_pred, (c0, c1, c2)
def socket_temp(x, y):
    # Include only the predicted y-vector
    y_pred = regress_2nd_order(x, y)[0]
    fig, ax = plt.subplots(1,1, figsize=(10,10))
    ax.plot(x, y_pred, ls='', marker='*', label='predicted')
    ax.plot(x, y, ls='', marker='.', label='iddt')
    ax.legend()

For my use, I had temperature as my y variable, and a dynamic supply current (IDDT) output as my x variable. Here are the values I used for x and y and the output vector, y_pred:

    'iddt'  'temp'  'y_pred'
0   627.23  20.0    20.42
1   627.38  19.9    20.52
2   612.39  9.7     10.32
3   612.48  9.7     10.38
4   597.73  -0.1    0.19
5   597.78  -0.2    0.22
6   583.24  -10.1   -9.99
7   583.24  -10.0   -9.99
8   569.22  -19.9   -19.99
9   568.97  -19.9   -20.17
10  555.62  -29.8   -29.83
11  555.99  -29.8   -29.56
12  541.15  -40.1   -40.45
13  541.50  -40.1   -40.19
14  639.88  30.3    28.89
15  640.11  30.2    29.05
16  656.61  40.2    39.92
17  656.62  40.3    39.92
18  672.73  50.4    50.33
19  672.68  50.4    50.30
20  689.56  60.6    61.01
21  689.26  60.6    60.82
22  705.30  70.8    70.79
23  705.22  70.8    70.75
24  722.27  81.0    81.14
25  722.24  81.0    81.12