<jupyter_start><jupyter_text>Lab 3: Cross validation and feature encoding. The data you will be using in this lab comes from a kaggle in class competition. The goal is to predict the number of shares an article will get on social media, from the article's topic, length, day of publication, and many other features.You are given labels, that is, number of shares, for 5000 of these articles. We will explore how to set up cross validation and do some feature encoding.<jupyter_code>import numpy as np
import pandas as pd
%pylab inline<jupyter_output><empty_output><jupyter_text>1. Model Selection: setting up a cross validationCross-validation is a good way to perform model selection empirically while avoiding overfitting. This procedure can be split into the following two steps: * the dataset is randomly split into K folds * the model is run K times, each run using K-1 folds as the training set and evaluating the performance on the remaining fold which is the test set. Prediction performance are averaged over all folds. When the model contains parameters that need to be tuned, the CV scheme is repeated for all considered values of the hyperparameters, and those leading to the best prediction performance averaged on all folds are retained.Depending on the size of the dataset, 5 or 10 folds are usualy considered. __Question:__ In a K-fold cross-validation, how many times does each sample appear in a test set? In a training set?   __Answer:__ **Question:** Implement a function which splits the _indices_ of the training data in K folds.<jupyter_code>def make_Kfolds(n_instance, n_folds):
    """
    set up a K-fold croos-validation.
    
    Parameters:
    -----------
    n_instances: int
        the number of instances in the dataset.
    n_folds: int
        the number of folds of the cross-validation scheme
        
    Outputs:
    --------
    fold_list: list
        list of folds, a fold is a tuple of 2 lists, 
            the first one containing the indices of instances of the training set,
            the second one containing the indices of instances of the test set
    """
    # Create a list of the n_instance indices [0, 1, ..., n_instance-1]
    list_indices = # TODO 
    # Shuffle the list with np.random.shuffle
    # TODO
    
    # Compute the number of instances per fold (i.e. in each test set)
    n_instance_per_fold = # TODO
    print(n_instance_per_fold)
    
    # For each of the first K-1 folds, create the list of train set and test set indices
    fold_list = []
    for ind_fold in range(n_folds-1):
        test_list = # TODO
        train_list = # TODO
        # TODO add the (train_list, test_list) tuple to fold_list

    # Process the last fold separately
    test_list = # TODO 
    train_list = # TODO 
    # TODO add the (train_list, test_list) tuple to fold_list    
    
    return fold_list
# Check whether your function does what is expected
perso_folds = make_Kfolds(1001, 5)

for ix, (tr, te) in enumerate(perso_folds):
    print("Fold %d" % ix)
    print("\t %d training points" % len(tr))
    print("\t %d test points" % len(te))
    if len(np.intersect1d(tr, te))>0:
        print('some instances are both in your training and test sets')<jupyter_output><empty_output><jupyter_text>In practice, when using scikit-learn, you will not implement your cross-validation yourself, but rather rely on the  library's functionalities for setting up cross-validation schemes. [Here](http://scikit-learn.org/stable/modules/classes.htmlmodule-sklearn.model_selection) is the list of available tools in the scikit-learn library.We list here one of the most important ones:* [K-fold](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.KFold.htmlsklearn.model_selection.KFold): Provides train/test indices to split data in train/test sets by dataset into k consecutive folds (without shuffling by default). * [stratified K-fold](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.StratifiedKFold.htmlsklearn.model_selection.StratifiedKFold) (to be used in case of classification): this cross-validation object is a variation of KFold that returns stratified folds. The folds are made by preserving the percentage of samples for each class.We will now explore the stratified K-fold on randomly generated data.<jupyter_code># Generate random data
n_instances, n_features = 1000, 7
# Design matrix
X = np.random.random((n_instances, n_features))
# Classification labels
y = np.where(np.random.random(n_instances) >=0.5, 1, 0)
print(y)<jupyter_output><empty_output><jupyter_text>**Question:** Using scikit-learn, set up a stratified 10-fold cross-validation for the above data.<jupyter_code>from sklearn import model_selection
# Initialize a StratifiedKFold object 
skf = # TODO
# Split the data using skf
sk_folds = # TODO
# This is one way to access the training and test points
for ix, (tr, te) in enumerate(sk_folds):
    print("Fold %d" % ix)
    print("\t %d training points" % len(tr))
    print("\t %d test points" % len(te))<jupyter_output><empty_output><jupyter_text>__Important note:__ `sk_folds` is a [_generator_](https://wiki.python.org/moin/Generators), meaning that once you are done looping through it, it will be empty. In practice it avoids storing all the indices (if you were doing 10-fold cross-validation on a million sample, you would have $10^7$ values to store). **Question:** Create a cross-validation function that takes a design matrix, label array, scikit-learn classifier, and scikit-learn cross_validation object and returns the corresponding list of cross-validated predictions. The function contains a loop that goes through all folds and for each fold:* trains a model on the training data* uses this model to make predictions on the test data. In this fashion you should be able to form *a single vector of predictions* `y_prob_cv` (as each point from the data appears once as a test point in the cross-validation).Make sure that you are returning the predictions in the correct order!Check the documentation of fit(X, y) and predict_proba(X) in [sklearn.naive_bayes.GaussianNB](http://scikit-learn.org/stable/modules/generated/sklearn.naive_bayes.GaussianNB.html). Every classifier implemented in scikit-learn has a fit(X,y) and a predict_proba(X) methods. Note that the predict_proba methods returns a 2 dimentional array, you must find a way to only keep the probability to belong to the positive class.<jupyter_code>def cross_validate(design_matrix, labels, classifier, cv_folds):
    """ Perform a cross-validation and returns the predictions.
    
    Parameters:
    -----------
    design_matrix: (n_samples, n_features) np.array
        Design matrix for the experiment.
    labels: (n_samples, ) np.array
        Vector of labels.
    classifier:  sklearn classifier object
        Classifier instance; must have the following methods:
        - fit(X, y) to train the classifier on the data X, y
        - predict_proba(X) to apply the trained classifier to the data X and return probability estimates 
    cv_folds: sklearn cross-validation object
        Cross-validation iterator.
        
    Return:
    -------
    pred: (n_samples, ) np.array
        Vectors of predictions (same order as labels).
    """
    pred = np.zeros(labels.shape)
    for tr, te in cv_folds:
        # TODO
    return pred
# To check whether your function runs properly, you can use the following

# import Gaussian Naive Bayes
from sklearn.naive_bayes import GaussianNB
from sklearn import metrics

# create a GNB classifier
gnb = GaussianNB()

# run your cross_validate function
y_prob_cv = cross_validate(X, y, gnb, sk_folds)

# check y and y_prob_cv have the same length (the number of instance)
print(len(y_prob_cv))
print(len(y))

# check the accuracy of your prediction (it should be close to 0.5 as we're considering random matrices). 
print(metrics.accuracy_score(y, np.where(y_prob_cv>=0.5, 1, 0)))<jupyter_output><empty_output><jupyter_text>** Extensions *** **Leave-one-out cross-validation: ** in this case, the number of folds is the number of available points in the dataset. To say it differently, the model is trained K times on K-1 points, and tested on the left out point. The LOO CV scheme is particularly convenient when the number of samples is very small. When the number of samples is large, it becomes computationally burdensome; moreover the cross-validated error tends to have a very large variance which makes it hard to interpret.* **Nested-cross-validation: ** The goal of the cross validation scheme is to assess the performance of the model on _new_ data which were not used to train or optimize the model. From that perspective, the CV scheme is not rigorous when optimizing hyperparameters. Indeed, the test data are both used to assess the performance and choosing the set of parameters which led to that best performance. To avoid selecting a possibly over-fitted set of parameters, we also used the so-called nested cross validation (_Nested CV_) scheme which consists in a cross validation (_inner-CV_) nested in a other cross validation (_outer-CV_). At each step of the _outer-CV_, the optimal parameters are found via the _inner-CV_ on the train set of the _outer-CV_, and the performance is assessed on the remaining test fold of the _outer-CV_ Therefore, in _Nested CV_, parameter optimization and performance assessment are performed on different _unseen_ data.  2. Data loading and visualizationDownload the data from the competition link and unzip it.<jupyter_code># we display the description of the features
!cat data/kaggle_data/features.txt
feature_data = pd.read_csv('data/kaggle_data/features.txt', header=None, sep="  ",
                           names=['feature_names', 'feature_description'])
feature_data<jupyter_output><empty_output><jupyter_text>** Now, let's load and look at the distribution of number of shares (output). **<jupyter_code>target_data = pd.read_csv('data/kaggle_data/train-targets.csv', sep=",")
target_data.head(5)
y_tr = target_data['Prediction'].values
plt.hist(y_tr,bins=2000)
plt.xlim((0,10000))<jupyter_output><empty_output><jupyter_text>** Now, let's load and visualize the features. **<jupyter_code>list_feature_names = list(feature_data['feature_names'])

train_data = pd.read_csv('data/kaggle_data/train.csv', header=None, sep=" ", names=list_feature_names)
train_data.head(5)
test_data = pd.read_csv('data/kaggle_data/test-val.csv', header=None, sep=" ", names=list_feature_names)
test_data.head(5)<jupyter_output><empty_output><jupyter_text>Let us use visualizations to explore the relationships between pairs of features, and between a feature an the output:<jupyter_code>from pandas import plotting
plotting.scatter_matrix(train_data.get(["nb_words_content", "pp_uniq_words"]), alpha=0.2,
               figsize=(6, 6), diagonal='kde')
import seaborn as sns
sns.set_style('whitegrid')

sns.jointplot("nb_words_content", "pp_uniq_words", data = train_data, 
              kind='reg', height=6, space=0)
sns.jointplot(train_data["pp_neg_words"], target_data['Prediction'], 
              kind='reg', ylim=5000, height=6, space=0)<jupyter_output><empty_output><jupyter_text>** Question: ** Change the features you display to explore relationships. What conclusions are you drawing from this exploratory analysis? Are you going to keep all the features in your predictors? __Answer:__  2. Data transformation  2.1 feature engineeringThis notion includes all kinds of manual modification and creation of features. All are of course problem dependant.* __Encoding categorical features:__ if a K-categorical feature is not ordered (categorie 1 is as far to categorie 2 as to categorie 3 etc), then it must not be encoded by a single integer specifying the categorie. We can encode such feature by creating K-1 binary features encoding the belonging to k-th category. (see [link](http://scikit-learn.org/stable/modules/preprocessing.htmlencoding-categorical-features))* __Feature binarization:__ some continuous features can gain predictive power when binarized. For exemple, in some prediction tasks, weekdays could be split into $working\ days$ and $not\ working\ days$. (see [link](http://scikit-learn.org/stable/modules/preprocessing.htmlbinarization))* __Imputation of missing values:__ there are multiple strategies to input missing values when required (see [link](http://scikit-learn.org/stable/modules/preprocessing.htmlimputation-of-missing-values)).* __Dealing with time features or other periodic features:__ when considering the hour of the day as a feature, we can't encode it by the an integer between 1 and 24 as midnigth is as close to 11pm to 1am. An easy strategy to encode periodic features is to apply this transformation $x \mapsto \sin(\frac{2\pi x}{T})$ (T is the period). In the case of the hour of the day, it is   $x \mapsto \sin(\frac{2\pi x}{24})$. * __Generating new features:__ you might want to combine the existing features into new ones that seem informative to you. It can be useful for exemple, notably when working with linear models, to generate polynomial features from the original ones. You can also use external data to transform your features; for instance, if one feature is a date, adding a feature that qualifies whether the day is a working day, a weekday or a holiday can be useful.  * ...In many practical cases, feature engineering is the key to obtaining a huge improvement in performance. ** Question: ** How do you want to engineer the features of the challenge (first, you can start encoding the categorical features)? Keep thinking of this question all along the challenge. Let us encode weekdays as a categorical feature rather than a periodic one. Remember this transformation later in the challenge: does it help your performance?<jupyter_code># Get the weekday data and encode it using a dummy categorical encoding
weekday_data = pd.get_dummies(train_data['weekday'], prefix='weekday', drop_first=True)

# Get the rest of the data
other_data = train_data.drop(['weekday'], axis=1)

# Create a new data set by concatenation of the new weekday data and the old rest of the data
training_data = pd.concat([weekday_data, other_data], axis=1)

# Print the created training data.
training_data.head(5)<jupyter_output><empty_output><jupyter_text>__Question:__ Repeat the process for the other categorical variable(s) in your data. Do not forget to apply your transformation to the test dataset as well!<jupyter_code># TODO<jupyter_output><empty_output><jupyter_text>2.2 Preprocessing data: standardizationYou might want to consider standardizing your data as seen in Lab 02.  2.3 Unsupervised projectionIf your number of features is high, it may be useful to reduce it with an unsupervised step prior to supervised steps. We have already worked on a widly used dimentionality reduction method in `Lab 1`, the Principal Component Analysis. We will discuss in `Lab 5` the combinaison of dimentionality reduction and a predictor.  2.4 Feature selectionSee [scikit-learn's documentation](http://scikit-learn.org/stable/modules/feature_selection.html).It may be useful to select a restricted number of important features to increase their predictive power. When the number of feature is particularly bigger than the number od instance, this issue of major importance.Multiple strategies can be considered depending on the problem such like:* considering the most varying features, condering the most correlated features to the output etc* using feed forward selection procedure: recursively adding features one by one by incresing improvement of performance* using embbeded feature selection like lasso or ElasticNet (see lab 5)* computing feature importance (via bagging procedure like [randomized lasso](https://stat.ethz.ch/~nicolai/stability.pdf) or bagging trees (see lab 5) for exemple) and thresholding the feature.* ...  3. Model evaluation and model selection  3.1 Our first classifier: Gaussian Naive BayesDocumentation: http://scikit-learn.org/stable/modules/generated/sklearn.naive_bayes.GaussianNB.html In order to start thinking about model evaluation and model selection, we will convert the regression problem of the KaggleInClass challenge into a classification task in order to to work with the first classifier we studied in class: the Gaussian Naive Bayes.Our goal here is to try to classify points between astonishingly and not astonishingly shared articles. Based on the distribution of the number of shares, the separation can be set at 1800 shares.<jupyter_code># Transform output into a classification task.
y_clf = np.where(y_tr >= 1800, 1, 0)
print(np.where(y_clf==0)[0].shape)
print(np.where(y_clf==1)[0].shape)
X_clf = training_data.values
# import Gaussian Naive Bayes
from sklearn.naive_bayes import GaussianNB
# create a Gaussian Naive Bayes classifier i.e. an instance of GaussianNB
gnb = GaussianNB()
# fit the classifier to the data
gnb.fit(X_clf, y_clf)
# predict on the same data
y_pred = gnb.predict(X_clf)
# compute the number of mislabeled articles
print("Number of mislabeled points out of a total %d points : %d" % \
      (X_clf.shape[0], (y_clf != y_pred).sum()))<jupyter_output><empty_output><jupyter_text>Note than all predictors implemented in the sklearn library are trained and applied via the `fit` and `predict` (or `predict_proba`) methods.**Question:** What are the parameters of the model we have trained? How many of them are they? How can you access them?<jupyter_code># Hint
gnb.__dict__<jupyter_output><empty_output><jupyter_text>3.2 Model Evaluation You must have a look at http://scikit-learn.org/stable/modules/model_evaluation.html which shows and details a list of metrics for evaluating regression or classification models. In the case of regression, the most commonly used metrics are :* `mean squared errors`* `mean absolute errors` which gives less importance to errors of very bad prediction and more importance to errors of good predictions as the following plot shows than `mean squared errors`* `R2` (coefficient of determination) which provides a measure of how well future samples are likely to be predicted by the model.<jupyter_code>x = np.arange(-2,2,0.01)
plt.plot(x,x*x, 'blue', label='$x^2$')
plt.plot(x,abs(x), 'orange', label='|x|')
plt.legend(loc="upper center")
plt.show()<jupyter_output><empty_output><jupyter_text>In the case of classification, lots of metrics are used depending on the considered problem:* `accuracy` is a default performance measure computing the proportion of missclassified tested instances* `sensitivity` or 'true positive rate' is the proportion of well classified positive samples* `specificity` or 'true negative rate' is the proportion of well classified negative samples* `precision` is the ability of the classifier not to label as positive a sample that is negative. Like in the case of cancer, we really want to avoid diagnose a cancer to somebody who does not have one.* `recall` is the ability of the classifier to find all the positive samples.    * `the area under the precision-recall curve`* `the area under the Receiver operating characteristic (ROC) curve`   **Question** Use the sklearn library to compute the accuracy score of the above prediction.<jupyter_code>from sklearn import metrics
# Score the predictions
print("Accuracy: %.3f" % # TODO
         )<jupyter_output><empty_output><jupyter_text>Building an ROC curve requires to use the probability estimates for the test data points *before* they are thresholded.<jupyter_code># Predict probability estimates instead of 0/1 class labels
y_prob = gnb.predict_proba(X_clf)
print(y_prob.shape)<jupyter_output><empty_output><jupyter_text>**Question:** `y_prob` returns two values for each data point because it returns one probability estimate per class for each data point. The order in which the classes appear are given by `gnb.classes_ `. How do you get the 1-dimensional array that only contains the estimated probability for each point to belong to the positive class?<jupyter_code>pos_index = list(gnb.classes_).index(1)

# ROC curve
fpr, tpr, thresholds = metrics.roc_curve(y_clf, y_prob[:, pos_index], pos_label=1)

# Area under the ROC curve
auc = metrics.auc(fpr, tpr)

# Plot the ROC curve
plt.plot(fpr, tpr, '-', color='orange', label='AUC = %0.3f' % auc)

plt.xlabel('False Positive Rate', fontsize=16)
plt.ylabel('True Positive Rate', fontsize=16)
plt.title('ROC curve: Gaussian Naive Bayes', fontsize=16)
plt.legend(loc="lower right")<jupyter_output><empty_output><jupyter_text>**Question:** What is it problematic to have evaluated our classifier on the training data?  __Answer:__  3.3 Model Selection: cross-validation We will now use the function `make_Kfolds` you have implemented in the first section to evaluate the accuracy of your model via a 5-fold cross-validation scheme. We will compare the results you obtained with those you get with scikit-learn's implementation of the cross-validation scheme.<jupyter_code># Set up a cross-validation with make_Kfolds
perso_folds = # TODO

# Set up a cross-validation with sklearn
sk_folds = # TODO
# Assess performance using the cross_validate function you have implemented
# On perso_folds
gnb = GaussianNB()
y_prob_cv_perso = # TODO use cross_validate and perso_folds
print("Your own cv-scheme: Accuracy: %.3f" % #TODO
     )

# On sk_folds
gnb = GaussianNB()
y_prob_cv_sk = # TODO use cross_validate and perso_folds
print("Sklearn cv-scheme: Accuracy: %.3f" % # TODO
     )<jupyter_output><empty_output><jupyter_text>We will now plot the ROC curve corresponding to your predictions.<jupyter_code># Compute the ROC curve corresponding to the y_prob_cv_perso predictions
fpr, tpr, thresholds = metrics.roc_curve(y_clf, y_prob_cv_perso, pos_label=1)

# Area under the ROC curve
auc = metrics.auc(fpr, tpr)
print("Your own cv-scheme: AUROC: %.3f" % auc)

# Plot the ROC curve
plt.plot(fpr, tpr, '-', color='orange', label='AUC = %0.3f' % auc)

# TODO: plot in blue the ROC curve corresponding to the y_prob_cv_sk predictions

plt.xlabel('False Positive Rate', fontsize=16)
plt.ylabel('True Positive Rate', fontsize=16)
plt.title('ROC curve: Gaussian Naive Bayes', fontsize=16)
plt.legend(loc="lower right")<jupyter_output><empty_output><jupyter_text>__Question:__ The `sklearn.cross_validation` module provides some utilities to make cross-validated predictions. Compare the results you obtained to what they return.Documentation: [cross_val_score](http://scikit-learn.org/stable/modules/generated/sklearn.model_selection.cross_val_score.html)<jupyter_code>gnb = GaussianNB()
skf = model_selection.StratifiedKFold(5, shuffle=True, random_state=91)

# Use model_selection.cross_val_score to compute the average cross-validated roc_auc score 
# of gnb on (X_clf, y_clf), using the skf iterator.
cv_aucs = # TODO

print(np.mean(cv_aucs))

# Note that averaging the AUCs obtained over 10 folds is not the same as 
# globally computing the AUC for the predictions made within the cross-validation loop.<jupyter_output><empty_output><jupyter_text>__Question__ Compare scikit-learn's implementation of the cross-validation with yours.<jupyter_code>gnb = GaussianNB()
skf = model_selection.StratifiedKFold(5, shuffle=True, random_state=91)

# Compute the cross-validation accuracy using model_selection.cross_val_predict
y_pred_sk = model_selection.cross_val_predict(# TODO
     )
print("Cross-validated accuracy: %.3f" %  metrics.accuracy_score(# TODO
     ))

# Compute the cross-validation accuracy using your own cross_validate function
y_prob_cv = cross_validate(# TODO
    )
# Transform y_prob_cv into a vector of binary predictions
y_pred_perso = # TODO 
print("Cross-validated accuracy: %.3f" %  metrics.accuracy_score(# TODO
     ))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0772985482852935, 'fraction_numerical': 0.007195455501788344, 'mean_word_length': 4.849372384937238, 'pattern_counts': {'":': 0, '<': 85, '<?xml version=': 0, '>': 89, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '254616', 'license': 'MIT', 'n_tokens_mistral': 6857, 'n_tokens_neox': 6545, 'n_words': 2847, 'path': '03-ProjectIntro.ipynb', 'repo_name': 'NelleV/2022-mines-HPC-AI-TD', 'size': '23765'}