<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Introduction to Machine Learning with scikit-learn
Lab 3
Step1:  Generating the data
As we did during the classification lab, we generate some data points uniformly using np.random.rand. In this example, we will generate a set of data points $x \in [0, 1]$ and their images $y$ under the function $f$ defined by
$$ y = f(x) = \cos \left( \dfrac{3\pi}{2} x \right)$$
which will be modeled by my_function.
We also add some noise following a Gaussian distribution (np.random.randn) centered in 0.
Step2:  Let us also generate the testing points needed to plot the result of our prediction. In this case, we generate 100 points uniformly spread between 0 and 1 using np.linspace.
Step3:  With this, we can plot $f$ and the generated points.
Step4:  Note that the samples are not on the function curve because of the additive noise.
Linear regression
Recall that the linear regression with a squared loss aims at finding the solution of the following optimization problem
Step5:  Defining polynomial features
The previous example is a simple linear regression and as we can see, the trained model fits the samples (and $f$) quite poorly. As we have seen during the lectures, we can use polynomial features to obtain a more accurate model.
Formally, it consists in defining a polynomial kernel $\phi$ such that of degree $d$ such that, for each $x \in [0, 1]$, we have
$$\phi(x) = \left[1, x, x^2, \dots, x^d\right]$$
In sklearn, this can be achieved with the sklearn.preprocessing.PolynomialFeatures function
Step6:  Note that here we use np.newaxis to create an array of single elements. Applying polynomial_features.fit_transform to X will not lead to the desired result (try it by yourself).
Step7:  So far, we have just created a feature of degree 1, which won't change much the data we have. However, the shape of the polynomial features has changed
Step8:  Let's take a closer look at what is inside polynomial_X.
Step9:  By default, PolynomialFeatures adds a constant term equal to 1. It corresponds to the intercept (or degree 0) of the model we are training.
Let's now define another linear regression and train it on the augmented features $\phi(x) = \left[1, x, x^2\right]$.
Step10:  And display the result
Step11:  (Note that we need to apply $\phi$ to X_test because my_regression has been trained on polynomial features!).
As we can see, a second order kernel gives a much better fit than a simple linear regression.
Putting everything together with pipelines
So far, our regression consists of 2 different steps
Step12:  Regularizing the linear regression
The Linear regression we have been working with solves the following problem
Step13:  Support vector regression
Support vector regression (sklearn.svm.SVR in sklearn) is another way to train a regression model. The difference with the linear regression with least squares we have used previously is in the loss function
Step14:  Note that chosing $\epsilon$ with an order of magnitude close to the estimated noise makes sense because we would not penalize errors smaller than the noise.
Step15:  Applying regression to standard datasets
As for classification, sklearn proposes standard datasets to evaluate models. These datasets can be found here
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning) 
Explanation: Introduction to Machine Learning with scikit-learn
Lab 3: Regression
The goal of this lab session is to discover a few regression tools from scikit-learn. As for the classification lab, we will start with generated data to have a more visual idea of the results.
End of explanation
import numpy as np
np.random.seed(0)
n_samples = 30
def my_function(X):
    return np.cos(1.5 * np.pi * X)
    
