<jupyter_start><jupyter_text>Selecting the Best ModelNow that we know the basics of validation and cross-validation, let's demonstrate with how it works with model selection and selection of hyperparameters.These issues are some of the most important aspects of the practice of machine learning.Of core importance is the following question: *if our estimator is underperforming, how should we move forward?*There are several possible answers:- Use a more complicated/more flexible model- Use a less complicated/less flexible model- Gather more training samples- Gather more data to add features to each sampleThe answer to this question is often counter-intuitive.In particular, sometimes using a more complicated model will give worse results, and adding more training samples may not improve your results!The ability to determine what steps will improve your model is what separates the successful machine learning practitioners from the unsuccessful.  Validation curves in Scikit-LearnLet's look at an example using cross-validation to compute the validation curve for a class of models to help us find the goldilocks zone between bias and variance.Here we will use a *polynomial regression* model: this is a generalized linear model in which the degree of the polynomial is a tunable parameter.For example, a degree-1 polynomial fits a straight line to the data; for model parameters $a$ and $b$:$$y = ax + b$$A degree-3 polynomial fits a cubic curve to the data; for model parameters $a, b, c, d$:$$y = ax^3 + bx^2 + cx + d$$We can generalize this to any number of polynomial features.In Scikit-Learn, we can implement this with a simple linear regression combined with the polynomial preprocessor.We will use a *pipeline* to string these operations together:<jupyter_code>from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline

def PolynomialRegression(degree=2, **kwargs):
    return make_pipeline(PolynomialFeatures(degree),
                         LinearRegression(**kwargs))<jupyter_output><empty_output><jupyter_text>Now let's create some data to which we will fit our model:<jupyter_code>import numpy as np

def make_data(N, err=1.0, rseed=1):
    # randomly sample the data
    rng = np.random.RandomState(rseed)
    X = rng.rand(N, 1) ** 2
    y = 10 - 1. / (X.ravel() + 0.1)
    if err > 0:
        y += err * rng.randn(N)
    return X, y

X, y = make_data(40)<jupyter_output><empty_output><jupyter_text>We can now visualize our data, along with polynomial fits of several degrees:<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import seaborn; seaborn.set()  # plot formatting

X_test = np.linspace(-0.1, 1.1, 500)[:, None]

plt.scatter(X.ravel(), y, color='black')
axis = plt.axis()
for degree in [1, 3, 5]:
    y_test = PolynomialRegression(degree).fit(X, y).predict(X_test)
    plt.plot(X_test.ravel(), y_test, label='degree={0}'.format(degree))
plt.xlim(-0.1, 1.0)
plt.ylim(-2, 12)
plt.legend(loc='best');<jupyter_output><empty_output><jupyter_text>The knob controlling model complexity in this case is the degree of the polynomial, which can be any non-negative integer.A useful question to answer is this: what degree of polynomial provides a suitable trade-off between bias (under-fitting) and variance (over-fitting)?We can make progress in this by visualizing the validation curve for this particular data and model; this can be done straightforwardly using the ``validation_curve`` convenience routine provided by Scikit-Learn.Given a model, data, parameter name, and a range to explore, this function will automatically compute both the training score and validation score across the range:<jupyter_code>from sklearn.model_selection import validation_curve
degree = np.arange(0, 21)
train_score, val_score = validation_curve(PolynomialRegression(), X, y,
                                          param_name='polynomialfeatures__degree',
                                          param_range=degree, cv=7,
                                          scoring='explained_variance')

plt.plot(degree, np.median(train_score, axis=1), color='blue', label='training score')
plt.plot(degree, np.median(val_score, axis=1), color='red', label='validation score')
plt.legend(loc='best')
plt.ylim(0, 1)
plt.xticks(np.arange(0, 20, 2))
plt.xlabel('degree')
plt.ylabel('score');<jupyter_output><empty_output><jupyter_text>This shows precisely the qualitative behavior we expect: the training score is everywhere higher than the validation score; the training score is monotonically improving with increased model complexity; and the validation score reaches a maximum before dropping off as the model becomes over-fit.From the validation curve, we can read-off that the optimal trade-off between bias and variance is found for a third-order polynomial; we can compute and display this fit over the original data as follows:<jupyter_code>plt.scatter(X.ravel(), y)
lim = plt.axis()
y_test = PolynomialRegression(3).fit(X, y).predict(X_test)
plt.plot(X_test.ravel(), y_test);
plt.axis(lim);<jupyter_output><empty_output><jupyter_text>Notice that finding this optimal model did not actually require us to compute the training score, but examining the relationship between the training score and validation score can give us useful insight into the performance of the model.  Learning curves in Scikit-LearnScikit-Learn offers a convenient utility for computing such learning curves from your models; here we will compute a learning curve for our original dataset with a second-order polynomial model and a ninth-order polynomial:<jupyter_code>from sklearn.model_selection import learning_curve

