In [1]:
import scipy
import numpy as np
import pandas as pd
import plotly.plotly as py
import visplots
from plotly.graph_objs import *
from plotly.offline import download_plotlyjs, init_notebook_mode, iplot
from sklearn import preprocessing, metrics
from sklearn.cross_validation import train_test_split, cross_val_score
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.grid_search import GridSearchCV, RandomizedSearchCV
from scipy.stats.distributions import randint
init_notebook_mode()
print("libraries all imported, ready to go")
The dataset we will be using throughout this workshop is an adapted version of the Wine Quality case study, available from the UCI Machine Learning repository (https://archive.ics.uci.edu/ml/datasets/Wine+Quality). The goal of this case study is to model the wine quality (into "low" or "high" quality) based on physicochemical tests (such as fixed and volatile acidity, citric acid, etc.).
The first thing you will need to do in order to work with the wine quality dataset is to read the contents from the provided wine_quality.csv data file using the read_csv command. You should also try to explore the first few rows of the imported wine DataFrame using the head function from the pandas package (http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.head.html):
In [8]:
# Import the data and explore the first few rows
df = pd.read_csv('data/wine_quality.csv')
#print df
df.head()
#df.tail()
Out[8]:
In order to feed the data into our classification models and sklearn, the imported wine quality DataFrame needs to be converted into a numpy array. For more information on numpy arrays, see http://scipy-lectures.github.io/intro/numpy/array_object.html.
In addition, it is always a good practice to check the dimensionality of the imported data using the shape command prior to constructing any classification model to make sure you have really imported all the data, and imported it in the correct way (e.g. one common mistake is to get the separator wrong and end up with only one column).
In [21]:
# Convert to numpy array and check the dimensionality
array = df.as_matrix()
print len(array)
print len(array[0])
print np.shape(array)
print type(array)
To select elements in an array, you specify their indices with square bracket notation. For a two-dimensional array, the first index indicates the row number and the second index indicates the column number. Try selecting the values of the first and second columns of the first sample in the npArray:
In [28]:
# Print the 1st row and 1st column of npArray
print array[0]
print array[:,0]
print np.shape(array[0])
print np.shape(array[:,0])
In [29]:
# Print the 1st row and 2nd column of npArray
print array[0]
print array[:,1]
print np.shape(array[0])
print np.shape(array[:,1])
To select ranges of elements, we use "index slicing". Index slicing is the technical name for the syntax A[lower:upper], where lower refers to the lower bound index that is included, and upper refers to the upper bound index that is not included. Try selecting the first three samples (rows):
In [51]:
# Print the first 3 rows of npArray
print array[0:3]
print np.shape(array[0:3])
and also the first three samples (rows) of the last column:
In [32]:
# Print the first 3 rows from the last column of npArray
print array[0:3,-1]
Subsequently, we need to split our initial dataset into the data matrix X (independent variable) and the associated class vector y (dependent or target variable). The input features, X, are the variables that you use to predict the outcome. In this data set, there are ten input features stored in columns 1-10 (index 0-9, although the upper bound is not included so the range for indexing is 0:10), all of which have continuous values. The output label, y, holds the information of whether the wine has been rated as high or low quality, and is stored in the final (eleventh) column (index 10). To split the data, we need to assign the columns of the input features and the columns of the output labels to different arrays:
In [65]:
# Split to input matrix X and class vector y
y = array[:,-1]
x = array[:,0:10]
Try printing the size of the input matrix X and class vector y using the "shape" command:
In [66]:
# Print the dimensions of X and y
print np.shape(y)
print np.shape(x)
Visualisation is an integral part of Data Science. Exploratory data analysis (EDA) is the field dealing with the analysis of data sets as a means of summarising their main characteristics, most often using visual methods.
Plotly is an online collaborative data analysis and graphing tool that we will use in order to construct fully interactive graphs. The Plotly API allows you to access all of the library's interactive functionality directly from Python (or other programming languages such as R, JavaScript and MATLAB, among others). Crucially, Plotly has recently been made open-source, which now enables plotting offline without requiring access to their API. Plotly Offline brings interactive Plotly graphs to the offline Jupyter (IPython) Notebook environment.
An important aspect to understand before applying any classification algorithms is how the output labels are distributed. Are they evenly distributed? Imbalances in distribution of labels can often lead to poor classification results for the minority class even if the classification results for the majority class are very good.
In [67]:
# Print the y frequencies
lows = 0
for e in y:
if e == 'low':
lows = lows + 1
highs = len(y) - lows
yFreq = scipy.stats.itemfreq(y)
print 'lows ' + str(lows)
print 'hights ' + str(highs)
print yFreq
In our current dataset, you can see that the y values are categorical (i.e. they can only take one of a discrete set of values) and have a non-numeric representation, "high" vs. "low". This can be problematic for scikit-learn and plotting functions in Python, since they assume numerical values, so we need to map the text categories to numerical representations using LabelEncoder and the fit_transform function from the preprocessing module:
In [68]:
# Convert the categorical to numeric values, and print the y frequencies
le = preprocessing.LabelEncoder()
y = le.fit_transform(y)
print y
Visualising the data in some way is a good way to get a feel for how the data is distributed. As a simple example, try plotting the frequencies of the class labels (held in yFreq), 1 and 0, and see how they are distributed using a barplot from Plotly:
In [69]:
# Display the y frequencies in a barplot with Plotly
py.sign_in('santiaago', 'dlekv8an96')
data = [
graph_objs.Bar(
x=['low', 'high'],
y=[lows, highs]
)
]
plot_url = py.plot(data, filename='basic-bar')
More examples on Plotly barplots can be found at https://plot.ly/python/bar-charts/. In addition, a full list of arguments on barplots can be found at https://plot.ly/python/reference/#bar/.
It is usually advisable to scale your data prior to fitting a classification model to avoid attributes with
greater numeric ranges dominating those with smaller numeric ranges. In order to investigate the range and descriptive statistics of our features, we can apply the describe() function from pandas to the original wineQ DataFrame (not the numpy array!). For instance:
In [70]:
# Print the descriptive statistics of the wineQ DataFrame
df.describe()
Out[70]:
Boxplots are a powerful visual aid, commonly used in order to investigate simultaneously the range differences of the input features. Boxplots are a standardised way of displaying the distribution of the data based on the "five number summary" (minimum, first quartile, median, third quartile, and maximum). For example, try and plot the features of the raw matrix X using the script for the boxplots:
In [83]:
# Create a boxplot of the raw data
import plotly.graph_objs as go
traces = []
t0 = go.Box(y=x[:,0])
t1 = go.Box(y=x[:,1])
t2 = go.Box(y=x[:,2])
t3 = go.Box(y=x[:,3])
t4 = go.Box(y=x[:,4])
t5 = go.Box(y=x[:,5])
t6 = go.Box(y=x[:,6])
t7 = go.Box(y=x[:,7])
t8 = go.Box(y=x[:,8])
t9 = go.Box(y=x[:,9])
traces = [t1,t2,t3,t4,t5,t6,t7,t8,t9]
data = traces
plot_url = py.plot(data, filename='basic-box-plot')
There are many ways of scaling but one common scaling mechanism is auto-scaling, where for each
column, the values are centred around the mean and divided by their standard deviation. This scaling
mechanism can be applied by calling the scale() function in scikit-learn’s preprocessing module.
In [84]:
# Auto-scale the data
x_scale = preprocessing.scale(x)
Try to re-run the previous plotting script and have a look at the outcome of the boxplot after scaling. Alternatively, if you feel more adventurous, you create a more enhanced version of the boxplot. You can find more online examples at https://plot.ly/python/box-plots/, and also a full list of boxplot arguments at https://plot.ly/python/reference/#box.
In [85]:
# Create a boxplot of the scaled data (simple or enhanced)
traces1 = []
t0 = go.Box(y=x_scale[:,0])
t1 = go.Box(y=x_scale[:,1])
t2 = go.Box(y=x_scale[:,2])
t3 = go.Box(y=x_scale[:,3])
t4 = go.Box(y=x_scale[:,4])
t5 = go.Box(y=x_scale[:,5])
t6 = go.Box(y=x_scale[:,6])
t7 = go.Box(y=x_scale[:,7])
t8 = go.Box(y=x_scale[:,8])
t9 = go.Box(y=x_scale[:,9])
traces1 = [t1,t2,t3,t4,t5,t6,t7,t8,t9]
data = traces1
plot_url = py.plot(data, filename='basic-box-plot')
You can visualise the relationship between two variables (features) using a simple scatter plot. This step can give you a good first indication of the ML model model to apply and its complexity (linear vs. non-linear). At this stage, let’s plot the first two variables against each other:
In [87]:
# Create a scatter plot of the first two features
trace1 = go.Scatter(
x = x[:,0],
y = x[:,1],
mode = 'markers'
)
data = [trace1]
# Plot and embed in ipython notebook!
py.iplot(data, filename='basic-scatter')
Out[87]:
We can also relate associations between features to their y classifications by making the colour of the points dependent on the corresponding y value:
In [ ]:
# Create an enhanced scatter plot of the first two features
Examples of Plotly scatterplots can be found at https://plot.ly/python/line-and-scatter/ (or for a list of arguments refer to https://plot.ly/python/reference/#scatter/).
The scatterplots we have seen so far investigated the relationship between two variables (features). A three-dimensional graph lets you introduce a third axis, typically called the z axis, and can help you understand the relationship between three variables. Plotly's fully interactive functionality allows you to plot, hover, zoom and rotate 3-dimensional scatterplots. For a full list of arguments on 3d plots in Plotly visit https://plot.ly/python/reference/#scatter3d. Other examples on 3D scatterplots using Plotly can be found at https://plot.ly/python/3d-scatter-plots/.
Hint: Investigate the Scatter3d object from Plotly
Axes in 3D Plotly plots work a bit differently than in 2D (axes are bound to a Scene object -- use help(Scene)).
In [ ]:
# Create a 3D scatterplot using the first three features
trace1 = go.Scatter(
x = x[:,0],
y = x[:,1],
mode = 'markers'
)
data = [trace1]
# Plot and embed in ipython notebook!
py.iplot(data, filename='basic-scatter')
A scatterplot matrix shows a grid of scatterplots where each attribute is plotted against all other attributes. For example, try to create a scatterplot matrix of the first four features. You can find further information on how to create and customise subplots with Plotly at https://plot.ly/python/subplots/.
Hints: You may want to use nested loops that iterate through the rows and columns of the grid, and also import and make use of the make_subplots() function from Plotly
In [ ]:
# Create a grid plot of scatterplots using a combination of features
Often, the different features (variables) in X are not completely independent from each other. For example,
fixed acidity is related to volatile acidity. To quickly identify which features are related and to
what extent, it is useful to see how they are correlated. You can do this by creating a correlation matrix
from X using corrcoef() in the numpy module:
In [ ]:
# Calculate the correlation coefficient
To search for linear relationships between features across all pairs of features, you can use a heatmap of correlations (directly from X), which is simply a matrix of subplots whose colours represent the sizes of the correlations:
In [ ]:
# Create a heatmap of the correlation coefficients
Training and testing a classification model on the same dataset is a methodological mistake: a model that would just repeat the labels of the samples that it has just seen would have a perfect score but would fail to predict anything useful on yet-unseen data (poor generalisation). To use different datasets for training and testing, we need to split the wine dataset into two disjoint sets: train and test (Holdout method) using the train_test_split function.
In [178]:
# Split into training and test sets
X = df.as_matrix()[:,0:10]
y = df.as_matrix()[:,-1]
xTrain, xTest, yTrain, yTest = train_test_split(X, y)
#xTrain, xTest = train_test_split(df)
#xTrain = xTrain.as_matrix()
#xTest = xTest.as_matrix()
#yTrain = xTrain[:,-1]
#yTest = xTest[:,-1]
#xTrain = xTrain[:,0:10]
#xTest = xTest[:,0:10]
XTrain and yTrain are the two arrays you use to train your model. XTest and yTest are the two arrays that you use to evaluate your model. By default, scikit-learn splits the data so that 25% of it is used for testing, but you can also specify the proportion of data you want to use for training and testing.
You can check the sizes of the different training and test sets by using the shape attribute:
In [179]:
# Print the dimensionality of the individual splits
print np.shape(xTrain)
print np.shape(xTest)
print np.shape(yTrain)
print np.shape(yTest)
You can also investigate how the class labels are distributed within the yTest vector by using the itemfreq function as previously
In [180]:
# Calculate the frequency of classes in yTest
scipy.stats.itemfreq(yTest)
Out[180]:
We can see that 129 random samples of class 0 (high quality) and 243 random samples of class 1 (low quality) are included in the yTest set.
To build KNN models using scikit-learn, you will be using the KNeighborsClassifier object, which allows you to set the value of K using the n_neighbors parameter (http://scikit-learn.org/stable/modules/generated/sklearn.neighbors.KNeighborsClassifier.html). The optimal choice for the value K is highly data-dependent: in general a larger K suppresses the effects of noise, but makes the classification boundaries less distinct.
For every classification model built with scikit-learn, we will follow four main steps: 1) Building the classification model (using either default, pre-defined or optimised parameters), 2) Training the model with data, 3) Testing the model, and 4) Performance evaluation using various metrics.
We are going to start by trying two pre-defined random values of K and compare their performance. Let us start with a small number of K such as K=3.
In [181]:
# Build a KNN classifier with 3 nearest neighbors
knn3 = KNeighborsClassifier(n_neighbors=3)
knn3.fit(xTrain, yTrain)
yPredK3 = knn3.predict(xTest)
print("Overall Accuracy:", round(metrics.accuracy_score(yTest, yPredK3), 2))
Let us try a larger value of K, for instance K = 99 or another number of your own choice; remember, it is good practice to select an odd number for K in a binary classification problem to avoid ties. Can you generate the KNN model and print the overall performance for a larger K (such as K=99) using as guidance the previous example?
In [182]:
# Build a KNN classifier with 99 nearest neighbors
knn99 = KNeighborsClassifier(n_neighbors=99)
knn99.fit(xTrain, yTrain)
yPredK99 = knn99.predict(xTest)
print("Overall Accuracy:", round(metrics.accuracy_score(yTest, yPredK99), 2))
In a classification task, once you have created your predictive model, you will need to evaluate it. Evaluation functions help you to do this by reporting the performance of the model through four main performance metrics: precision, recall and specificity for the different classes, and overall accuracy. To understand these metrics, it is useful to create a confusion matrix, which records all the true positive, true negative, false positive and false negative values.
We can compute the confusion matrix for our classifier using the confusion_matrix function in the metrics module.
In [183]:
# Get the confusion matrix for your classifier using metrics.confusion_matrix
mat = metrics.confusion_matrix(yTest, yPredK3)
print(mat)
Because performance metrics are such an important step of model evaluation, scikit-learn offers a wrapper around these functions, metrics.classification_report, to facilitate their computation. It also offers the function metrics.accuracy_score that we tried before to compute the overall accuracy.
In [184]:
# Report the metrics using metrics.classification_report
We can visualise the classification boundary created by the KNN classifier using the built-in function visplots.knnDecisionPlot. For easier visualisation, you can (interactively) select to view only the test samples from the plot. Remember though that the decision boundary has been built using the training data!
In [185]:
# Check the arguments of the function
help(visplots.knnDecisionPlot)
# Visualise the boundaries
yTrain = le.fit_transform(yTrain)
yTest = le.fit_transform(yTest)
header = df.columns
visplots.knnDecisionPlot(xTrain, yTrain, xTest, yTest, header, n_neighbors = 3)
visplots.knnDecisionPlot(xTrain, yTrain, xTest, yTest, header, n_neighbors = 99)
For smaller values of K the decision boundaries present many "creases". In this case the models may suffer from instances of overfitting. For larger values of K, we can see that the decision boundaries are less distinct and tend towards linearity. In these cases the boundaries may be too simple and unable to learn thus leading to cases of underfitting.
Under some circumstances, it is better to give more importance ("weight" in computing terms) to nearer neighbors. This can be accomplished through the weights parameter. When weights = 'distance', weights are assigned to the training data points in a way that is proportional to the inverse of the distance from the query point. In other words, nearer neighbors contribute more to the fit.
What if we use weights based on distance? Does it improve the overall performance?
In [186]:
# Build the classifier with two pre-defined parameters (n_neighbors and weights)
knn = KNeighborsClassifier(n_neighbors=3, weights= 'distance')
knn.fit(xTrain, yTrain)
yPredK = knn.predict(xTest)
print(metrics.classification_report(yTest, yPredK))
print("accuracy: ", round(metrics.accuracy_score(yTest, yPredK), 2))
# Visualise the boundaries of a KNN model with weights equal to "distance"
visplots.knnDecisionPlot(xTrain, yTrain, xTest, yTest, header, n_neighbors = 3, weights = 'distance')
Let us estimate the accuracy of the classifier on the wine quality dataset by splitting the data 5 consecutive times (the parameter cv gives the number of samples the data is split into) using the cross_val_score function. For example, try to implement cross-validation for knn3, your KNN model with k=3:
In [187]:
# Implement cross-validation for knn3
knn3scores = cross_val_score(knn3, X, y, cv = 5)
print(knn3scores)
print("Mean of scores KNN3", knn3scores.mean())
knn10 = KNeighborsClassifier(n_neighbors=10)
knn10.fit(xTrain, yTrain)
knn10scores = cross_val_score(knn10, X, y, cv = 5)
print(knn10scores)
print("Mean of scores KNN10", knn10scores.mean())
What happens if we change the number of K?
Rather than trying one-by-one predefined values of K, we can automate this process. The scikit-learn library provides the grid search function GridSearchCV (http://scikit-learn.org/stable/modules/generated/sklearn.grid_search.GridSearchCV.html), which allows us to exhaustively search for the optimum combination of parameters by evaluating models trained with a particular algorithm with all provided parameter combinations. Further details and examples on grid search with scikit-learn can be found at http://scikit-learn.org/stable/modules/grid_search.html
You can use the GridSearchCV function with the validation technique of your choice (in this example, 10-fold cross-validation has been applied) to search for a parametisation of the KNN algorithm that gives a more optimal model:
In [188]:
# Grid search with 10-fold cross-validation using a dictionary of parameters
n_neighbors = np.arange(1, 51, 2)
weights = ['uniform','distance']
params = [{'n_neighbors': n_neighbors, 'weights': weights}]
gridCV = GridSearchCV(KNeighborsClassifier(), params, cv=10)
gridCV.fit(xTrain, yTrain)
bestN = gridCV.best_params_['n_neighbors']
bestW = gridCV.best_params_['weights']
print("Best parameters: n_neighbors=", bestN, "and weight=", bestW)
We can also graphically represent the results of the grid search using a heatmap:
In [189]:
# Visualise the grid search results using a heatmap
scores = np.zeros((len(n_neighbors), len(weights)))
# grid_scores_ contains parameter settings and scores
for score in gridCV.grid_scores_:
ne = score[0]['n_neighbors']
i = np.argmax(n_neighbors == ne)
j = 0 if (score[0]['weights'] == 'uniform') else 1
scores[i,j] = score[1]
# Make a heatmap with the performance
data = [
Heatmap(
x = n_neighbors,
y = weights,
z = scores.T,
colorscale='Jet',
reversescale=True,
colorbar = dict(
title = "Classification Accuracy",
len = 4,
nticks=10
)
)
]
layout = Layout(
xaxis = dict(title = "Number of K nearest neighbors", tickvals = n_neighbors),
yaxis = dict(title = "Weights"),
height= 250
)
fig = dict(data = data, layout = layout)
iplot(fig)
When evaluating the resulting model it is important to do it on held-out samples that were not seen during the grid search process (XTest).
So, we are testing our independent XTest dataset using the optimised model:
In [191]:
# Build the classifier using the optimal parameters detected by grid search
knnBest = KNeighborsClassifier(n_neighbors=bestN, weights=bestW)
knnBest.fit(xTrain, yTrain)
yPredKnnBest = knnBest.predict(xTest)
# Report the test accuracy and performance
print(metrics.classification_report(yTest, yPredKnnBest))
print("Overall Accuracy:", round(metrics.accuracy_score(yTest, yPredKnnBest), 2))
We can graphically represent and investigate how the systematic increase of the number of K neighbors influences the validation, train and test accuracy (and attempt to detect cases of over- or under-fitting).
In [192]:
# Explore the benefit of cross-validated results vs. simple training and test data separation
train_scores = []
test_scores = []
cv_scores = [x[0] for x in scores]
for n in n_neighbors:
knn = KNeighborsClassifier(n_neighbors=n)
knn.fit(xTrain, yTrain)
train_scores.append(metrics.accuracy_score(yTrain, knn.predict(xTrain)))
test_scores.append(metrics.accuracy_score(yTest, knn.predict(xTest)))
In [193]:
# Plot the train, test and validation accuracies
trace0 = Scatter(
x = n_neighbors,
y = train_scores,
mode = "lines+markers",
name = "Training Scores"
)
trace1 = Scatter(
x = n_neighbors,
y = test_scores,
mode = "lines+markers",
name = "Test Scores"
)
trace2 = Scatter(
x = n_neighbors,
y = cv_scores,
mode = "lines+markers",
name = "CV Scores"
)
layout = Layout(
xaxis = dict(title = 'number of neighbors'),
yaxis = dict(title = 'prediction accuracy')
)
fig = Figure(data=[trace0, trace1, trace2], layout=layout)
iplot(fig)
Unlike GridSearchCV, RandomizedSearchCV does not exhaustively try all the parameter settings. Instead, it samples a fixed number of parameter settings based on the distributions you specify (e.g. you might specify that one parameter should be sampled uniformly while another is sampled following a Gaussian distribution). The number of parameter settings that are tried is given by n_iter. If all parameters are presented as a list, sampling without replacement is performed. If at least one parameter is given as a distribution, sampling with replacement is used. You should use continuous distributions for continuous parameters. Further details can be found at http://scikit-learn.org/stable/modules/grid_search.html
In [194]:
# Conduct a randomised search on hyperparameters
parameters = {'n_neighbors': randint(1,200)}
randomCV = RandomizedSearchCV(KNeighborsClassifier(),
param_distributions=parameters, n_iter=20)
randomCV.fit(xTrain, yTrain)
Out[194]:
In [195]:
# Print the optimal n_neighbors detected by randomised search
yRandomCV = randomCV.predict(xTest)
# Report the test accuracy and performance
print(metrics.classification_report(yTest, yRandomCV))
print("Overall Accuracy:", round(metrics.accuracy_score(yTest, yRandomCV), 2))
We can also graphically represent the results of the randomised search using a scatterplot:
In [202]:
# Visualise the randomised search results using a scatterplot
neighbor = [score_tuple[0]['n_neighbors'] for score_tuple in randomCV.grid_scores_]
result = [score_tuple[1] for score_tuple in randomCV.grid_scores_]
print neighbor
print result
trace = go.Scatter(
x = neighbor,
y = result,
mode = 'markers'
)
data = [trace]
# Plot and embed in ipython notebook!
py.iplot(data, filename='basic-scatter')
Out[202]:
In [205]:
# Build the classifier using the optimal parameters detected by randomised search
#print neighbor[33]
knn33 = KNeighborsClassifier(n_neighbors=33)
knn33.fit(xTrain, yTrain)
y33 = knn33.predict(xTest)
print(metrics.classification_report(yTest, y33))
print("Overall Accuracy:", round(metrics.accuracy_score(yTest, y33), 2))
Decision Tree classifiers construct classification models in the form of a tree structure. A decision tree progressively splits the training set into smaller subsets. Each node of the tree represents a subset of the data. Once a new sample is presented to the data, it is classified according to the test condition generated for each node of the tree.
Let us build a simple decision tree with 3 layers. We will first evaluate the accuracy, then plot the decision boundaries just as we did for knn. (See http://scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeClassifier.html for the documentation of the classifier.)
In [208]:
# Build a Decision Tree classifier with 3 layers
dtc = DecisionTreeClassifier(max_depth=3)
dtc.fit(xTrain, yTrain)
predDT = dtc.predict(xTest)
print(metrics.classification_report(yTest, predDT))
print("Overall Accuracy:", round(metrics.accuracy_score(yTest, predDT),2))
visplots.dtDecisionPlot(xTrain, yTrain, xTest, yTest, header, max_depth=3)
The random forests model is an ensemble method since it aggregates a group of decision trees into an ensemble (http://scikit-learn.org/stable/modules/ensemble.html). Ensemble learning involves the combination of several models to solve a single prediction problem. It works by generating multiple classifiers/models which learn and make predictions independently. Those predictions are then combined into a single (mega) prediction that should be as good or better than the prediction made by any one classifer. Unlike single decision trees which are likely to suffer from high Variance or high Bias (depending on how they are tuned) Random Forests use averaging to find a natural balance between the two extremes.
Let us start by building a simple Random Forest model which consists of 100 independently trained decision trees. For further details and examples on how to construct a Random Forest, see http://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html
In [228]:
# Build a Random Forest classifier with 100 decision trees
rf = RandomForestClassifier(n_estimators=100, random_state=0)
rf.fit(xTrain, yTrain)
predRF = rf.predict(xTest)
print(metrics.classification_report(yTest, predRF))
print("Overall Accuracy:", round(metrics.accuracy_score(yTest, predRF),2))
In [229]:
# Visualise the average accuracy
visplots.rfAvgAcc(rfModel=rf, XTest=xTest, yTest=yTest)
Random forests allow you to compute a heuristic for determining how “important” a feature is in predicting a target. This heuristic measures the change in prediction accuracy if you take a given feature and permute (scramble) it across the datapoints in the training set. The more the accuracy drops when the feature is permuted, the more “important” we can conclude the feature is.
We can use the feature_importances_ attribute of the RF classifier to obtain the relative importance of each feature, which we can then visualise using a simple bar plot.
In [215]:
# Display the importance of the features in a barplot
data = [
go.Bar(
x=rf.feature_importances_,
y=header[0:10],
orientation = 'h',
)
]
plot_url = py.plot(data, filename='horizontal-bar')
We can visualise the classification boundary created by the Random Forest using the visplots.rfDecisionPlot function. You can check the arguments passed in this function by using the help command. For easier visualisation, only the test samples have been included in the plot. And remember that the decision boundary has been built using the training data!
In [216]:
# Check the arguments of the function
help(visplots.rfDecisionPlot)
# Visualise the boundaries
visplots.rfDecisionPlot(XTrain, yTrain, XTest, yTest, header)
In [217]:
# View the list of arguments to be optimised
help(RandomForestClassifier())
Create a dictionary of allowed parameter ranges for n_estimators and max_depth (or include more of the parameters you would like to tune) and conduct a grid search with cross validation using the GridSearchCV function as before:
In [219]:
# Conduct a grid search with 10-fold cross-validation using the dictionary of parameters
n_estimators = np.arange(1, 30, 5)
max_depth = np.arange(1, 100, 5)
parameters = [{'n_estimators': n_estimators, 'max_depth': max_depth}]
gridCV = GridSearchCV(RandomForestClassifier(), param_grid=parameters, cv=10)
gridCV.fit(xTrain, yTrain)
# Print the optimal parameters
best_n_estim = gridCV.best_params_['n_estimators']
best_max_depth = gridCV.best_params_['max_depth']
print ("Best parameters: n_estimators=", best_n_estim,", max_depth=", best_max_depth)
Finally, testing our independent XTest dataset using the optimised model:
In [220]:
# Build the classifier using the optimal parameters detected by grid search
clfRDF = RandomForestClassifier(n_estimators=best_n_estim, max_depth=best_max_depth)
clfRDF.fit(XTrain, yTrain)
predRF = clfRDF.predict(XTest)
print(metrics.classification_report(yTest, predRF))
print("Overall Accuracy:", round(metrics.accuracy_score(yTest, predRF),2))
We can also graphically represent the results of the grid search using a heatmap:
In [224]:
# Create a heatmap like the one you made when you applied GridSearchCV to KNN
scores = np.zeros((len(n_estimators), len(max_depth)))
for score in gridCV.grid_scores_:
ne = score[0]['n_estimators']
md = score[0]['max_depth']
i = np.argmax(n_estimators == ne)
j = np.argmax(max_depth == md)
scores[i,j] = score[1]
print scores
data = [
Heatmap(
x = n_neighbors,
y = weights,
z = scores.T,
colorscale='Jet',
reversescale=True,
colorbar = dict(
title = "Create a Plotly heatmap based on the scores of gridCV for the combination of n_estimators and max_depth",
len = 4,
nticks=10
)
)
]
layout = Layout(
xaxis = dict(title = "Number of estimators", tickvals = n_neighbors),
yaxis = dict(title = "Max depth"),
height= 750
)
fig = dict(data = data, layout = layout)
iplot(fig)
The scikit-learn implementation of Random Forests also features the parallel construction of the trees and the parallel computation of the predictions through the n_jobs parameter.
If n_jobs=k then computations are partitioned into k jobs, and run on k cores of the machine.
If n_jobs=-1 then all cores available on the machine are used.
In [226]:
# Change the value of n_jobs and estimate the excution time of fit
import timeit
n_jobs_to_try = np.arange(1,9)
elapsed_times = []
for jobs in n_jobs_to_try:
start_time = timeit.default_timer()
rf = RandomForestClassifier(n_estimators=500, random_state=0, n_jobs=jobs)
rf.fit(XTrain, yTrain)
elapsed = timeit.default_timer() - start_time
elapsed_times.append(elapsed)
In [227]:
# Plot a graph
data = [
Bar(
x=n_jobs_to_try,
y=elapsed_times
)
]
fig = dict(data=data, layout=Layout(yaxis=dict(title='seconds')))
iplot(fig)
In [ ]: