<jupyter_start><jupyter_text>Under- and overfitting, model selection PreliminariesIn the first set of exercises you had to implement the training and evaluation of the linear regression and $k$-NN methods from scratch in order to practice your `numpy` skills. From this set of exercises onward, you can use the implementations provided in `scikit-learn` or other higher-level libraries. We start this set of exercises by demonstrating some of the features of `scikit-learn`.For example, implementation of linear regression model fitting with an analytical solution for the parameters is provided by the class `sklearn.linar_model.LinearRegression`. You can train a linear regression model in the following way:<jupyter_code>import numpy as np
from sklearn import datasets, linear_model

# load the diabetes dataset
diabetes = datasets.load_diabetes()

# use only one feature
X = diabetes.data[:, np.newaxis, 2]
y = diabetes.target

# split the data into training/testing sets
X_train = X[:-20]
X_test = X[-20:]

# split the targets into training/testing sets
y_train = y[:-20]
y_test = y[-20:]

# create linear regression object
model = linear_model.LinearRegression()

# train the model using the training dataset
model.fit(X_train, y_train)<jupyter_output><empty_output><jupyter_text>Let's visualize the training dataset and the learned regression model.<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt

fig = plt.figure()
plt.plot(X_train, y_train, 'r.', markersize=12)
X_edge = np.array([np.min(X_train, 0), np.max(X_train, 0)])
plt.plot(X_edge, model.predict(X_edge), 'b-')
plt.legend(('Data', 'Linear regression'), loc='lower right')
plt.title('Linear regression')
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()<jupyter_output><empty_output><jupyter_text>Once trained, the model can be used to make predictions on the test data:<jupyter_code># Make predictions using the testing dataset
prediction = model.predict(X_test)<jupyter_output><empty_output><jupyter_text>The next step (not shown here) is to evaluate the performance of the trained model.Note that the `scikit-learn` interface works by first initializing an object from the class that implements the machine learning model (linear regression in this case) and then fitting the initialized model using the data in the training set. Finally, the trained (fitted) model can be used to make predictions on unseen data. In fact, all models implemented in this library follow the same *initialize-fit-predict* programming interface. For example, a $k$-NN classifier can be trained in the following way:<jupyter_code>from sklearn.model_selection import train_test_split
from sklearn import datasets, neighbors

breast_cancer = datasets.load_breast_cancer()

X = breast_cancer.data
y = breast_cancer.target

# make use of the train_test_split() utility function instead
# of manually dividing the data
X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=40)

# initialize a 3-NN classifier
model = neighbors.KNeighborsClassifier(n_neighbors=3)

# train the model using the training dataset
model.fit(X_train, y_train)

# make predictions using the testing dataset
prediction = model.predict(X_test)<jupyter_output><empty_output><jupyter_text>Note that the features in the breast cancer dataset have different scales (some have on average very small absolute values, and some very large), which means that the distance metric used by $k$-NN will me dominated by the features with large values. You can use any of the number of feature transformation methods implemented in `scikit-learn` to scale the features. For example, you can use the `sklearn.preprocessing.StandardScaler` method to transform all features to a have a zero mean and unit variance:<jupyter_code>from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
scaler.fit(X_train)

X_train_scaled = scaler.transform(X_train)<jupyter_output><empty_output><jupyter_text>The scaler has its own parameters which are the means and standard deviations of the features estimated from the training set. If you train a model with the scaled features, you will have to remember to also apply the scaling transformation every time you make a prediction on new unseen and unscaled data. This is somewhat prone to error. One option for making the code more robust is to create a processing pipeline that includes the scaling and $k$-NN models in a sequence:<jupyter_code>from sklearn.pipeline import Pipeline

knn = neighbors.KNeighborsClassifier(n_neighbors=3)

model = Pipeline([
                 ("scaler", scaler),
                 ("knn", knn)
                ])

# train the model using the training dataset
model.fit(X_train, y_train)

# make predictions using the testing dataset
prediction = model.predict(X_test)<jupyter_output><empty_output><jupyter_text>If you are curious, more information about the design of the `scikit-learn` application programming interface (API) can be found [in this paper](https://arxiv.org/pdf/1309.0238.pdf).   Exercises Bias-variance decompositionShow that the mean squared error of the estimate of a parameter can be decomposed into an expression that includes both the bias and variance (Eq. 5.53-5.54 in "Deep learning" by Goodfellow et al.).***Answer***: $MSE = \mathbb{E} [(\hat{\theta}_{m} - \theta)^2] $let $\mu = \mathbb{E} [\hat{\theta}_{m}]$add and substract $\mu$$\mathbb{E} [(\hat{\theta}_{m} - \theta)^2] = \mathbb{E}[((\hat{\theta}_{m} - \mu ) + ( \mu - \theta))^2]$$\mathbb{E} [(\hat{\theta}_{m} - \theta)^2] = \mathbb{E}[(\hat{\theta}_{m} - \mu )^2 + 2(\hat{\theta}_{m} - \mu)(\mu - \theta)+ (\mu - \theta)^2]$$\mathbb{E}(\hat{\theta}_{m} - \mu ) = 0$$\mathbb{E} [(\hat{\theta}_{m} - \theta)^2] = \mathbb{E}[(\hat{\theta}_{m} - \mu )^2 + (\mu - \theta)^2]$$\mathbb{E} [(\hat{\theta}_{m} - \theta)^2] = \mathbb{E}[(\hat{\theta}_{m} - \mathbb{E} [\hat{\theta}_{m}] )^2 + (\mathbb{E} [\hat{\theta}_{m}] - \theta)^2]$Bias of an estimator is defined as: $bias(\boldsymbol{\hat{\theta}_{m}}) = \mathbb{E} (\boldsymbol{\hat{\theta}_{m}}) - \boldsymbol{\theta}$Variance is defined as: $var(\boldsymbol{\hat{\theta}_{m}}) = \mathbb{E}(\hat{\theta}_{m} - \mathbb{E} [\hat{\theta}_{m}] )^2$$MSE = bias(\boldsymbol{\hat{\theta}_{m}})^2 + Var ({\hat{\theta}_{m}}) $  Polynomial regressionFor this exercise we will be using generated data to better show the effects of the different polynomial orders.The data is created using the make_polynomial_regression function.<jupyter_code>%matplotlib inline

def generate_dataset(n=100, degree=1, noise=1, factors=None):
    # Generates a dataset by adding random noise to a randomly
    # generated polynomial function.
    
    x = np.random.uniform(low=-1, high=1, size=n)
    
    factors = np.random.uniform(0, 10, degree+1)
    
    y = np.zeros(x.shape)
    
    for idx in range(degree+1):
        y += factors[idx] * (x ** idx)

    # add noise
    y += np.random.normal(-noise, noise, n)
    
    return x, y

# load generated data
np.random.seed(0)

X, y = generate_dataset(n=100, degree=4, noise=1.5)

plt.plot(X, y, 'r.', markersize=12)<jupyter_output><empty_output><jupyter_text>Implement polynomial regression using the `sklearn.preprocessing.PolynomialFeatures` transformation. Using the `sklearn.grid_search.GridSearchCV` class, perform a grid search of the polynomial order hyperparameter space with cross-validation and report the performance on an independent test set.Plot a learning curve that show the validation accuracy as a function of the polynomial order. Which models have a high bias, and which models have high variance? Motivate your answer.    Repeat this experiment, this time using the diabetes dataset instead of the generated data.<jupyter_code>from sklearn.preprocessing import PolynomialFeatures
from sklearn.model_selection import GridSearchCV, cross_val_score
from sklearn.metrics import accuracy_score
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline

# independent test set
X_test, y_test = generate_dataset(n=100, degree=4, noise=1.5)
# reshape because otherwise X is an 1D array instead of 2D:
X_test = X_test.reshape(-1, 1)
X = X.reshape(-1, 1)
x_plot = np.linspace(-1, 1, 100)

# polynomial order hyperparameter for grid search
max_order = 5
parameters = {'order':np.arange(1,max_order+1)}

x_plot = np.linspace(-1, 1, 100)

# create matrix versions of these arrays

X_plot = x_plot[:, np.newaxis]


plt.scatter(X, y, color='navy', s=30, marker='o')
colors = ['teal', 'yellowgreen', 'gold', 'orange', 'purple']

for count, degree in enumerate([ 2, 3, 4, 5]):
    model = make_pipeline(PolynomialFeatures(degree), Ridge())
    model.fit(X, y)
    cvs = cross_val_score(model,X,y)
    
    
    y_pred = model.predict(X_plot)
    plt.plot(x_plot, y_pred, color=colors[count], label = 'degree = %i' % degree)
    
plt.legend(loc = 'upper left')
plt.show()




# cv is the number of folds for which is fitted, verbose means how much will be printed
poly_grid = GridSearchCV(estimator = PolynomialFeatures(), param_grid = parameters, cv=5, scoring='accuracy', verbose=3)

poly_grid.fit(X,y) # error: 'PolynomialFeatures' object has no attribute 'predict'

#print((cross_val_score(poly_grid.fit(X,y)))
#y_pred = poly_grid.predict(X_test) 
#acc = poly_grid.score(X_test,y_test)

#acc = accuracy_score(y_test, y_pred, *, normalize=True, sample_weight=None)


res = poly_grid.cv_results_
mean_test_score = res['mean_test_score']

# find best model 
model_poly = poly_grid.best_estimator_
model_poly.score(X_test, y_test) # performance of independent test set
y_pred = model_poly.fit(X, y).predict(X)<jupyter_output><empty_output><jupyter_text>Unfortunately GridSearchCV gave the error that the 'PolynomialFeatures' object has no attribute 'predict'. We could not make this work so we were not possible to do the cross-validation. However we can look at the graph and see that the polynomials with a 2nd degree is not a good option for this dataset. It has a relatively high bias compared to the polynomial of a 3rd or higher degree. The variance is also relatively high in the lower degree polynomials in this case. This would make you believe that the higher the degree, the better the model will perform. This is not the case because in this case overfitting has not yet occurred, when this does happen the scores will go back to being relatively high for the variance as well as the bias.  ROC curve analysisA common method to evaluate binary classifiers is the receiver operating characteristic (ROC) curve. Similar to the week one practicals, implement a $k$-NN classifier on the breast cancer dataset, however, his time use the $k$-NN pipeline from the preliminary. Train the model for different values of $k$ and evaluate their respective performance with an ROC curve, use the `sklearn.metrics.roc_curve` function.<jupyter_code>import matplotlib.pyplot as plt
from sklearn.datasets import load_breast_cancer
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_curve
from sklearn.metrics import auc
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

X, y = load_breast_cancer(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=40)

scaler = StandardScaler()
scaler.fit(X_train)
X_train = scaler.transform(X_train)
X_test = scaler.transform(X_test)

# kNN for k = 1 until k = 10
fig = plt.figure(figsize = (13,20)) #subplots(5,2, figsize = (14,20))
for k in range(1,11):
    knn = KNeighborsClassifier(n_neighbors=k)
    model = Pipeline([
                     ("scaler", scaler),
                     ("knn", knn)
                    ])

    # train the model using the training dataset
    model.fit(X_train, y_train)
    prediction = model.predict(X_test)

    y_scores = model.predict_proba(X_test)
    fpr, tpr, threshold = roc_curve(y_test, y_scores[:, 1])
    roc_auc = auc(fpr, tpr)

    ax = fig.add_subplot(5,2,k)
    ax.plot(fpr, tpr, 'b', label = 'AUC = %0.3f' % roc_auc)
    ax.legend(loc = 'lower right')
    ax.plot([0, 1], [0, 1],'r--')
    ax.set_xlim([0, 1])
    ax.set_ylim([0, 1])
    ax.set_ylabel('True Positive Rate')
    ax.set_xlabel('False Positive Rate')
    ax.set_title('k = %i' % k)

plt.suptitle("ROC Curves with Pipeline",fontsize=14)
plt.subplots_adjust(hspace=0.35)
    
from sklearn.datasets import load_breast_cancer
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_curve
from sklearn.metrics import auc
import matplotlib.pyplot as plt

# new split so not exactly the same data as above
X, y = load_breast_cancer(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X, y, stratify=y, random_state=40)

fig = plt.figure(figsize = (13,20)) 
for k in range(1,11):
    knn = KNeighborsClassifier(n_neighbors = k)
    knn.fit(X_train,y_train)

    y_scores = knn.predict_proba(X_test)
    fpr, tpr, threshold = roc_curve(y_test, y_scores[:, 1])
    roc_auc = auc(fpr, tpr)

    ax = fig.add_subplot(5,2,k)
    ax.plot(fpr, tpr, 'b', label = 'AUC = %0.3f' % roc_auc)
    ax.legend(loc = 'lower right')
    ax.plot([0, 1], [0, 1],'r--')
    ax.set_xlim([0, 1])
    ax.set_ylim([0, 1])
    ax.set_ylabel('True Positive Rate')
    ax.set_xlabel('False Positive Rate')
    ax.set_title('k = %i' % k)

plt.suptitle("ROC Curves without Pipeline",fontsize=14)
plt.subplots_adjust(hspace=0.35)<jupyter_output><empty_output><jupyter_text>$F_1$ score and Dice similarity coefficientThe Dice similarity coefficient is a very popular evaluation measure for image segmentation applications. Assuming that $A$ is the ground truth segmentation of an object represented as a binary image, and $B$ is the binary output of an image segmentation method, the Dice similarity coefficient is computed as:$\text{Dice}(A,B) = \frac{2|A\cap B|}{|A| + |B|}$where $|\cdot|$ represents the cardinality of the objects (e.g. $|A|$ is the number of non-zero pixels in the ground truth segmentation).For example, the Dice similarity can be computed in the following way:<jupyter_code># generate some test objecys
A = np.zeros((32, 32))
A[10:-10, 10:-10] = 1
B = np.zeros((32, 32))
B[5:-15, 5:-15] = 1

dice = 2*np.sum(A*B)/(np.sum(A)+np.sum(B))

# display the results
plt.plot()
plt.imshow(A)
plt.imshow(B, alpha=0.7)

print(dice)<jupyter_output>0.3402777777777778<jupyter_text>Show that the $F_1$ score, which is the harmonic mean of precision and recall, is equivalent to the Dice similarity coefficient<jupyter_code>from sklearn.metrics import f1_score

A = A.flatten()
B = B.flatten()

y_true = A
y_pred = B

f = f1_score(y_true, y_pred)
print (f)

#basic functions:
#precision = tp / (tp + fp)
#recall = tp / (tp + fn)
#f = 2 * ((precision * recall) / (precision + recall))<jupyter_output>0.3402777777777778