fig, ax = plt.subplots(1, 2, figsize=(16, 6))
fig.subplots_adjust(left=0.0625, right=0.95, wspace=0.1)

for i, degree in enumerate([2, 9]):
    N, train_lc, val_lc = learning_curve(PolynomialRegression(degree),
                                         X, y, cv=7,
                                         train_sizes=np.linspace(0.3, 1, 25))

    ax[i].plot(N, np.mean(train_lc, 1), color='blue', label='training score')
    ax[i].plot(N, np.mean(val_lc, 1), color='red', label='validation score')
    ax[i].hlines(np.mean([train_lc[-1], val_lc[-1]]), N[0], N[-1],
                 color='gray', linestyle='dashed')

    ax[i].set_ylim(0, 1)
    ax[i].set_xlim(N[0], N[-1])
    ax[i].set_xlabel('training size')
    ax[i].set_ylabel('score')
    ax[i].set_title('degree = {0}'.format(degree), size=14)
    ax[i].legend(loc='best')<jupyter_output><empty_output><jupyter_text>This is a valuable diagnostic, because it gives us a visual depiction of how our model responds to increasing training data.In particular, when your learning curve has already converged (i.e., when the training and validation curves are already close to each other) *adding more training data will not significantly improve the fit!*This situation is seen in the left panel, with the learning curve for the degree-2 model.The only way to increase the converged score is to use a different (usually more complicated) model.We see this in the right panel: by moving to a much more complicated model, we increase the score of convergence (indicated by the dashed line), but at the expense of higher model variance (indicated by the difference between the training and validation scores).If we were to add even more data points, the learning curve for the more complicated model would eventually converge.Plotting a learning curve for your particular choice of model and dataset can help you to make this type of decision about how to move forward in improving your analysis. estimator.get_params() Validation in Practice: Grid SearchThe preceding discussion is meant to give you some intuition into the trade-off between bias and variance, and its dependence on model complexity and training set size.In practice, models generally have more than one knob to turn, and thus plots of validation and learning curves change from lines to multi-dimensional surfaces.In these cases, such visualizations are difficult and we would rather simply find the particular model that maximizes the validation score.Scikit-Learn provides automated tools to do this in the grid search module.Here is an example of using grid search to find the optimal polynomial model.We will explore a three-dimensional grid of model features; namely the polynomial degree, the flag telling us whether to fit the intercept, and the flag telling us whether to normalize the problem.This can be set up using Scikit-Learn's ``GridSearchCV`` meta-estimator:<jupyter_code># Identify potential parameters to test
PolynomialRegression().get_params()
from sklearn.model_selection import GridSearchCV

param_grid = {'polynomialfeatures__degree': np.arange(21),
              'linearregression__fit_intercept': [True, False],
              'polynomialfeatures__interaction_only': [True, False]}

grid = GridSearchCV(PolynomialRegression(), param_grid, cv=7)<jupyter_output><empty_output><jupyter_text>Notice that like a normal estimator, this has not yet been applied to any data.Calling the ``fit()`` method will fit the model at each grid point, keeping track of the scores along the way:<jupyter_code>grid.fit(X, y);<jupyter_output><empty_output><jupyter_text>Now that this is fit, we can ask for the best parameters as follows:<jupyter_code>grid.best_params_<jupyter_output><empty_output><jupyter_text>Finally, if we wish, we can use the best model and show the fit to our data using code from before:<jupyter_code>model = grid.best_estimator_