X = np.random.rand(n_samples)
y = my_function(X) + np.random.randn(n_samples) * 0.1
X.shape
Explanation: Generating the data
As we did during the classification lab, we generate some data points uniformly using np.random.rand. In this example, we will generate a set of data points $x \in [0, 1]$ and their images $y$ under the function $f$ defined by
$$ y = f(x) = \cos \left( \dfrac{3\pi}{2} x \right)$$
which will be modeled by my_function.
We also add some noise following a Gaussian distribution (np.random.randn) centered in 0.
End of explanation
X_test = np.linspace(0, 1, 100)
X_test
Explanation: Let us also generate the testing points needed to plot the result of our prediction. In this case, we generate 100 points uniformly spread between 0 and 1 using np.linspace.
End of explanation
import matplotlib.pyplot as plt
plt.plot(X_test, my_function(X_test), label="True function")
plt.scatter(X, y, label="Samples")
plt.legend(loc="best")
Explanation: With this, we can plot $f$ and the generated points.
End of explanation
from sklearn.linear_model import LinearRegression
linear_regression = LinearRegression()
my_regression = linear_regression.fit(X[:, np.newaxis], y)
X_test = np.linspace(0, 1, 100)
plt.plot(X_test, my_regression.predict(X_test[:, np.newaxis]), label="Model")
# my_regression.predict(X_test[:, np.newaxis] = y_test (the prediction of X_test)
plt.plot(X_test, my_function(X_test), label="True function")
plt.scatter(X, y, label="Samples")
plt.xlabel("x")
plt.ylabel("y")
plt.xlim((0, 1))
plt.ylim((-2, 2))
plt.legend(loc="best")
Explanation: Note that the samples are not on the function curve because of the additive noise.
Linear regression
Recall that the linear regression with a squared loss aims at finding the solution of the following optimization problem:
$$\min_w \sum_i \left(w^T \phi\left(x^{(i)}\right) - y^{(i)}\right)^2$$
Linear regression is implemented in sklearn.linear_model.LinearRegression.
End of explanation
from sklearn.preprocessing import PolynomialFeatures
polynomial_features = PolynomialFeatures(degree=1)
polynomial_X = polynomial_features.fit_transform(X[:, np.newaxis])
Explanation: Defining polynomial features
The previous example is a simple linear regression and as we can see, the trained model fits the samples (and $f$) quite poorly. As we have seen during the lectures, we can use polynomial features to obtain a more accurate model.
Formally, it consists in defining a polynomial kernel $\phi$ such that of degree $d$ such that, for each $x \in [0, 1]$, we have
$$\phi(x) = \left[1, x, x^2, \dots, x^d\right]$$
In sklearn, this can be achieved with the sklearn.preprocessing.PolynomialFeatures function:
End of explanation
X[:, np.newaxis]
Explanation: Note that here we use np.newaxis to create an array of single elements. Applying polynomial_features.fit_transform to X will not lead to the desired result (try it by yourself).
End of explanation
polynomial_X.shape
X[:, np.newaxis].shape
Explanation: So far, we have just created a feature of degree 1, which won't change much the data we have. However, the shape of the polynomial features has changed:
End of explanation
polynomial_X
Explanation: Let's take a closer look at what is inside polynomial_X.
End of explanation
linear_regression = LinearRegression()
polynomial_features = PolynomialFeatures(degree=2)
polynomial_X = polynomial_features.fit_transform(X[:, np.newaxis])
my_regression = linear_regression.fit(polynomial_X, y)
Explanation: By default, PolynomialFeatures adds a constant term equal to 1. It corresponds to the intercept (or degree 0) of the model we are training.
Let's now define another linear regression and train it on the augmented features $\phi(x) = \left[1, x, x^2\right]$.
End of explanation
X_test = np.linspace(0, 1, 100)
polynomial_X_test = polynomial_features.fit_transform(X_test[:, np.newaxis])
plt.plot(X_test, my_regression.predict(polynomial_X_test), label="Model")
plt.plot(X_test, my_function(X_test), label="True function")
plt.scatter(X, y, label="Samples")
plt.xlabel("x")
plt.ylabel("y")
plt.xlim((0, 1))
plt.ylim((-2, 2))
plt.legend(loc="best")
Explanation: And display the result:
End of explanation
from sklearn.pipeline import Pipeline
from sklearn.cross_validation import cross_val_score
degrees = [1, 4, 15]
plt.figure(figsize=(14, 5))
for i in range(len(degrees)):
    ax = plt.subplot(1, len(degrees), i + 1)
    plt.setp(ax, xticks=(), yticks=())
    polynomial_features = PolynomialFeatures(degree=degrees[i])
    linear_regression = LinearRegression()
    pipeline = Pipeline([("polynomial_features", polynomial_features),
                         ("linear_regression", linear_regression)])
    pipeline.fit(X[:, np.newaxis], y)
    # Evaluate the models using crossvalidation
    scores = cross_val_score(pipeline, X[:, np.newaxis], y,
                             scoring="mean_squared_error", cv=10)
    X_test = np.linspace(0, 1, 100)
    plt.plot(X_test, pipeline.predict(X_test[:, np.newaxis]), label="Model")
    plt.plot(X_test, my_function(X_test), label="True function")
    plt.scatter(X, y, label="Samples")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.xlim((0, 1))
    plt.ylim((-2, 2))
    plt.legend(loc="best")
    plt.title("Degree {}\nMSE = {:.2e}(+/- {:.2e})".format(
        degrees[i], -scores.mean(), scores.std()))