plt.scatter(X.ravel(), y)
lim = plt.axis()
y_test = model.fit(X, y).predict(X_test)
plt.plot(X_test.ravel(), y_test);
plt.axis(lim);<jupyter_output><empty_output><jupyter_text>The grid search provides many more options, including the ability to specify a custom scoring function, to parallelize the computations, to do randomized searches, and more.For information, see the examples in [In-Depth: Kernel Density Estimation](05.13-Kernel-Density-Estimation.ipynb) and [Feature Engineering: Working with Images](05.14-Image-Features.ipynb), or refer to Scikit-Learn's [grid search documentation](http://Scikit-Learn.org/stable/modules/grid_search.html).  Model Evaluation - ClassificationFirst we'll set up a classification prediction with our Iris dataset using the Gaussian Naive Bayes model, then we'll use some of the common classification metrics to demonstrate how to use those functions.1. Start with model using a single predictor2. Use all available predictors<jupyter_code>import pandas as pd
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.naive_bayes import GaussianNB

# Read data
df = pd.read_csv('https://datahub.io/machine-learning/iris/r/iris.csv')

# Create new column to make numeric values of the three classes of irises
reclass = {'Iris-setosa': 1, 'Iris-versicolor': 3, 'Iris-virginica': 5}
df['reclass'] = df['class'].map(reclass)

# First split data into feature (X) and target (y) matrix
X_iris = df['sepallength'].values[:, np.newaxis]
#X_iris = df.drop(columns=['class', 'reclass'])
y_iris = df['class']

# Split data into training and testing sets
Xtrain, Xtest, ytrain, ytest = train_test_split(X_iris, y_iris,
                                                random_state=1)


model = GaussianNB()
model.fit(Xtrain, ytrain)
y_model = model.predict(Xtest)
from sklearn.metrics import accuracy_score, precision_score, recall_score, f1_score, confusion_matrix, roc_auc_score

print(f'Accuracy Score: {accuracy_score(ytest, y_model)}')
print(f'Precision Score: {precision_score(ytest, y_model, average="macro")}')
print(f'Recall Score: {recall_score(ytest, y_model, average="macro")}')
print(f'F1 Score: {f1_score(ytest, y_model, average="macro")}')
plt.figure(figsize=(8, 8))
sns.heatmap(confusion_matrix(ytest, y_model).T,
            xticklabels=model.classes_,
            yticklabels=model.classes_,
            annot=True, cmap=plt.cm.OrRd)

plt.show()<jupyter_output><empty_output><jupyter_text>Model Metrics - Regression<jupyter_code>from sklearn.linear_model import LinearRegression

df = pd.read_csv('https://raw.githubusercontent.com/kgoebber/data_151/main/notebooks/petrol_consumption.csv')

X = df[['Petrol_tax', 'Average_income', 'Paved_Highways',
       'Population_Driver_licence(%)']]
y = df['Petrol_Consumption']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1)

regressor = LinearRegression()
regressor.fit(X_train, y_train)

y_pred = regressor.predict(X_test)
from sklearn import metrics

print(f'Explained Variance Score: {metrics.explained_variance_score(y_test, y_pred)}')
print('Mean Absolute Error:', metrics.mean_absolute_error(y_test, y_pred))
print('Mean Squared Error:', metrics.mean_squared_error(y_test, y_pred))
print('Root Mean Squared Error:', np.sqrt(metrics.mean_squared_error(y_test, y_pred)))<jupyter_output><empty_output><jupyter_text>Now apply ridge regression and cross validation to generate a better model.<jupyter_code>from sklearn.linear_model import RidgeCV
cv_regressor = RidgeCV(alphas=[5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1] , cv=5)
cv_regressor.fit(X_train, y_train)
y_pred = cv_regressor.predict(X_test)
print(f'Explained Variance Score: {metrics.explained_variance_score(y_test, y_pred)}')
print('Mean Absolute Error:', metrics.mean_absolute_error(y_test, y_pred))
print('Mean Squared Error:', metrics.mean_squared_error(y_test, y_pred))
print('Root Mean Squared Error:', np.sqrt(metrics.mean_squared_error(y_test, y_pred)))
cv_regressor.alpha_<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09122266255687655, 'fraction_numerical': 0.009540584177308088, 'mean_word_length': 5.118994162550516, 'pattern_counts': {'":': 0, '<': 53, '<?xml version=': 0, '>': 54, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '209617', 'license': 'MIT', 'n_tokens_mistral': 4128, 'n_tokens_neox': 3913, 'n_words': 1438, 'path': 'notebooks/bias_variance_model_evaluation.ipynb', 'repo_name': 'kgoebber/data_151', 'size': '13626'}