plt.show()
Explanation: (Note that we need to apply $\phi$ to X_test because my_regression has been trained on polynomial features!).
As we can see, a second order kernel gives a much better fit than a simple linear regression.
Putting everything together with pipelines
So far, our regression consists of 2 different steps:
- A polynomial transformation on X which returns polynomial_X
- A linear regression on polynomial_X and y
It is possible to write this in a more compact way with a Pipeline object, as the following:
End of explanation
from sklearn.linear_model import Ridge
polynomial_features = PolynomialFeatures(degree=15)
polynomial_X = polynomial_features.fit_transform(X[:, np.newaxis])
ridge_regression = Ridge(alpha = .1)
my_regression = ridge_regression.fit(polynomial_X, y)
X_test = np.linspace(0, 1, 100)
X_test_polynomial = polynomial_features.fit_transform(X_test[:, np.newaxis])
plt.plot(X_test, my_regression.predict(X_test_polynomial), label="Model")
plt.plot(X_test, my_function(X_test), label="True function")
plt.scatter(X, y, label="Samples")
plt.xlabel("x")
plt.ylabel("y")
plt.xlim((0, 1))
plt.ylim((-2, 2))
plt.legend(loc="best")
plt.title("Ridge regression")
Explanation: Regularizing the linear regression
The Linear regression we have been working with solves the following problem:
$$\min_w \sum_i \left(w^T \phi\left(x^{(i)}\right) - y^{(i)}\right)^2$$
We can add a regularization term to the optimization problem in order to penalise big coefficients in $w$ and avoid overfitting as a consequence:
$$\min_w \sum_i \left(w^T \phi\left(x^{(i)}\right) - y^{(i)}\right)^2 + \alpha ||w||_2^2$$
Note that this leads to adding an extra parameter, $\alpha$, to the model. As we will see in the following, this parameter has a big impact on how much we overfit or underfit the training data.
End of explanation
from sklearn.svm import SVR
epsilon = .1
sv_regression = SVR(epsilon = epsilon)
my_svr = sv_regression.fit(polynomial_X, y)
Explanation: Support vector regression
Support vector regression (sklearn.svm.SVR in sklearn) is another way to train a regression model. The difference with the linear regression with least squares we have used previously is in the loss function:
$$\min_w \sum_i \ell \left(w^T \phi\left(x^{(i)}\right), y^{(i)}\right)$$
For ordinary least squares, we had 
$$\ell \left(w^T \phi\left(x^{(i)}\right), y^{(i)}\right) = \left(w^T \phi\left(x^{(i)}\right) - y^{(i)}\right)^2.$$ For SVR, we have
$$\ell \left(w^T \phi\left(x^{(i)}\right), y^{(i)}\right) = 
\begin{cases}
    0,& \text{if } -\epsilon \leq w^T \phi\left(x\right) - y \leq +\epsilon\
    \left| w^T \phi\left(x\right) - y \right| - \epsilon,              & \text{otherwise}
\end{cases}.
$$
Having this in mind, we can apply SVR by chosing an apprioriate value for $\epsilon$.
End of explanation
plt.plot(X_test, my_regression.predict(X_test_polynomial), label="Model")
plt.plot(X_test, my_function(X_test), label="True function")
plt.scatter(X, y, label="Samples")
plt.xlabel("x")
plt.ylabel("y")
plt.xlim((0, 1))
plt.ylim((-2, 2))
plt.legend(loc="best")
plt.title("Support vector regression")
Explanation: Note that chosing $\epsilon$ with an order of magnitude close to the estimated noise makes sense because we would not penalize errors smaller than the noise.
End of explanation
from sklearn import datasets
boston = datasets.load_boston()
print boston.DESCR
boston.feature_names
print boston.data.shape, boston.target.shape
polynomial_features = PolynomialFeatures(degree=3)
polynomial_X = polynomial_features.fit_transform(boston.data)
alpha = .1
ridge_regression = Ridge(alpha = alpha)
my_regression = ridge_regression.fit(polynomial_X, boston.target)
predicted_prices = my_regression.predict(polynomial_X)
predicted_prices
from sklearn.cross_validation import train_test_split
features_train, features_test, labels_train, labels_test = train_test_split(polynomial_X, boston.target, test_size = 0.5)
print features_train.shape, features_test.shape, labels_train.shape, labels_test.shape
from sklearn import metrics
predicted_labels_train = my_regression.predict(features_train)
print "MSE on the train set: ", metrics.mean_squared_error(labels_train, predicted_labels_train)
predicted_labels_test = my_regression.predict(features_test)
print "MSE on the test set: ", metrics.mean_squared_error(labels_test, predicted_labels_test)
Explanation: Applying regression to standard datasets
As for classification, sklearn proposes standard datasets to evaluate models. These datasets can be found here: http://scikit-learn.org/stable/datasets/.
Two dataset are suited for regression problem:
- The Boston house-prices dataset
- diabetes dataset
End of explanation
<END_TASK>